SCHUR_2D Numerical solution of the 2D b.v.p. -Delta u + gam u -Schur complement
-Delta u + gam u = f x in Omega
+ Dirichlet bc
by Schur Complement Matrix (Sect. 6.4.3, CHQZ3)
and by SEM Numerical Integration with LGL quadrature formulas.
[xy,un,param]=schur_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 : No preconditioner (P=I)
2 : P= Neumann-Neumann (Alg. 6.4.3, CHQZ3)
3 : P= balancing Neumann-Neumann (Alg. 6.4.4, CHQZ3)
param(2) : NOT USED
param(3) : NOT USED
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 stopping test
param(11)= maximum number of iterations for PCG stopping test
param(12:20) NOT USED
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]=schur_2d(xa,xb,ya,yb,gam,... 0002 uex,uex_x,uex_y,ff,g,h,cb,nex,nx,ney,ny,gammax,gammay,param); 0003 % SCHUR_2D Numerical solution of the 2D b.v.p. -Delta u + gam u -Schur complement 0004 % 0005 % -Delta u + gam u = f x in Omega 0006 % 0007 % + Dirichlet bc 0008 % 0009 % by Schur Complement Matrix (Sect. 6.4.3, CHQZ3) 0010 % and by SEM Numerical Integration with LGL quadrature formulas. 0011 % 0012 % [xy,un,param]=schur_2d(xa,xb,ya,yb,gam,... 0013 % uex,uex_x,uex_y,ff,nex,nx,ney,ny,gammax,gammay,param); 0014 % 0015 % Omega=(xa,xb) x (ya,yb) 0016 % 0017 % side 3 0018 % V4 __________ V3 0019 % | | 0020 % | | 0021 % side 4 | Omega | side 2 0022 % | | 0023 % |__________| 0024 % V1 V2 0025 % side 1 0026 % 0027 % 0028 % __________________________ 0029 % | | | | | 0030 % | 3 | 6 | 9 | 12 | Omega and spectral elements 0031 % | | | | | ordering 0032 % __________________________ 0033 % | | | | | 0034 % | 2 | 5 | 8 | 11 | 0035 % | | | | | 0036 % __________________________ 0037 % | | | | | 0038 % | 1 | 4 | 7 | 10 | 0039 % | | | | | 0040 % __________________________ 0041 % 0042 % 0043 % Input: xa= abscissa of either vertex V1 or vertex V4 0044 % xb= abscissa of either vertex V2 or vertex V3 0045 % ya= ordinate of either vertex V1 or vertex V2 0046 % yb= ordinate of either vertex V3 or vertex V4 0047 % gam = coefficient of zeroth order term (constant>=0) 0048 % uex = exact solution (uex=@(x,y)[uex(x,y)], with .*, .^, ./) 0049 % uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)], with .*, .^, ./) 0050 % uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)], with .*, .^, ./) 0051 % ff = r.h.s. solution (ff=@(x,y)[ff(x,y)], with .*, .^, ./) 0052 % g = @(x,y)[....] function handle to the expression of Dirichlet 0053 % boundary data 0054 % h =@(x,y)[....] function handle to the expression of Neumann 0055 % boundary data. It is a vector of 4 functions, each for any 0056 % side. 0057 % cb= a string of four elements, each for any side of \partial\Omega. 0058 % cb(i) refers to side number "i" of \partial\Omega 0059 % 'd' character stands for Dirichlet boundary conditions on that side 0060 % 'n' character stands for Neumann boundary conditions on that side 0061 % It works only if cb='dddd'; 0062 % nex = number of elements (equally spaced) along x-direction 0063 % nx = polynomial degree in each element (the same in each element) 0064 % along x-direction 0065 % ney = number of elements (equally spaced) along y-direction 0066 % ny = polynomial degree in each element (the same in each element) 0067 % along y-direction% gammax = column or row array of length nex-1. 0068 % If the deomposition of Omega is not 0069 % uniform, gammax is the vector of position of interfaces between 0070 % spectral elements along x-direction. If gammax=[], uniform 0071 % decomposition is used. 0072 % gammay = column or row array of length ney-1. 0073 % If the deomposition of Omega is not 0074 % uniform, gammay is the vector of position of interfaces between 0075 % spectral elements along y-direction. If gammay=[], uniform 0076 % decomposition is used. 0077 % param = array of parameters 0078 % param(1) = 1 : No preconditioner (P=I) 0079 % 2 : P= Neumann-Neumann (Alg. 6.4.3, CHQZ3) 0080 % 3 : P= balancing Neumann-Neumann (Alg. 6.4.4, CHQZ3) 0081 % param(2) : NOT USED 0082 % param(3) : NOT USED 0083 % param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm) 0084 % on the exact solution 0085 % 2: no errors are computed 0086 % param(5) = 0: LG quadrature formulas with high precision degree are 0087 % used to compute norms (exact norms) 0088 % 1: LGL quadrature formulas with npdx,npdy nodes are 0089 % used to compute norms (discrete norms) 0090 % (used only if param(4) == 1) 0091 % param(6) = number of nodes for high degree quadrature formula, 0092 % (used only if param(5) == 0 & param(4) == 1) 0093 % param(7) = 0: absolute errors are computed 0094 % 1: relative errors are computed 0095 % (used only if param(4) == 1) 0096 % param(8) = 0: not plot the solution 0097 % 1: plot the solution (mesh) 0098 % 2: plot the solution (surf) 0099 % 3: plot the solution (contour) 0100 % param(9) = number of nodes in each element and along each direction 0101 % for plotting the solution 0102 % (used only if param(8) > 0) 0103 % param(10)= tolerance for PCG stopping test 0104 % param(11)= maximum number of iterations for PCG stopping test 0105 % param(12:20) NOT USED 0106 % 0107 % Output: xy = 2-indexes array wiht coordinates of 2D LGL mesh 0108 % un = numerical solution 0109 % param(21) = (if param(3)==1) iterations of PCG 0110 % param(22) = (if param(3)==1) final residual 0111 % param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of 0112 % the preconditioned Schur complement 0113 % param(25,26,27) = (if param(4)==1 & param(3)==1) 0114 % errors on the exact solution: 0115 % L^inf-norm, L2-norm, H1-norm 0116 % 0117 % 0118 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0119 % "Spectral Methods. Fundamentals in Single Domains" 0120 % Springer Verlag, Berlin Heidelberg New York, 2006. 0121 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0122 % "Spectral Methods. Evolution to Complex Geometries 0123 % and Applications to Fluid DynamicsSpectral Methods" 0124 % Springer Verlag, Berlin Heidelberg New York, 2007. 0125 0126 % Written by Paola Gervasio 0127 % $Date: 2007/04/01$ 0128 0129 xy=[]; un=[]; 0130 0131 npdx=nx+1; npdy=ny+1; ldnov=npdx*npdy; mn=ldnov; ne=nex*ney; 0132 0133 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0134 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0135 0136 % Generation of nov: its implements the 0137 % restriction map R_m from \overline\Omega to overline\Omega_m 0138 % (CHQZ, pag. 394) 0139 0140 [nov]=cosnov_2d(npdx,nex,npdy,ney); 0141 noe=nov(ldnov,ne); 0142 0143 % Mesh generation 0144 0145 [xx,yy,jacx,jacy,xy,ww,ifro]=mesh_2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0146 nov,x,wx,y,wy,gammax,gammay); 0147 0148 0149 % Generation of lists of internal, boundary, interface nodes. All lists are referred to 0150 % global ordering on Omega 0151 % 0152 % ldir: list of Dirichlet boundary nodes 0153 % lint: list of internal nodes 0154 % lintint: list of internal nodes which are internal to spectral elements 0155 % lgamma: list of internal nodes which are on the interfaces between spectral elements 0156 0157 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0158 [novi,nvli]=cosnovi(nov,ifro,lint); 0159 0160 % Rigth Hand Side 0161 f=ff(xy(:,1),xy(:,2)).*ww; 0162 ub=ifro.*uex(xy(:,1),xy(:,2)); 0163 0164 % Generation of local matrices and local r.h.s 0165 0166 [AGG,Amm,AGm,Lmm,LGG,Am,f]=schur_loc(ifro,nov,wx,dx,jacx,... 0167 wy,dy,jacy,nvli,gam,f,ub,param); 0168 0169 % Restrict r.h.s to internal nodes 0170 f=f(lint); 0171 0172 noei=length(lint); 0173 ifroi=ifro(lint); 0174 xyi=[xy(lint,1),xy(lint,2)]; 0175 wwi=ww(lint); 0176 0177 % Generation of novg: its implements the 0178 % restriction map R_{\Gamma_m} from \Gamma to \Gamma_m (= 0179 % \partial\Omega_m\cap\Gamma) (CHQZ, pag. 394) 0180 0181 ngamma=length(lgamma); 0182 [novg]=cosnovg(xyi,noei,ifroi,lgamma,ldnov,novi,nvli); 0183 0184 % Rgamma = matrix of size (ne, ngamma) (defined in CHQZ3, pag. 399) 0185 % (R_\Gamma)_ij: = 1/n_j if x_j \in \Gamma_i 0186 % = 0 otherwise 0187 0188 [Rgamma]=cosrgam(novg,LGG,ne,ngamma); 0189 fsigma=f(lgamma); 0190 0191 % compute f_sigma with formula (6.4.31) pag. 394 CHQZ3 0192 for ie=1:ne 0193 ngammam=length(LGG{ie}); 0194 agm=AGm{ie};amm=Amm{ie};lmm=Lmm{ie};lgg=LGG{ie}; 0195 floc=agm*(amm*f(novi(lmm,ie))); 0196 fsigma(novg(lgg,ie))=fsigma(novg(lgg,ie))-floc; 0197 end 0198 0199 0200 u0=zeros(ngamma,1); 0201 if (param(1)==1) 0202 0203 % identity preconditioner 0204 0205 [usigma,iter,res]=schur_pcg(u0, fsigma, param(10),param(11),param,... 0206 AGG,Amm,AGm,LGG,Am,nvli,novg); 0207 0208 elseif param(1)==2 % neumann-neumann preconditioner 0209 0210 % Computes the diagonal weighting matrix D relative to the interface 0211 % unknowns: 0212 % D_ii=1/n_i where n_i is the number of subdomains x_i belongs to. 0213 % D_m=D(novg(LGG_m,m)) is the restriction of D to Gamma_m 0214 % (used in (6.4.40), pag. 397, CHQZ3 0215 0216 D=partition(Rgamma); 0217 [usigma,iter,res]=schur_pcg(u0, fsigma, param(10),param(11),param,... 0218 AGG,Amm,AGm,LGG,Am,nvli,novg,D); 0219 0220 elseif param(1)==3 % balancing neumann-neumann preconditioner 0221 0222 % Computes the diagonal weighting matrix D relative to the interface 0223 % unknowns: 0224 % D_ii=1/n_i where n_i is the number of subdomains x_i belongs to. 0225 % D_m=D(novg(LGG_m,m)) is the restriction of D to Gamma_m 0226 % (used in (6.4.40), pag. 397, CHQZ3 0227 0228 D=partition(Rgamma); 0229 0230 % Computes PSH= pseudoinverse of Sigma_H 0231 0232 [PSH]=pinv_sigma(AGG,Amm,AGm,LGG,novg,Rgamma); 0233 0234 [usigma,iter,res]=schur_pcg(u0, fsigma, param(10),param(11),param,... 0235 AGG,Amm,AGm,LGG,Am,nvli,novg,D,Rgamma,PSH); 0236 end 0237 0238 0239 param(21)=iter; 0240 param(22)=res; 0241 0242 % computes local solutions in Omega_ie 0243 0244 [un]=local_solver(Amm,AGm,Lmm,LGG,novi,nvli,nov,novg,lint,lgamma,usigma,f,ub); 0245 0246 0247 if (param(4)==1) 0248 % errors computing on the exact solution 0249 [err_inf,err_h1,err_l2]=errors_2d(x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0250 xy,ww,nov,un,uex,uex_x,uex_y,param); 0251 param(25:27)=[err_inf;err_h1;err_l2]; 0252 end 0253 0254 if (param(8)>0) 0255 if param(8)==1 0256 command='mesh'; 0257 elseif param(8)==2 0258 command='surf'; 0259 elseif param(8)==3; 0260 command='contour' 0261 end 0262 % plot the solution 0263 fig=figure(..., 0264 'Name','SEM-NI solution of -Delta u + gam u=f, Dirichlet bc',... 0265 'Visible','on'); 0266 [hf,ha]=plot_sem(fig,command,nex,ney,x,xx,jacx,y,yy,jacy,xy,ww,nov,un,param); 0267 0268 end