EIG_SCHWARZ_2D Eigenvalues computation for the matrix associated to 2D b.v.p. -Delta u + gam u = f + Dirichlet b.c.
-Delta u + gam u = f x in Omega
+ Dirichlet bc
SEM Numerical Integration with LGL quadrature formulas, with
additive Schwarz preconditioner with overlapped elements and
coarse mesh (Sect. 6.3.3, CHQZ3)
[param]=eig_schwarz_2d(xa,xb,ya,yb,gam,...
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)
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) : maximum number of iterations in eigs
param(4) : tolerance for stopping test in eigs
param(5:20) NOT USED or INTERNAL USE
Output:
param(11,12) extrema eigenvalues (max,min) of
the preconditioned stiffness matrix
param(13) iterative condition number (|lam_max|/|lam_min|) of
the preconditioned stiffness matrix
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 [param]=eig_schwarz_2d(xa,xb,ya,yb,gam,... 0002 cb,nex,nx,ney,ny,gammax,gammay,param); 0003 % EIG_SCHWARZ_2D Eigenvalues computation for the matrix associated to 2D b.v.p. -Delta u + gam u = f + Dirichlet b.c. 0004 % 0005 % 0006 % -Delta u + gam u = f x in Omega 0007 % 0008 % + Dirichlet bc 0009 % 0010 % SEM Numerical Integration with LGL quadrature formulas, with 0011 % additive Schwarz preconditioner with overlapped elements and 0012 % coarse mesh (Sect. 6.3.3, CHQZ3) 0013 % 0014 % [param]=eig_schwarz_2d(xa,xb,ya,yb,gam,... 0015 % nex,nx,ney,ny,gammax,gammay,param); 0016 % 0017 % Omega=(xa,xb) x (ya,yb) 0018 % 0019 % side 3 0020 % V4 __________ V3 0021 % | | 0022 % | | 0023 % side 4 | Omega | side 2 0024 % | | 0025 % |__________| 0026 % V1 V2 0027 % side 1 0028 % 0029 % 0030 % __________________________ 0031 % | | | | | 0032 % | 3 | 6 | 9 | 12 | Omega and spectral elements 0033 % | | | | | ordering 0034 % __________________________ 0035 % | | | | | 0036 % | 2 | 5 | 8 | 11 | 0037 % | | | | | 0038 % __________________________ 0039 % | | | | | 0040 % | 1 | 4 | 7 | 10 | 0041 % | | | | | 0042 % __________________________ 0043 % 0044 % 0045 % Input: xa= abscissa of either vertex V1 or vertex V4 0046 % xb= abscissa of either vertex V2 or vertex V3 0047 % ya= ordinate of either vertex V1 or vertex V2 0048 % yb= ordinate of either vertex V3 or vertex V4 0049 % gam = coefficient of zeroth order term (constant>=0) 0050 % cb= a string of four elements, each for any side of \partial\Omega. 0051 % cb(i) refers to side number "i" of \partial\Omega 0052 % 'd' character stands for Dirichlet boundary conditions on that side 0053 % 'n' character stands for Neumann boundary conditions on that side 0054 % It works only if cb='dddd' 0055 % nex = number of elements (equally spaced) along x-direction 0056 % nx = polynomial degree in each element (the same in each element) 0057 % along x-direction 0058 % ney = number of elements (equally spaced) along y-direction 0059 % ny = polynomial degree in each element (the same in each element) 0060 % along y-direction% gammax = column or row array of length nex-1. 0061 % If the deomposition of Omega is not 0062 % uniform, gammax is the vector of position of interfaces between 0063 % spectral elements along x-direction. If gammax=[], uniform 0064 % decomposition is used. 0065 % gammay = column or row array of length ney-1. 0066 % If the deomposition of Omega is not 0067 % uniform, gammay is the vector of position of interfaces between 0068 % spectral elements along y-direction. If gammay=[], uniform 0069 % decomposition is used. 0070 % param = array of parameters 0071 % param(1) =1 P=I 0072 % =2 P=Additive Schwarz with overlap and coarse mesh 0073 % param(2) = number of added layers for extending spectral elements 0074 % inside additive Schwarz preconditioner. 0075 % If param(2)==1, the preconditioner is P^{m}_{as,H} 0076 % (minimum overlap), pag. 377 CHQZ3) 0077 % If param(2)==2, the preconditioner is P^{s}_{as,H} 0078 % (small overlap), pag. 377 CHQZ3) 0079 % param(2) is a positive integer less than min(nx,ny) 0080 % param(3) : maximum number of iterations in eigs 0081 % param(4) : tolerance for stopping test in eigs 0082 % param(5:20) NOT USED or INTERNAL USE 0083 % 0084 % Output: 0085 % param(11,12) extrema eigenvalues (max,min) of 0086 % the preconditioned stiffness matrix 0087 % param(13) iterative condition number (|lam_max|/|lam_min|) of 0088 % the preconditioned stiffness matrix 0089 % 0090 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0091 % "Spectral Methods. Fundamentals in Single Domains" 0092 % Springer Verlag, Berlin Heidelberg New York, 2006. 0093 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0094 % "Spectral Methods. Evolution to Complex Geometries 0095 % and Applications to Fluid DynamicsSpectral Methods" 0096 % Springer Verlag, Berlin Heidelberg New York, 2007. 0097 0098 % Written by Paola Gervasio 0099 % $Date: 2007/04/01$ 0100 0101 xy=[]; un=[]; 0102 0103 npdx=nx+1; npdy=ny+1; ldnov=npdx*npdy; mn=ldnov; ne=nex*ney; 0104 0105 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0106 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0107 0108 % Generation of nov: its implements the 0109 % restriction map R_m from \overline\Omega to overline\Omega_m 0110 % (CHQZ, pag. 394) 0111 0112 [nov]=cosnov_2d(npdx,nex,npdy,ney); 0113 noe=nov(ldnov,ne); 0114 0115 % Mesh generation 0116 0117 [xx,yy,jacx,jacy,xy,ww,ifro]=mesh_2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0118 nov,x,wx,y,wy,gammax,gammay); 0119 0120 0121 % Generation of lists of internal, boundary, interface nodes. All lists are referred to 0122 % global ordering on Omega 0123 % 0124 % ldir: list of Dirichlet boundary nodes 0125 % lint: list of internal nodes 0126 % lintint: list of internal nodes which are internal to spectral elements 0127 % lgamma: list of internal nodes which are on the interfaces between spectral elements 0128 0129 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0130 0131 % Matrix assembling 0132 0133 A=stiff_2d_se(npdx,nex,npdy,ney,nov,wx,dx,jacx,wy,dy,jacy); 0134 M=spdiags(ww,0,noe,noe); 0135 0136 if gam ~=0 0137 A=A+gam*M; 0138 end 0139 0140 0141 % 0142 0143 A=A(lint,lint); 0144 noei=length(lint); 0145 wwi=ww(lint); 0146 0147 % extended elements construction 0148 0149 [nove,nvle]=cosnovenew(nx,nex,ny,ney,nov,ifro,param(2)); 0150 0151 % Construction of local stiffness Q1 matrices on extended elements 0152 0153 [Aq1,Abq1,wwq1,linte,ldire,nove]=stiffq1(ifro,nov,xy,nove,nvle); 0154 0155 % Construction of stiffness Q1 matrix on the coarse grid. 0156 0157 [Ac,Acb,wwc,lista_coarse,noec,novc,lintc,ldirc]=stiffq1H(nx,nex,ny,ney,xy,nov,ifro); 0158 r0t=matr0t(nx,ny,xy,nov,novc,noec,lista_coarse); 0159 0160 % Unity partition 0161 0162 nvl=[mn*ones(ne,1),npdx*ones(ne,1),npdy*ones(ne,1)]; 0163 [p_unity]=partition_e(nov,nvl,noe); 0164 param(13:16)=[nx,ny,nex,ney]; 0165 0166 % Matrix construction 0167 0168 Pas=zeros(noei,noei); 0169 for i=1:noei 0170 r=A(:,i); 0171 [z]=precoasc(r,param,noei,lint,p_unity,... 0172 xy, ww, nov, x,wx,y,wy,xx,jacx,yy,jacy,... 0173 Aq1, wwq1, linte, nove, nvle,... 0174 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0175 Pas(:,i)=z; 0176 end 0177 0178 % Eigenvalues computation 0179 0180 if param(4)==0 0181 param(4)=1000; 0182 end 0183 if param(5)==0 0184 param(5)=1.d-12; 0185 end 0186 opts.issym=1;opts.disp=0; 0187 opts.maxit=param(4); opts.tol=param(5); 0188 lambda_max=eigs(Pas,1,'lm',opts); 0189 lambda_min=eigs(Pas,1,'sm',opts); 0190 param(11)=lambda_max; 0191 param(12)=lambda_min; 0192 param(13)=lambda_max/lambda_min; 0193 0194 0195 return