SCHWARZ_2D Numerical solution of the 2D boundary value problem
-Delta u + gam u = f x in Omega
+ Dirichlet bc
by SEM Numerical Integration with LGL quadrature formulas, with
additive Schwarz preconditioner with overlapping elements and
coarse mesh (Sect. 6.3.3, CHQZ3)
[xy,un,param]=schwarz_2d(xa,xb,ya,yb,gam,...
uex,uex_x,uex_y,ff,nex,nx,ney,ny,gammax,gammay,param);
Omega=(xa,xb) x (ya,yb)
side 3
V4 __________ V3
| |
| |
side 4 | Omega | side 2
| |
|__________|
V1 V2
side 1
__________________________
| | | | |
| 3 | 6 | 9 | 12 | Omega and spectral elements
| | | | | ordering
__________________________
| | | | |
| 2 | 5 | 8 | 11 |
| | | | |
__________________________
| | | | |
| 1 | 4 | 7 | 10 |
| | | | |
__________________________
Input: xa= abscissa of either vertex V1 or vertex V4
xb= abscissa of either vertex V2 or vertex V3
ya= ordinate of either vertex V1 or vertex V2
yb= ordinate of either vertex V3 or vertex V4
gam = coefficient of zeroth order term (constant>=0)
uex = exact solution (uex=@(x,y)[uex(x,y)], with .*, .^, ./)
uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)], with .*, .^, ./)
uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)], with .*, .^, ./)
ff = r.h.s. solution (ff=@(x,y)[ff(x,y)], with .*, .^, ./)
g = @(x,y)[....] function handle to the expression of Dirichlet
boundary data
h =@(x,y)[....] function handle to the expression of Neumann
boundary data. It is a vector of 4 functions, each for any
side.
cb= a string of four elements, each for any side of \partial\Omega.
cb(i) refers to side number "i" of \partial\Omega
'd' character stands for Dirichlet boundary conditions on that side
'n' character stands for Neumann boundary conditions on that side
It works only if cb='dddd'
nex = number of elements (equally spaced) along x-direction
nx = polynomial degree in each element (the same in each element)
along x-direction
ney = number of elements (equally spaced) along y-direction
ny = polynomial degree in each element (the same in each element)
along y-direction
gammax = column or row array of length nex-1.
If the deomposition of Omega is not
uniform, gammax is the vector of position of interfaces between
spectral elements along x-direction. If gammax=[], uniform
decomposition is used.
gammay = column or row array of length ney-1.
If the deomposition of Omega is not
uniform, gammay is the vector of position of interfaces between
spectral elements along y-direction. If gammay=[], uniform
decomposition is used.
param = array of parameters
param(1) =1 P=I
=2 P=Additive Schwarz with overlap and coarse mesh
param(2) = number of added layers for extending spectral elements
inside additive Schwarz preconditioner.
If param(2)==1, the preconditioner is P^{m}_{as,H}
(minimum overlap), pag. 377 CHQZ3)
If param(2)==2, the preconditioner is P^{s}_{as,H}
(small overlap), pag. 377 CHQZ3)
param(2) is a positive integer less than min(nx,ny)
param(3) = 1 PCG is used to solve linear system
param(3) = 2 PBiCGStab is used to solve linear system
param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm)
on the exact solution
2: no errors are computed
param(5) = 0: LG quadrature formulas with high precision degree are
used to compute norms (exact norms)
1: LGL quadrature formulas with npdx,npdy nodes are
used to compute norms (discrete norms)
(used only if param(4) == 1)
param(6) = number of nodes for high degree quadrature formula,
(used only if param(5) == 0 & param(4) == 1)
param(7) = 0: absolute errors are computed
1: relative errors are computed
(used only if param(4) == 1)
param(8) = 0: not plot the solution
1: plot the solution (mesh)
2: plot the solution (surf)
3: plot the solution (contour)
param(9) = number of nodes in each element and along each direction
for plotting the solution
(used only if param(8) > 0)
param(10)= tolerance for PCG/PBCGSTAB stopping test
param(11)= maximum number of iterations for PCG/PBCGSTAB stopping test
param(12:20) INTERNAL USE
Output: xy = 2-indexes array wiht coordinates of 2D LGL mesh
un = numerical solution
param(21) = (if param(3)==1) iterations of PCG
param(22) = (if param(3)==1) final residual
param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of
the preconditioned Schur complement
param(25,26,27) = (if param(4)==1 & param(3)==1)
errors on the exact solution:
L^inf-norm, L2-norm, H1-norm
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 [xy,un,param]=schwarz_2d(xa,xb,ya,yb,gam,... 0002 uex,uex_x,uex_y,ff,g,h,cb,nex,nx,ney,ny,gammax,gammay,param); 0003 % SCHWARZ_2D Numerical solution of the 2D boundary value problem 0004 % 0005 % -Delta u + gam u = f x in Omega 0006 % 0007 % + Dirichlet bc 0008 % 0009 % by SEM Numerical Integration with LGL quadrature formulas, with 0010 % additive Schwarz preconditioner with overlapping elements and 0011 % coarse mesh (Sect. 6.3.3, CHQZ3) 0012 % 0013 % [xy,un,param]=schwarz_2d(xa,xb,ya,yb,gam,... 0014 % uex,uex_x,uex_y,ff,nex,nx,ney,ny,gammax,gammay,param); 0015 % 0016 % Omega=(xa,xb) x (ya,yb) 0017 % 0018 % side 3 0019 % V4 __________ V3 0020 % | | 0021 % | | 0022 % side 4 | Omega | side 2 0023 % | | 0024 % |__________| 0025 % V1 V2 0026 % side 1 0027 % 0028 % 0029 % __________________________ 0030 % | | | | | 0031 % | 3 | 6 | 9 | 12 | Omega and spectral elements 0032 % | | | | | ordering 0033 % __________________________ 0034 % | | | | | 0035 % | 2 | 5 | 8 | 11 | 0036 % | | | | | 0037 % __________________________ 0038 % | | | | | 0039 % | 1 | 4 | 7 | 10 | 0040 % | | | | | 0041 % __________________________ 0042 % 0043 % 0044 % Input: xa= abscissa of either vertex V1 or vertex V4 0045 % xb= abscissa of either vertex V2 or vertex V3 0046 % ya= ordinate of either vertex V1 or vertex V2 0047 % yb= ordinate of either vertex V3 or vertex V4 0048 % gam = coefficient of zeroth order term (constant>=0) 0049 % uex = exact solution (uex=@(x,y)[uex(x,y)], with .*, .^, ./) 0050 % uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)], with .*, .^, ./) 0051 % uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)], with .*, .^, ./) 0052 % ff = r.h.s. solution (ff=@(x,y)[ff(x,y)], with .*, .^, ./) 0053 % g = @(x,y)[....] function handle to the expression of Dirichlet 0054 % boundary data 0055 % h =@(x,y)[....] function handle to the expression of Neumann 0056 % boundary data. It is a vector of 4 functions, each for any 0057 % side. 0058 % cb= a string of four elements, each for any side of \partial\Omega. 0059 % cb(i) refers to side number "i" of \partial\Omega 0060 % 'd' character stands for Dirichlet boundary conditions on that side 0061 % 'n' character stands for Neumann boundary conditions on that side 0062 % It works only if cb='dddd' 0063 % nex = number of elements (equally spaced) along x-direction 0064 % nx = polynomial degree in each element (the same in each element) 0065 % along x-direction 0066 % ney = number of elements (equally spaced) along y-direction 0067 % ny = polynomial degree in each element (the same in each element) 0068 % along y-direction 0069 % gammax = column or row array of length nex-1. 0070 % If the deomposition of Omega is not 0071 % uniform, gammax is the vector of position of interfaces between 0072 % spectral elements along x-direction. If gammax=[], uniform 0073 % decomposition is used. 0074 % gammay = column or row array of length ney-1. 0075 % If the deomposition of Omega is not 0076 % uniform, gammay is the vector of position of interfaces between 0077 % spectral elements along y-direction. If gammay=[], uniform 0078 % decomposition is used. 0079 % param = array of parameters 0080 % param(1) =1 P=I 0081 % =2 P=Additive Schwarz with overlap and coarse mesh 0082 % param(2) = number of added layers for extending spectral elements 0083 % inside additive Schwarz preconditioner. 0084 % If param(2)==1, the preconditioner is P^{m}_{as,H} 0085 % (minimum overlap), pag. 377 CHQZ3) 0086 % If param(2)==2, the preconditioner is P^{s}_{as,H} 0087 % (small overlap), pag. 377 CHQZ3) 0088 % param(2) is a positive integer less than min(nx,ny) 0089 % param(3) = 1 PCG is used to solve linear system 0090 % param(3) = 2 PBiCGStab is used to solve linear system 0091 % param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm) 0092 % on the exact solution 0093 % 2: no errors are computed 0094 % param(5) = 0: LG quadrature formulas with high precision degree are 0095 % used to compute norms (exact norms) 0096 % 1: LGL quadrature formulas with npdx,npdy nodes are 0097 % used to compute norms (discrete norms) 0098 % (used only if param(4) == 1) 0099 % param(6) = number of nodes for high degree quadrature formula, 0100 % (used only if param(5) == 0 & param(4) == 1) 0101 % param(7) = 0: absolute errors are computed 0102 % 1: relative errors are computed 0103 % (used only if param(4) == 1) 0104 % param(8) = 0: not plot the solution 0105 % 1: plot the solution (mesh) 0106 % 2: plot the solution (surf) 0107 % 3: plot the solution (contour) 0108 % param(9) = number of nodes in each element and along each direction 0109 % for plotting the solution 0110 % (used only if param(8) > 0) 0111 % param(10)= tolerance for PCG/PBCGSTAB stopping test 0112 % param(11)= maximum number of iterations for PCG/PBCGSTAB stopping test 0113 % param(12:20) INTERNAL USE 0114 % 0115 % Output: xy = 2-indexes array wiht coordinates of 2D LGL mesh 0116 % un = numerical solution 0117 % param(21) = (if param(3)==1) iterations of PCG 0118 % param(22) = (if param(3)==1) final residual 0119 % param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of 0120 % the preconditioned Schur complement 0121 % param(25,26,27) = (if param(4)==1 & param(3)==1) 0122 % errors on the exact solution: 0123 % L^inf-norm, L2-norm, H1-norm 0124 % 0125 % 0126 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0127 % "Spectral Methods. Fundamentals in Single Domains" 0128 % Springer Verlag, Berlin Heidelberg New York, 2006. 0129 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0130 % "Spectral Methods. Evolution to Complex Geometries 0131 % and Applications to Fluid DynamicsSpectral Methods" 0132 % Springer Verlag, Berlin Heidelberg New York, 2007. 0133 0134 % Written by Paola Gervasio 0135 % $Date: 2007/04/01$ 0136 0137 xy=[]; un=[]; 0138 0139 npdx=nx+1; npdy=ny+1; ldnov=npdx*npdy; mn=ldnov; ne=nex*ney; 0140 0141 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0142 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0143 0144 % Generation of nov: its implements the 0145 % restriction map R_m from \overline\Omega to overline\Omega_m 0146 % (CHQZ, pag. 394) 0147 0148 [nov]=cosnov_2d(npdx,nex,npdy,ney); 0149 noe=nov(ldnov,ne); 0150 0151 % Mesh generation 0152 0153 [xx,yy,jacx,jacy,xy,ww,ifro]=mesh_2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0154 nov,x,wx,y,wy,gammax,gammay); 0155 0156 0157 % Generation of lists of internal, boundary, interface nodes. All lists are referred to 0158 % global ordering on Omega 0159 % 0160 % ldir: list of Dirichlet boundary nodes 0161 % lint: list of internal nodes 0162 % lintint: list of internal nodes which are internal to spectral elements 0163 % lgamma: list of internal nodes which are on the interfaces between spectral elements 0164 0165 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0166 0167 % Matrix assembling 0168 0169 A=stiff_2d_se(npdx,nex,npdy,ney,nov,wx,dx,jacx,wy,dy,jacy); 0170 M=spdiags(ww,0,noe,noe); 0171 0172 if gam ~=0 0173 A=A+gam*M; 0174 end 0175 0176 % Rigth Hand Side 0177 f=ff(xy(:,1),xy(:,2)).*ww; 0178 0179 % Dirichlet boundary conditions 0180 ub=uex(xy(:,1),xy(:,2)); 0181 0182 % 0183 0184 Ab=A(lint,ldir); 0185 A=A(lint,lint); 0186 f=f(lint)-Ab*ub(ldir); 0187 noei=length(lint); 0188 wwi=ww(lint); 0189 0190 % extended elements construction 0191 0192 [nove,nvle]=cosnovenew(nx,nex,ny,ney,nov,ifro,param(2)); 0193 0194 % Construction of local stiffness Q1 matrices on extended elements 0195 0196 [Aq1,Abq1,wwq1,linte,ldire,nove]=stiffq1(ifro,nov,xy,nove,nvle); 0197 0198 % Construction of Q1 stiffness matrix on the coarse grid. 0199 0200 [Ac,Acb,wwc,lista_coarse,noec,novc,lintc,ldirc]=stiffq1H(nx,nex,ny,ney,xy,nov,ifro); 0201 r0t=matr0t(nx,ny,xy,nov,novc,noec,lista_coarse); 0202 0203 % Unity partition 0204 0205 nvl=[mn*ones(ne,1),npdx*ones(ne,1),npdy*ones(ne,1)]; 0206 [p_unity]=partition_e(nov,nvl,noe); 0207 0208 % PCG or PBCGstab with Schwarz preconditioner 0209 0210 u0=zeros(noei,1); 0211 param(13:16)=[nx,ny,nex,ney]; 0212 if param(3)==1 0213 [un,iter,res]=schwarz_pcg(u0, f, param, p_unity,... 0214 xy, ww, A, nov, noei, lint, x,wx, y,wy, xx,jacx,yy,jacy,... 0215 Aq1, wwq1,linte,nove,nvle,... 0216 Ac, Acb,wwc,r0t,lista_coarse,noec,novc,lintc,ldirc); 0217 elseif param(3)==2 0218 [un,iter,res]=schwarz_pbcgstab(u0, f, param, p_unity,... 0219 xy, ww, A, nov, noei, lint, x,wx, y,wy, xx,jacx,yy,jacy,... 0220 Aq1, wwq1,linte,nove,nvle,... 0221 Ac, Acb,wwc,r0t,lista_coarse,noec,novc,lintc,ldirc); 0222 end 0223 param(21:22)=[iter,res]; 0224 ub(lint)=un; 0225 un=ub; 0226 0227 if (param(4)==1) 0228 % errors computing on the exact solution 0229 [err_inf,err_h1,err_l2]=errors_2d(x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0230 xy,ww,nov,un,uex,uex_x,uex_y,param); 0231 param(25:27)=[err_inf;err_h1;err_l2]; 0232 end 0233 0234 if (param(8)>0) 0235 if param(8)==1 0236 command='mesh'; 0237 elseif param(8)==2 0238 command='surf'; 0239 elseif param(8)==3; 0240 command='contour' 0241 end 0242 % plot the solution 0243 fig=figure(..., 0244 'Name','SEM-NI solution of -Delta u + gam u=f, Dirichlet bc',... 0245 'Visible','on'); 0246 [ha]=plot_sem_2d(fig,command,nex,ney,x,xx,jacx,y,yy,jacy,xy,ww,nov,un,param(9)); 0247 0248 end