SCHWARZ_PBCGSTAB BiCGStab method with additive Schwarz preconditioner with overlap and coarse mesh
[u,iter,res]=schwarz_pbcgstab(u0, f, param,...
p_unity,...
xy, ww, A, nov, noei, lint,x,wx,y,wy,xx,jacx,yy,jacy,...
Aq1, wwq1, linte,nove,nvle,...
Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc);
Input : u0 = column array, initial datum for pbcgstab
f = column array, r.h.s. of linear system to solve
param = array of parameters
p_unity = unity partition associated to extended elements
xy = 2-indexes array wiht coordinates of 2D LGL mesh
ww = mass matrix (in column array)
nov = local -global map, previously generated by cosnov_2d
noei = number of nodes internal to Omega
lint = list of internal nodes (set in liste.n)
x = npdx LGL nodes in [-1,1],
(produced by calling [x,wx]=xwlgl(npdx))
wx= npdx LGL weigths in [-1,1],
(produced by calling [x,wx]=xwlgl(npdx))
y = npdx LGL nodes in [-1,1],
(produced by calling [y,wy]=xwlgl(npdy))
wy= npdx LGL weigths in [-1,1],
(produced by calling [y,wy]=xwlgl(npdy))
xx = 2-indexes array of size (4,ne)
xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie]
(ne=nex*ney is the global number of spectral elements)
jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
yy = 2-indexes array of size (4,ne):
yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie]
jacy = array (length(jacy)=ne); jacy(ie)= (x_V3_ie-x_V1_ie)/2
Aq1 = cell array with Q1 stiffness matrices (internal/internal nodes)
on extended elements
wwq1 = cell array with Q1 mass matrices (internal nodes)
on extended elements
linte = cell array with list of internal nodes of extended
elements
nove = 2-index array of "extended local" to global map, permuted.
nvle = 2-index array:
nvle (:,1)=number of nodes of the extendend elements
nvle (:,2)=number of Q1 elements of the extendend elements
along x-direction
nvle (:,3)=number of Q1 elements of the extendend elements
along y-direction
Ac = Cholesky factor of the Q1 stiffness matrix
(internal/internal nodes) on the coarse mesh
Acb = Q1 stiffness matrix (internal/boundary nodes)
on the coarse mesh
wwc = Q1 mass matrix (internal nodes)
on the coarse mesh
r0t = matrix R_H^T
lista_coarse = list of nodes of Omega wich are nodes of the coarse mesh
noec = number of nodes of the coarse mesh
novc(4,ne) = it maps each macro-Q1 element on the coarse mesh Q1
lintc = lista of internal nodes of the coarse mesh
ldirc = lista of boudanry nodes of the coarse mesh
Output: u = column array of PCG solution
iter = number of iterations required to satisfy stopping test
||r^(k)||/||b|| < tol
res = value of ||r^(k)||/||b|| at last iteration.
References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Fundamentals in Single Domains"
Springer Verlag, Berlin Heidelberg New York, 2006.
CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Evolution to Complex Geometries
and Applications to Fluid DynamicsSpectral Methods"
Springer Verlag, Berlin Heidelberg New York, 2007.0001 function [u,iter,res]=schwarz_pbcgstab(u0, f, param,... 0002 p_unity,... 0003 xy, ww, A, nov, noei, lint,x,wx,y,wy,xx,jacx,yy,jacy,... 0004 Aq1, wwq1, linte,nove,nvle,... 0005 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0006 % SCHWARZ_PBCGSTAB BiCGStab method with additive Schwarz preconditioner with overlap and coarse mesh 0007 % 0008 % [u,iter,res]=schwarz_pbcgstab(u0, f, param,... 0009 % p_unity,... 0010 % xy, ww, A, nov, noei, lint,x,wx,y,wy,xx,jacx,yy,jacy,... 0011 % Aq1, wwq1, linte,nove,nvle,... 0012 % Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0013 % 0014 % Input : u0 = column array, initial datum for pbcgstab 0015 % f = column array, r.h.s. of linear system to solve 0016 % param = array of parameters 0017 % p_unity = unity partition associated to extended elements 0018 % xy = 2-indexes array wiht coordinates of 2D LGL mesh 0019 % ww = mass matrix (in column array) 0020 % nov = local -global map, previously generated by cosnov_2d 0021 % noei = number of nodes internal to Omega 0022 % lint = list of internal nodes (set in liste.n) 0023 % x = npdx LGL nodes in [-1,1], 0024 % (produced by calling [x,wx]=xwlgl(npdx)) 0025 % wx= npdx LGL weigths in [-1,1], 0026 % (produced by calling [x,wx]=xwlgl(npdx)) 0027 % y = npdx LGL nodes in [-1,1], 0028 % (produced by calling [y,wy]=xwlgl(npdy)) 0029 % wy= npdx LGL weigths in [-1,1], 0030 % (produced by calling [y,wy]=xwlgl(npdy)) 0031 % xx = 2-indexes array of size (4,ne) 0032 % xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie] 0033 % (ne=nex*ney is the global number of spectral elements) 0034 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0035 % yy = 2-indexes array of size (4,ne): 0036 % yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie] 0037 % jacy = array (length(jacy)=ne); jacy(ie)= (x_V3_ie-x_V1_ie)/2 0038 % Aq1 = cell array with Q1 stiffness matrices (internal/internal nodes) 0039 % on extended elements 0040 % wwq1 = cell array with Q1 mass matrices (internal nodes) 0041 % on extended elements 0042 % linte = cell array with list of internal nodes of extended 0043 % elements 0044 % nove = 2-index array of "extended local" to global map, permuted. 0045 % nvle = 2-index array: 0046 % nvle (:,1)=number of nodes of the extendend elements 0047 % nvle (:,2)=number of Q1 elements of the extendend elements 0048 % along x-direction 0049 % nvle (:,3)=number of Q1 elements of the extendend elements 0050 % along y-direction 0051 % Ac = Cholesky factor of the Q1 stiffness matrix 0052 % (internal/internal nodes) on the coarse mesh 0053 % Acb = Q1 stiffness matrix (internal/boundary nodes) 0054 % on the coarse mesh 0055 % wwc = Q1 mass matrix (internal nodes) 0056 % on the coarse mesh 0057 % r0t = matrix R_H^T 0058 % lista_coarse = list of nodes of Omega wich are nodes of the coarse mesh 0059 % noec = number of nodes of the coarse mesh 0060 % novc(4,ne) = it maps each macro-Q1 element on the coarse mesh Q1 0061 % lintc = lista of internal nodes of the coarse mesh 0062 % ldirc = lista of boudanry nodes of the coarse mesh 0063 % 0064 % Output: u = column array of PCG solution 0065 % iter = number of iterations required to satisfy stopping test 0066 % ||r^(k)||/||b|| < tol 0067 % res = value of ||r^(k)||/||b|| at last iteration. 0068 % 0069 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0070 % "Spectral Methods. Fundamentals in Single Domains" 0071 % Springer Verlag, Berlin Heidelberg New York, 2006. 0072 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0073 % "Spectral Methods. Evolution to Complex Geometries 0074 % and Applications to Fluid DynamicsSpectral Methods" 0075 % Springer Verlag, Berlin Heidelberg New York, 2007. 0076 0077 % Written by Paola Gervasio 0078 % $Date: 2007/04/01$ 0079 0080 tol=param(10); itmax=param(11); 0081 [ldnov,ne]=size(nov); 0082 u=u0; 0083 bb=norm(f); 0084 res=[]; 0085 r=f-A*u0; 0086 rtilde=r; 0087 p=r; 0088 rr0=r'*r; 0089 normar0=sqrt(r'*r); 0090 err=normar0/bb; iter=0; 0091 0092 0093 while err > tol & iter < itmax 0094 0095 z=precoasc(p,param,noei,lint,p_unity,... 0096 xy, ww, nov, x,wx,y,wy,xx, jacx,yy,jacy,... 0097 Aq1, wwq1, linte, nove,nvle,... 0098 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0099 %v=A*p; 0100 v=A*z; 0101 alpha=rr0/(rtilde'*v); 0102 s=r-alpha*v; 0103 ss=precoasc(s,param,noei,lint,p_unity,... 0104 xy, ww, nov, x,wx,y,wy,xx, jacx,yy,jacy,... 0105 Aq1, wwq1, linte, nove,nvle,... 0106 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0107 t=A*ss; 0108 w=(t'*s)/(t'*t); 0109 u=u+alpha*z+w*ss; 0110 r=s-w*t; 0111 rr=rtilde'*r; 0112 beta=rr/rr0; 0113 p=r+(rr*alpha)/(rr0*w)*(p-w*v); 0114 rr0=rr; 0115 err=sqrt(r'*r)/normar0; 0116 if(abs(alpha)<eps) err=0; 0117 disp('schwarz_pcg: a break occurred') 0118 end 0119 0120 %fprintf('iter=%d alpha=%12.6e err=%12.6e \n',iter,alpha,err) 0121 %res=[res;err]; 0122 iter=iter+1; 0123 rr0=rr; 0124 0125 0126 %fprintf('It. Bicgstab %d, err %13.6e \n',iter,err) 0127 end 0128 if err>tol & iter>=itmax 0129 fprintf('The method does not converge in %d iterations\n',iter) 0130 end 0131 res=err; 0132 return