52 integer,
parameter :: VERB = 2
53 logical,
parameter :: WITH_MUMPS = .false.
66 integer ,
parameter :: nev = 6
67 integer ,
parameter :: itMax = 100
68 real(RP),
parameter :: tol = 1
e-7_rp
69 real(RP),
parameter :: sigma =-1.0_rp
72 integer ,
parameter :: nev_exact = 6
73 real(RP),
dimension(nev_exact) :: lambda
93 character(len=100),
parameter :: &
94 & mesh_file= trim(GMSH_DIR)//
"square/square_3.msh" 114 procedure(rntornxl),
pointer :: Kinv => null()
120 character(LEN=300) :: shell_com
133 write(*,*)
'eigen_Laplacian_2D: start' 135 write(*,*)
" RESOLUTION OF THE EIGEN-PROBLEM" 137 write(*,*)
" -Delta u = lambda u on \Omega" 139 write(*,*)
" WITH A HOMOGENEOUS NEUMANN BOUNDARY CONDITION " 141 write(*,*)
" grad u . n = 0 on \partial\Omega" 145 write(*,*)
"================================== ASSEMBLING THE MESH " 146 msh = mesh(mesh_file,
'gmsh')
147 if (verb>1)
call print(msh)
150 write(*,*)
"================================== ASSEMBLING THE& 151 & FINITE ELEMENT SPACE" 155 if (verb>1)
call print(x_h)
158 write(*,*)
"================================== ASSEMBLING THE& 159 & INTEGRATION METHOD ON \Omega" 163 if (verb>1)
call print(qdm)
166 write(*,*)
"================================== ASSEMBLING THE& 167 & MASS AND STIFFNESS MATRICES" 170 call diffusion_massmat(m , one_r3, x_h, qdm)
173 call diffusion_stiffmat(s, emetric, x_h, qdm)
176 write(*,*)
"================================== COMPUTE THE SYSTEM & 177 &MATRIX K = S - sigma * M" 179 call add(k, s, 1._rp, m, -sigma)
182 write(*,*)
"================================== LINEAR SOLVER K*x = y" 185 write(*,*)
" K = LDLT, decomposition with MUMPS" 191 write(*,*)
" Krylov CG settings" 192 kry = krylov(
kry_cg, tol=tol*1.0
e-2_rp, itmax=1000, verb=0)
202 write(*,*)
"================================== EIGEN-SOLVER DEF." 206 eig = eigen(
size =
size, nev = nev, sym=.true., std=.false.)
209 call set(eig, itmax=itmax, tol=tol, verb=2)
212 if (verb>1)
call print(eig)
217 write(*,*)
"================================== NUMERICAL RESOLUTION" 220 call solve(eig,
b=
pmv_m, op=kinv)
225 write(*,*)
"================================== POST-PROCESSING" 227 write(*,*)
"Computed eigenvalues" 229 write(*,*)
" i =", int(ii,1),
", lambda = ", &
230 &
real(eig%lambda(ii)/PI**2, SP),
"* Pi**2" 234 write(*,*)
"Numerical errors |lambda - lambda_h| " 241 lambda(4) = 2.0_rp * pi**2
242 lambda(5) = 4.0_rp * pi**2
243 lambda(6) = 5.0_rp * pi**2
245 do ii=1, min(nev, nev_exact)
246 err = abs( eig%lambda(ii) - lambda(ii))
247 write(*,*)
" i =", int(ii,1),
", error = ", &
256 shell_com =
'eig_Laplace.msh' 259 call write(x_h, trim(shell_com),
'gmsh')
264 call gmsh_addview(x_h, eig%v(:,ii), &
266 &
'Eigen-functions',
real(ii, RP), ii)
270 shell_com =
'gmsh -option '//trim(gmsh_dir)//
'gmsh-options-view'&
271 & //
' '//trim(shell_com)
272 call system(shell_com)
282 subroutine pmv_m(y, x)
283 real(RP),
dimension(:),
intent(out) :: y
284 real(RP),
dimension(:),
intent(in) :: x
286 call matvecprod(y, m, x)
293 real(RP),
dimension(:),
intent(inout) :: y
294 logical ,
intent(out) :: ierr
295 real(RP),
dimension(:),
intent(in) :: x
297 call solve(y, kry, x, k)
307 real(RP),
dimension(:),
intent(inout) :: x
308 logical ,
intent(out) :: ierr
309 real(RP),
dimension(:),
intent(in) :: y
319 real(RP),
dimension(size) :: u, v, w
324 write(*,*)
"Residual | Ku - lambda Mu |_l2 / | u |_l2" 330 call matvecprod(v, k , u)
333 w = v - eig%lambda(ii)*w
334 res = sqrt( sum(w*w) / sum(u*u) )
336 print*,
" i = ", int(ii,1),
" : res = ",
real(res, sp)
341 write(*,*)
"Residual | \nu u - K^{-1} M u |_l2 / | u |_l2" 347 call kinv(w, ierr, v)
350 nu = 1.0_rp / ( eig%lambda(ii) - sigma )
352 res = sqrt( sum(w*w) / sum(u*u) )
354 print*,
" i = ", int(ii,1),
" : res = ",
real(res, sp)
integer, parameter eig_shift_invert
Shift-invert mode.
subroutine kinv_mumps(x, ierr, y)
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
integer, parameter fe_p2_2d
subroutine pmv_m(y, x)
Matrix-vector product x –> M*x.
real(rp) function e(x, v1, v2)
program eigen_laplacian_2d
Solves the Laplacian eigenvalue problem
integer, parameter pc_jacobi
Jacobi diagonal preconditioning.
subroutine kinv_pcg(y, ierr, x)
integer, parameter mumps_ldlt_sdp
MUMPS LDLT-SDP, sym def pos mat.
integer, parameter eig_which_sm
Computation of eigenvalues xith smallest magnitude.
integer, parameter quad_gauss_trg_12
integer, parameter kry_cg
CG linear solver.
real(rp) function b(x, w1, w2)