SCHWARZ_PCG Conjugate Gradient method with additive Schwarz prec with overlap and coarse mesh
[u,iter,res]=schwarz_pcg(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_pcg(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_PCG Conjugate Gradient method with additive Schwarz prec with overlap and coarse mesh 0007 % 0008 % [u,iter,res]=schwarz_pcg(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 tol=param(10); itmax=param(11); 0080 [ldnov,ne]=size(nov); 0081 u=u0; 0082 bb=norm(f); 0083 res=[]; 0084 r=f-A*u0; 0085 if(param(1)==2) 0086 % chiamo il preco Schwarz additivo coarse su mesh Q1 0087 z=precoasc(r,param,noei,lint,p_unity,... 0088 xy, ww, nov, x,wx,y,wy,xx,jacx,yy,jacy,... 0089 Aq1, wwq1, linte, nove,nvle,... 0090 Ac, Acb,wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0091 elseif param(1)==1 0092 z=r; 0093 end 0094 rr0=z'*r; 0095 normar0=sqrt(r'*r); 0096 p=z; 0097 err=normar0/bb; iter=0; 0098 0099 0100 while err > tol & iter < itmax 0101 0102 %v=A*p; 0103 v=A*p; 0104 alpha=(r'*z)/(p'*v); 0105 u=u+alpha*p; 0106 r=r-alpha*v; 0107 if param(1)==2 0108 z=precoasc(r,param,noei,lint,p_unity,... 0109 xy, ww, nov, x,wx,y,wy,xx, jacx,yy,jacy,... 0110 Aq1, wwq1, linte, nove,nvle,... 0111 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0112 elseif param(1)==1 0113 z=r; 0114 end 0115 rr=z'*r; 0116 beta=rr/rr0; 0117 p=z+beta*p; 0118 err=sqrt(r'*r)/normar0; 0119 if(abs(alpha)<eps) err=0; 0120 disp('schwarz_pcg: a break occurred') 0121 end 0122 0123 %fprintf('iter=%d alpha=%12.6e err=%12.6e \n',iter,alpha,err) 0124 %res=[res;err]; 0125 iter=iter+1; 0126 rr0=rr; 0127 0128 0129 %fprintf('It. GC %d, err %13.6e \n',iter,err) 0130 end 0131 if err>tol & iter>=itmax 0132 fprintf('The method does not converge in %d iterations\n',iter) 0133 end 0134 res=err; 0135 return