Choral
eigen_Laplacian_2D.f90
Go to the documentation of this file.
1 
32 
34 
35  use choral
37 
38  implicit none
39 
40  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42  !!
43  !! VARAIBLE DEFINITION : BEGIN
44  !!
45  !!
46 
47  !! GENERAL SETTINGS
48  !!
49  !! verb = verbosity level
50  !! WITH_MUMPS = use MUMPS direct solver for linear system inversions
51  !!
52  integer, parameter :: VERB = 2
53  logical, parameter :: WITH_MUMPS = .false.
54 
55  !! Arnoldi algorithm parameters
56  !! nev = number of computed eigenvalues
57  !! size = problem size (number of DOF)
58  !! sigma = shift factor
59  !! tol = tolerance for resolution
60  !! itMax = maximal number of iterations
61  !!
62  !! lambda = exact eigenvalues (analytic)
63  !!
64  type(eigen) :: eig
65  !!
66  integer , parameter :: nev = 6
67  integer , parameter :: itMax = 100
68  real(RP), parameter :: tol = 1e-7_rp
69  real(RP), parameter :: sigma =-1.0_rp
70  integer :: size
71  !!
72  integer , parameter :: nev_exact = 6
73  real(RP), dimension(nev_exact) :: lambda
74 
75 
76  !! SPACE DISCRETISATION
77  !!
78  !! fe = finite element method
79  !! quad = quadrature method
80  !!
81  integer, parameter :: fe = fe_p2_2d
82  integer, parameter :: quad = quad_gauss_trg_12
83  !!
84  !! msh = mesh
85  !! X_h = finite element space
86  !! qdm = integration method
87  !!
88  type(mesh) :: msh
89  type(fespace) :: X_h
90  type(quadmesh) :: qdm
91  !!
92  !! msh_file = mesh file
93  character(len=100), parameter :: &
94  & mesh_file= trim(GMSH_DIR)//"square/square_3.msh"
95 
96 
97  !! LINEAR SYSTEM
98  !!
99  !! M = mass matrix
100  !! stiff = stiffness matrix
101  !! K = matrix S - sigma*M, sigma the shift factor
102  !!
103  !! lu = mumps decomposition for K ( if WITH_MUMPS )
104  !! kry = krylov method def. ( else )
105  !! pc = preconditioning for the CG krylov
106  !!
107  !! Kinv = pointer to the system inversion K x = y
108  !!
109  type(csr) :: M, S, K
110  type(mumps) :: LU
111  type(krylov) :: kry
112  type(precond) :: pc
113  !!
114  procedure(rntornxl), pointer :: Kinv => null()
115 
116  !! OTHER VARS
117  !!
118  real(RP) :: err
119  integer :: ii
120  character(LEN=300) :: shell_com
121 
122  !!
123  !! VARAIBLE DEFINITION : END
124  !!
125  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
126  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
127 
128 
129  !! !!!!!!!!!!!!!!!!!!!!! INITIALISATION
130  !!
131  call choral_init(verb=verb)
132  write(*,*)
133  write(*,*)'eigen_Laplacian_2D: start'
134  write(*,*)""
135  write(*,*)" RESOLUTION OF THE EIGEN-PROBLEM"
136  write(*,*)""
137  write(*,*)" -Delta u = lambda u on \Omega"
138  write(*,*)""
139  write(*,*)" WITH A HOMOGENEOUS NEUMANN BOUNDARY CONDITION "
140  write(*,*)""
141  write(*,*)" grad u . n = 0 on \partial\Omega"
142  write(*,*)""
143 
144  write(*,*) ""
145  write(*,*)"================================== ASSEMBLING THE MESH "
146  msh = mesh(mesh_file, 'gmsh')
147  if (verb>1) call print(msh)
148 
149  write(*,*) ""
150  write(*,*)"================================== ASSEMBLING THE&
151  & FINITE ELEMENT SPACE"
152  x_h = fespace(msh)
153  call set(x_h, fe)
154  call assemble(x_h)
155  if (verb>1) call print(x_h)
156 
157  write(*,*) ""
158  write(*,*)"================================== ASSEMBLING THE&
159  & INTEGRATION METHOD ON \Omega"
160  qdm = quadmesh(msh)
161  call set(qdm, quad)
162  call assemble(qdm)
163  if (verb>1) call print(qdm)
164 
165  write(*,*) ""
166  write(*,*)"================================== ASSEMBLING THE&
167  & MASS AND STIFFNESS MATRICES"
168  !!
169  !! 'one_R3' is the density function equal to 1 on \Omega
170  call diffusion_massmat(m , one_r3, x_h, qdm)
171  !!
172  !! 'EMetric' is the Euclidian metric on \Omega
173  call diffusion_stiffmat(s, emetric, x_h, qdm)
174 
175  write(*,*) ""
176  write(*,*)"================================== COMPUTE THE SYSTEM &
177  &MATRIX K = S - sigma * M"
178  !!
179  call add(k, s, 1._rp, m, -sigma)
180 
181  write(*,*) ""
182  write(*,*)"================================== LINEAR SOLVER K*x = y"
183  !!
184  if (with_mumps) then
185  write(*,*) " K = LDLT, decomposition with MUMPS"
186  lu = mumps(k, mumps_ldlt_sdp, verb = 0)
187 
188  kinv => kinv_mumps
189 
190  else
191  write(*,*) " Krylov CG settings"
192  kry = krylov(kry_cg, tol=tol*1.0e-2_rp, itmax=1000, verb=0)
193 
194  pc = precond(k, pc_jacobi)
195 
196  kinv => kinv_pcg
197 
198  end if
199 
200 
201  write(*,*) ""
202  write(*,*)"================================== EIGEN-SOLVER DEF."
203  !!
204  !! create the eigen-solver 'eig'
205  size = x_h%nbDof
206  eig = eigen(size = size, nev = nev, sym=.true., std=.false.)
207  !!
208  !! set the eigen-solver 'eig'
209  call set(eig, itmax=itmax, tol=tol, verb=2)
210  call set(eig, which = eig_which_sm)
211  call set(eig, mode = eig_shift_invert, shift = sigma)
212  if (verb>1) call print(eig)
213  write(*,*) ""
214 
215 
216  write(*,*) ""
217  write(*,*)"================================== NUMERICAL RESOLUTION"
218  !!
219  !! Numerical resolution
220  call solve(eig, b=pmv_m, op=kinv)
221  !!
222  call residual()
223 
224  write(*,*) ""
225  write(*,*)"================================== POST-PROCESSING"
226  print*
227  write(*,*) "Computed eigenvalues"
228  do ii=1, nev
229  write(*,*) " i =", int(ii,1), ", lambda = ", &
230  & real(eig%lambda(ii)/PI**2, SP), "* Pi**2"
231  end do
232 
233  write(*,*) ""
234  write(*,*) "Numerical errors |lambda - lambda_h| "
235  !!
236  !! Exact eigen values
237  !!
238  lambda(1) = 0._rp
239  lambda(2) = pi**2
240  lambda(3) = pi**2
241  lambda(4) = 2.0_rp * pi**2
242  lambda(5) = 4.0_rp * pi**2
243  lambda(6) = 5.0_rp * pi**2
244  !!
245  do ii=1, min(nev, nev_exact)
246  err = abs( eig%lambda(ii) - lambda(ii))
247  write(*,*) " i =", int(ii,1), ", error = ", &
248  & real(err, sp)
249  end do
250 
251  if (verb>1) then
252  !!
253  !! graphical display
254  !!
255  !! output file def.
256  shell_com = 'eig_Laplace.msh'
257  !!
258  !! first write the 'finite element mesh'
259  call write(x_h, trim(shell_com), 'gmsh')
260  !!
261  !! second save the eigenvectors u_h
262  !! each eigenvector is a finite element function
263  do ii=1, nev
264  call gmsh_addview(x_h, eig%v(:,ii), &
265  & trim(shell_com), &
266  & 'Eigen-functions', real(ii, RP), ii)
267  end do
268  !!
269  !! shell command to visualise with gmsh
270  shell_com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'&
271  & //' '//trim(shell_com)
272  call system(shell_com)
273 
274  end if
275 
276 contains
277 
278  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
279  !!
282  subroutine pmv_m(y, x)
283  real(RP), dimension(:), intent(out) :: y
284  real(RP), dimension(:), intent(in) :: x
285 
286  call matvecprod(y, m, x)
287 
288  end subroutine pmv_m
289 
290  !! Solve K x = y
291  !!
292  subroutine kinv_pcg(y, ierr, x)
293  real(RP), dimension(:), intent(inout) :: y
294  logical , intent(out) :: ierr
295  real(RP), dimension(:), intent(in) :: x
296 
297  call solve(y, kry, x, k)
298 
299  ierr = kry%ierr
300 
301  end subroutine kinv_pcg
302 
303  !! Solve LU*x = y, LU = LDLT MUMPS decomposition
304  !! of K = (S1 + S2 - sigma*M2)
305  !!
306  subroutine kinv_mumps(x, ierr, y)
307  real(RP), dimension(:), intent(inout) :: x
308  logical , intent(out) :: ierr
309  real(RP), dimension(:), intent(in) :: y
310 
311  call solve(x, lu, y)
312  ierr = .false.
313  end subroutine kinv_mumps
314 
315  !! Compute and display the residuals on the
316  !! computed eigenvalues and eigenvectors
317  !!
318  subroutine residual()
319  real(RP), dimension(size) :: u, v, w
320  real(RP) :: res, nu
321  logical :: ierr
322 
323  write(*,*)
324  write(*,*) "Residual | Ku - lambda Mu |_l2 / | u |_l2"
325  !!
326  do ii=1, nev
327 
328  u = eig%v(:,ii)
329 
330  call matvecprod(v, k , u)
331  call pmv_m(w, u)
332 
333  w = v - eig%lambda(ii)*w
334  res = sqrt( sum(w*w) / sum(u*u) )
335 
336  print*, " i = ", int(ii,1), " : res = ", real(res, sp)
337 
338  end do
339 
340  write(*,*)
341  write(*,*) "Residual | \nu u - K^{-1} M u |_l2 / | u |_l2"
342  !!
343  do ii=1, nev
344  u = eig%v(:,ii)
345 
346  call pmv_m(v , u)
347  call kinv(w, ierr, v)
348 
349 
350  nu = 1.0_rp / ( eig%lambda(ii) - sigma )
351  w = w - nu*u
352  res = sqrt( sum(w*w) / sum(u*u) )
353 
354  print*, " i = ", int(ii,1), " : res = ", real(res, sp)
355 
356  end do
357 
358  end subroutine residual
359 
360 
361 end program eigen_laplacian_2d
362 
363 
integer, parameter eig_shift_invert
Shift-invert mode.
subroutine kinv_mumps(x, ierr, y)
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
Definition: choral.F90:166
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
Definition: choral.F90:27
integer, parameter fe_p2_2d
subroutine pmv_m(y, x)
Matrix-vector product x –> M*x.
real(rp) function e(x, v1, v2)
CHORAL CONSTANTS
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.
subroutine residual()
integer, parameter quad_gauss_trg_12
integer, parameter kry_cg
CG linear solver.
real(rp) function b(x, w1, w2)