Choral
diffusion_Dirichlet_Neumann.f90
Go to the documentation of this file.
1 
94 
96 
98  use choral
99 
100  implicit none
101 
102  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
103  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
104  !!
105  !! VARAIBLE DEFINITION : BEGIN
106  !!
107  !!
108 
109  !! verb = verbosity level
110  !!
111  integer , parameter :: verb = 2
112 
113  !! penal = penalty coeff for the Dirichlet condition
114  !! tol = linear system inversion tolerance
115  !!
116  real(RP), parameter :: penal = 1e10_rp
117  real(RP), parameter :: tol = 1e-15_rp
118 
119 
120  !! SPACE DISCRETISATION
121  !!
122  !! fe_v = finite element method (volumic)
123  !! fe_s = finite element method (surfacic)
124  !! qd_v = quadrature method (volumic)
125  !! qd_s = quadrature method (surfacic)
126  !! mesh_idx = index for the mesh
127  !!
128  integer , parameter :: fe_v = fe_p2_2d
129  integer , parameter :: fe_s = fe_p2_1d
130  integer , parameter :: qd_v = quad_gauss_trg_12
131  integer , parameter :: qd_s = quad_gauss_edg_4
132  character(LEN=1), parameter :: mesh_idx = "2"
133  !!
134  !! m = mesh
135  !! X_h = finite element space
136  !! qdm = integration method
137  !!
138  type(mesh) :: msh
139  type(fespace) :: X_h
140  type(quadmesh) :: qdm
141  !!
142  !! msh_file
143  character(len=100), parameter :: mesh_file= &
144  & trim(GMSH_DIR)//"square/square_"//mesh_idx//".msh"
145 
146 
147  !! LINEAR SYSTEM
148  !!
149  !!
150  !! kry = krylov method def.
151  !! mass = mass matrix
152  !! stiff = stiffness matrix
153  !! K = linear system matrix ( = mass + stiff )
154  !! rhs = right hand side
155  !! u_h = numerical solution
156  !!
157  type(krylov) :: kry
158  type(csr) :: mass, stiff, K
159  real(RP), dimension(:), allocatable :: rhs, u_h, aux
160 
161 
162  !! NUMERICAL ERRORS
163  !!
164  real(RP) :: err_L2, err_H1
165 
166  character(LEN=300) :: shell_com
167 
168  !!
169  !! VARAIBLE DEFINITION : END
170  !!
171  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
172  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173 
174 
175  !! !!!!!!!!!!!!!!!!!!!!! INITIALISATION
176  !!
177  call choral_init(verb=verb)
178  write(*,*)
179  write(*,*)'diffusion_Dirichlet_Neumann: start'
180  write(*,*)""
181  write(*,*)" EXAMPLE OF AN ELLIPTIC PROBLEM"
182  write(*,*)""
183  write(*,*)" -\Delta(u) + u = f on \Omega"
184  write(*,*)""
185  write(*,*)" WITH A MIXED DIRICHLET / NEUMANN BOUNDARY CONDITION "
186  write(*,*)""
187  write(*,*)""
188 
189  write(*,*) ""
190  write(*,*)"================================== MESH ASSEMBLING"
191  msh = mesh(mesh_file, 'gmsh')
192 
193  call print(msh)
194 
195  write(*,*) ""
196  write(*,*)"================================== FINITE ELEMENT&
197  & SPACE ASSEMBLING"
198  x_h =fespace(msh)
199  call set(x_h, fe_v)
200  call set(x_h, fe_s)
201  call assemble(x_h)
202  call print(x_h)
203 
204  write(*,*) ""
205  write(*,*)"================================== INTEGRATION METHOD&
206  & ON \Omega ASSEMBLING"
207  qdm = quadmesh(msh)
208  call set(qdm, qd_v)
209  call assemble(qdm)
210  call print(qdm)
211 
212  write(*,*) ""
213  write(*,*)"================================== MASS AND STIFFNESS&
214  & MATRIX ASSEMBLING"
215  call diffusion_massmat(mass , one_r3, x_h, qdm)
216  call diffusion_stiffmat(stiff, emetric, x_h, qdm)
217 
218 
219  write(*,*) ""
220  write(*,*)"================================== LINEAR SYSTEM MATRIX"
221  !!
222  !! K = mass + stiff
223  call add(k, stiff, 1._rp, mass, 1._rp)
224  !!
225  !! cleaning
226  call clear(mass)
227  call clear(stiff)
228 
229  write(*,*) ""
230  write(*,*)"================================== RIGHT HAND SIDE"
231  write(*,*)" 1- COMPUTATION OF &
232  & THE VOLUMIC SOURCE TERM F"
233  !!
234  !! \int_\Omega f v_i dx
235  call l2_product(rhs, f, x_h, qdm)
236 
237 
238  write(*,*) ""
239  write(*,*)"================================== RIGHT HAND SIDE"
240  write(*,*)" 2- COMPUTATION OF &
241  & THE SURFACIC SOURCE TERM G (NEUMANN CONDITION)"
242  !!
243  !! \int_{Gamma, x=1} dx_u_1 v_i dx
244  call diffusion_neumann_rhs(aux, dx_u_1, x_h, qd_s, x_ge_1)
245  rhs = rhs + aux
246  !!
247  !! \int_{Gamma, y=1} dy_u_1 v_i dx
248  call diffusion_neumann_rhs(aux, dy_u_1, x_h, qd_s, y_ge_1)
249  rhs = rhs + aux
250 
251  write(*,*) ""
252  write(*,*)"================================== RIGHT HAND SIDE"
253  write(*,*)" 3- PENALISATION OF &
254  & THE SYSTEM FOR THE DIRICHLET CONDITION"
255  !!
256  !! Dirichlet condition on {Gamma, x=0}
257  call diffusion_dirichlet(k, rhs, u, x_h, penal, x_le_0)
258  !!
259  !! Dirichlet condition on {Gamma, y=0}
260  call diffusion_dirichlet(k, rhs, u, x_h, penal, y_le_0)
261 
262  write(*,*) ""
263  write(*,*)"================================== SOLVE LINEAR SYSTEM&
264  & K U = RHS "
265  !!
266  !! Solver setting
267  kry = krylov(kry_cg, tol=tol, itmax=1000, verb=2)
268  call print(kry)
269  !!
270  !! initial guess
271  call interp_scal_func(u_h, u, x_h)
272  u_h = 0.0_rp
273  !!
274  !! linear system inversion
275  call solve(u_h, kry, rhs, k)
276  !!
277  !! cleaning
278  call clear(k)
279  call freemem(rhs)
280 
281  write(*,*) ""
282  write(*,*)"================================== POST-TREATMENT"
283  !!
284  !! graphical display
285  !!
286  !! output file def.
287  shell_com = 'diff_Neumann_u.msh'
288  !!
289  !! first write the 'finite element mesh'
290  call write(x_h, trim(shell_com), 'gmsh')
291  !!
292  !! second save the finite element solution u_h
293  call gmsh_addview(x_h, u_h, &
294  & trim(shell_com), &
295  & 'Diffusion problem: numerical solution', 0.0_rp, 1)
296  !!
297  !! shell command to visualise with gmsh
298  shell_com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'&
299  & //' '//trim(shell_com)
300  call system(shell_com)
301  !!
302  !! L2 and H1_0 numerical errors
303  write(*,*) ""
304  write(*,*) "Numerical errors"
305  err_l2 = l2_dist(u, u_h, x_h, qdm)
306  err_h1 = l2_dist_grad(grad_u, u_h, x_h, qdm)
307  write(*,*)" Error | u - u_h|_L2 =", err_l2
308  write(*,*)" Error |\nabla (u - u_h)|_L2 =", err_h1
309 
310  call freemem(u_h)
311  call freemem(aux)
312 
313 contains
314 
317  function u(x)
318  real(RP) :: u
319  real(RP), dimension(3), intent(in) :: x
320 
321  u = cos( x(1) ) * (1._rp + x(2)**2)
322 
323  end function u
324 
325 
328  function grad_u(x)
329  real(RP), dimension(3) :: grad_u
330  real(RP), dimension(3), intent(in) :: x
331 
332  grad_u(1) =-sin(x(1)) * (1._rp + x(2)**2)
333  grad_u(2) = cos(x(1)) * 2._rp * x(2)
334  grad_u(3) = 0._rp
335 
336  end function grad_u
337 
338 
341  function f(x)
342  real(RP) :: f
343  real(RP), dimension(3), intent(in) :: x
344 
345  f = 2._rp * u(x) - 2._rp * cos( x(1) )
346 
347 
348  end function f
349 
350 
353  function x_le_0(x) result(r)
354  real(RP) :: r
355  real(RP), dimension(3), intent(in) :: x
356 
357  r = -x(1)
358 
359  end function x_le_0
360 
361 
364  function x_ge_1(x) result(r)
365  real(RP) :: r
366  real(RP), dimension(3), intent(in) :: x
367 
368  r = x(1) - 1.0_rp
369 
370  end function x_ge_1
371 
374  function y_le_0(x) result(r)
375  real(RP) :: r
376  real(RP), dimension(3), intent(in) :: x
377 
378  r = -x(2)
379 
380  end function y_le_0
381 
382 
385  function y_ge_1(x) result(r)
386  real(RP) :: r
387  real(RP), dimension(3), intent(in) :: x
388 
389  r = x(2) - 1.0_rp
390 
391  end function y_ge_1
392 
393 
396  function dx_u_1(x) result(r)
397  real(RP) :: r
398  real(RP), dimension(3), intent(in) :: x
399 
400  r = -sin( 1.0_rp ) * ( 1.0_rp + x(2)**2 )
401 
402  end function dx_u_1
403 
406  function dy_u_1(x) result(r)
407  real(RP) :: r
408  real(RP), dimension(3), intent(in) :: x
409 
410  r = cos( x(1) ) * 2.0_rp
411 
412  end function dy_u_1
413 
414 end program diffusion_dirichlet_neumann
real(rp) function x_ge_1(x)
To characterise {x=1}.
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
integer, parameter quad_gauss_edg_4
real(rp) function x_le_0(x)
To characterise {x=0}.
real(rp) function e(x, v1, v2)
real(rp) function y_le_0(x)
To characterise {x=0}.
real(rp) function, dimension(3) grad_u(x)
gradient of the exact solution
CHORAL CONSTANTS
real(rp) function f(x)
right hand side 'f'
real(rp) function u(x)
exact solution 'u'
real(rp) function y_ge_1(x)
To characterise {x=1}.
integer, parameter quad_gauss_trg_12
real(rp) function dx_u_1(x)
u on {x=1}
program diffusion_dirichlet_neumann
SOLVES THE DIFFUSION PROBLEM with mixed Dirichlet / Neumann boundary conditions ...
integer, parameter fe_p2_1d
integer, parameter kry_cg
CG linear solver.
real(rp) function dy_u_1(x)
u on {y=1}