All Classes Files Functions Variables Pages
module_init_medium.f90
Go to the documentation of this file.
1 !=====================================================================================================================================
2 ! EFISPEC3D !
3 ! (Elements FInis SPECtraux 3D) !
4 ! !
5 ! http://efispec.free.fr !
6 ! !
7 ! !
8 ! This file is part of EFISPEC3D !
9 ! Please refer to http://efispec.free.fr if you use it or part of it !
10 ! !
11 ! !
12 !1 ---> French License: CeCILL V2 !
13 ! !
14 ! Copyright BRGM 2009 contributeurs : Florent DE MARTIN !
15 ! David MICHEA !
16 ! Philippe THIERRY !
17 ! !
18 ! Contact: f.demartin at brgm.fr !
19 ! !
20 ! Ce logiciel est un programme informatique servant a resoudre l'equation du !
21 ! mouvement en trois dimensions via une methode des elements finis spectraux. !
22 ! !
23 ! Ce logiciel est regi par la licence CeCILL soumise au droit francais et !
24 ! respectant les principes de diffusion des logiciels libres. Vous pouvez !
25 ! utiliser, modifier et/ou redistribuer ce programme sous les conditions de la !
26 ! licence CeCILL telle que diffusee par le CEA, le CNRS et l'INRIA sur le site !
27 ! "http://www.cecill.info". !
28 ! !
29 ! En contrepartie de l'accessibilite au code source et des droits de copie, de !
30 ! modification et de redistribution accordes par cette licence, il n'est offert !
31 ! aux utilisateurs qu'une garantie limitee. Pour les memes raisons, seule une !
32 ! responsabilite restreinte pese sur l'auteur du programme, le titulaire des !
33 ! droits patrimoniaux et les concedants successifs. !
34 ! !
35 ! A cet egard l'attention de l'utilisateur est attiree sur les risques associes !
36 ! au chargement, a l'utilisation, a la modification et/ou au developpement et a !
37 ! la reproduction du logiciel par l'utilisateur etant donne sa specificite de !
38 ! logiciel libre, qui peut le rendre complexe a manipuler et qui le reserve donc !
39 ! a des developpeurs et des professionnels avertis possedant des connaissances !
40 ! informatiques approfondies. Les utilisateurs sont donc invites a charger et !
41 ! tester l'adequation du logiciel a leurs besoins dans des conditions permettant !
42 ! d'assurer la securite de leurs systemes et ou de leurs donnees et, plus !
43 ! generalement, a l'utiliser et l'exploiter dans les memes conditions de !
44 ! securite. !
45 ! !
46 ! Le fait que vous puissiez acceder a cet en-tete signifie que vous avez pris !
47 ! connaissance de la licence CeCILL et que vous en avez accepte les termes. !
48 ! !
49 ! !
50 !2 ---> International license: GNU GPL V3 !
51 ! !
52 ! EFISPEC3D is a computer program that solves the three-dimensional equations of !
53 ! motion using a finite spectral-element method. !
54 ! !
55 ! Copyright (C) 2009 Florent DE MARTIN !
56 ! !
57 ! Contact: f.demartin at brgm.fr !
58 ! !
59 ! This program is free software: you can redistribute it and/or modify it under !
60 ! the terms of the GNU General Public License as published by the Free Software !
61 ! Foundation, either version 3 of the License, or (at your option) any later !
62 ! version. !
63 ! !
64 ! This program is distributed in the hope that it will be useful, but WITHOUT ANY !
65 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A !
66 ! PARTICULAR PURPOSE. See the GNU General Public License for more details. !
67 ! !
68 ! You should have received a copy of the GNU General Public License along with !
69 ! this program. If not, see http://www.gnu.org/licenses/. !
70 ! !
71 ! !
72 !3 ---> Third party libraries !
73 ! !
74 ! EFISPEC3D uses the following of third party libraries: !
75 ! !
76 ! --> METIS 5.1.0 !
77 ! see http://glaros.dtc.umn.edu/gkhome/metis/metis/overview !
78 ! !
79 ! --> Lib_VTK_IO !
80 ! see S. Zaghi's website: https://github.com/szaghi/Lib_VTK_IO !
81 ! !
82 ! --> INTERP_LINEAR !
83 ! see J. Burkardt website: http://people.sc.fsu.edu/~jburkardt/ !
84 ! !
85 ! --> SAC !
86 ! http://ds.iris.edu/files/sac-manual/ !
87 ! !
88 ! --> EXODUS II !
89 ! http://sourceforge.net/projects/exodusii/ !
90 ! !
91 ! --> NETCDF !
92 ! http://www.unidata.ucar.edu/software/netcdf/ !
93 ! !
94 ! --> HDF5 !
95 ! http://www.hdfgroup.org/HDF5/ !
96 ! !
97 ! Some of these libraries are located in directory lib and pre-compiled !
98 ! with intel compiler for x86_64 platform. !
99 ! !
100 ! !
101 !4 ---> Related Articles !
102 ! !
103 ! De Martin, F., Matsushima, M., Kawase, H. (BSSA, 2013) !
104 ! Impact of geometric effects on near-surface Green's functions !
105 ! doi:10.1785/0120130039 !
106 ! !
107 ! Aochi, H., Ducellier, A., Dupros, F., Delatre, M., Ulrich, T., De Martin, F., !
108 ! Yoshimi, M., (Pure Appl. Geophys. 2013) !
109 ! Finite Difference Simulations of Seismic Wave Propagation for the 2007 Mw 6.6 !
110 ! Niigata-ken Chuetsu-Oki Earthquake: Validity of Models and Reliable !
111 ! Input Ground Motion in the Near-Field !
112 ! doi:10.1007/s00024-011-0429-5 !
113 ! !
114 ! De Martin, F. (BSSA, 2011) !
115 ! Verification of a Spectral-Element Method Code for the Southern California !
116 ! Earthquake Center LOH.3 Viscoelastic Case !
117 ! doi:10.1785/0120100305 !
118 ! !
119 !=====================================================================================================================================
120 
123 
126 module mod_init_medium
127 
128  use mpi
129 
130  implicit none
131 
132  private
133 
134  public :: init_hexa_medium
135  public :: init_quadp_medium
136  private :: init_medium_from_mat_file
137  private :: write_medium_vtk_cell_xml
138  private :: write_medium_vtk_node_xml
139  private :: write_medium_collection
140 
141  contains
142 
143 
144 !
145 !
157 !***********************************************************************************************************************************************************************************
158  subroutine init_hexa_medium()
159 !***********************************************************************************************************************************************************************************
160 
161  use mpi
162 
163  use mod_global_variables, only :&
164  ig_nmaterial&
165  ,tg_elastic_material&
166  ,ig_myrank&
167  ,cg_prefix&
168  ,cg_myrank&
169  ,ig_ngll&
170  ,ig_material_type&
171  ,ig_hexa_material_number&
172  ,lg_visco&
173  ,lg_output_medium_vtk&
174  ,ig_nhexa&
175  ,rg_hexa_gll_rho&
176  ,rg_hexa_gll_rhovs2&
177  ,rg_hexa_gll_rhovp2&
178  ,rg_hexa_gll_wkqs&
179  ,rg_hexa_gll_wkqp&
180  ,rg_hexa_gll_ksixx&
181  ,rg_hexa_gll_ksiyy&
182  ,rg_hexa_gll_ksizz&
183  ,rg_hexa_gll_ksixy&
184  ,rg_hexa_gll_ksixz&
185  ,rg_hexa_gll_ksiyz&
186  ,rg_mem_var_exp&
187  ,tg_viscoelastic_material&
188  ,rg_relax_coeff&
189  ,ig_nrelax&
190  ,rg_pi&
191  ,rg_dt&
192  ,error_stop&
193  ,ig_lst_unit
194 
195  use mod_init_memory
196 
197  implicit none
198 
199  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL,ig_nhexa) :: var1
200  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL,ig_nhexa) :: var2
201 
202  complex :: ctmpqs
203  complex :: ctmpqp
204 
205  real :: dtmpqs
206  real :: dtmpqp
207 
208  integer :: ios
209  integer :: iel
210  integer :: imat
211  integer :: imem_var
212  integer :: k
213  integer :: l
214  integer :: m
215  integer :: n
216 
217  logical :: is_file_exist
218 
219  character(len=255) :: info
220 
221 
222  !
223  !
224  !************************************************************************************************
225  !initialize memory
226  !************************************************************************************************
227 
228  ios = init_array_real(rg_hexa_gll_rho,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_rho") !deallocated in subroutine init_mass_matrix.f90
229 
230  ios = init_array_real(rg_hexa_gll_rhovs2,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_rhovs2")
231 
232  ios = init_array_real(rg_hexa_gll_rhovp2,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_rhovp2")
233 
234 
235  if (lg_visco) then
236 
237  ios = init_array_real(rg_hexa_gll_ksixx,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksixx")
238 
239  ios = init_array_real(rg_hexa_gll_ksiyy,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksiyy")
240 
241  ios = init_array_real(rg_hexa_gll_ksizz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksizz")
242 
243  ios = init_array_real(rg_hexa_gll_ksixy,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksixy")
244 
245  ios = init_array_real(rg_hexa_gll_ksixz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksixz")
246 
247  ios = init_array_real(rg_hexa_gll_ksiyz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_ksiyz")
248 
249  ios = init_array_real(rg_hexa_gll_wkqs,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_wkqs")
250 
251  ios = init_array_real(rg_hexa_gll_wkqp,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,ig_nrelax,"rg_hexa_gll_wkqp")
252 
253 !
254 !------->precompute exponential used for computing internal forces
255  do imem_var = 1,ig_nrelax
256  rg_mem_var_exp(imem_var) = exp(-rg_dt/rg_relax_coeff(imem_var,1))
257  enddo
258 
259  endif
260 
261 
262  !
263  !
264  !***********************************************************************************************************
265  !->if files *.mat exists, then the medium (ig_hexa_material_number) is read from such files
266  !***********************************************************************************************************
267  is_file_exist = .false.
268 
269  inquire(file=trim(cg_prefix)//".cpu."//trim(cg_myrank)//".mat", exist=is_file_exist)
270  if (is_file_exist) then
271  if (ig_myrank == 0) then
272  write(ig_lst_unit,'(a)') " --> medium generated from geomodeler's material (file *.mat)"
273  endif
274  call init_medium_from_mat_file(ig_hexa_material_number)
275  endif
276 
277  if ( .not.(is_file_exist) .and. (ig_myrank == 0) ) then
278  write(ig_lst_unit,'(a)') " --> medium generated from cubit block's information"
279  endif
280 
281 
282  !
283  !
284  !************************************************************************************************
285  !->check if the medium (ig_hexa_material_number) is entirely defined
286  !************************************************************************************************
287  do iel = 1,ig_nhexa
288 
289  imat = ig_hexa_material_number(iel)
290 
291  if ( (imat <= 0) .or. (imat > ig_nmaterial) ) then
292  write(info,'(a)') "error in subroutine init_hexa_medium: undefined material for hexa"
293  call error_stop(info)
294  endif
295 
296  enddo
297 
298 
299  !
300  !
301  !************************************************************************************************
302  !Compute elastic and viscoelastic parameter pour media_type = 0
303  !************************************************************************************************
304  do imat = 1,ig_nmaterial
305 !
306 !
307 !------->elastic
308  tg_elastic_material(imat)%lam2 = tg_elastic_material(imat)%dens*tg_elastic_material(imat)%svel**2
309  tg_elastic_material(imat)%lam1 = tg_elastic_material(imat)%dens*tg_elastic_material(imat)%pvel**2 - 2.0*tg_elastic_material(imat)%lam2
310  tg_elastic_material(imat)%pois = (tg_elastic_material(imat)%pvel**2 - 2.0*tg_elastic_material(imat)%svel**2)/(2.0*(tg_elastic_material(imat)%pvel**2-tg_elastic_material(imat)%svel**2))
311  tg_elastic_material(imat)%youn = 2.0*tg_elastic_material(imat)%lam2*(1.0+tg_elastic_material(imat)%pois)
312  tg_elastic_material(imat)%bulk = tg_elastic_material(imat)%lam1+2.0/3.0*tg_elastic_material(imat)%lam2
313  tg_elastic_material(imat)%rvel = (0.862+1.14*tg_elastic_material(imat)%pois)/(1.0+tg_elastic_material(imat)%pois)*tg_elastic_material(imat)%svel
314  tg_elastic_material(imat)%pwmo = tg_elastic_material(imat)%lam1 + 2.0*tg_elastic_material(imat)%lam2
315 
316  if (ig_myrank == 0) then
317  if (ig_material_type(imat) == 1) write(ig_lst_unit,'(" ",/,a,i0,a)') "material ",imat," is elastic"
318  if (ig_material_type(imat) == 2) write(ig_lst_unit,'(" ",/,a,i0,a)') "material ",imat," is viscoelactic"
319  write(ig_lst_unit,'(a,f15.8 ,a)') "s-wave velocity = ",tg_elastic_material(imat)%svel," m/s"
320  write(ig_lst_unit,'(a,f15.8 ,a)') "p-wave velocity = ",tg_elastic_material(imat)%pvel," m/s"
321  write(ig_lst_unit,'(a,f15.8 ,a)') "poisson ratio = ",tg_elastic_material(imat)%pois," "
322  write(ig_lst_unit,'(a,f15.8 ,a)') "density = ",tg_elastic_material(imat)%dens," kg/m3"
323  write(ig_lst_unit,'(a,e15.7e3,a)') "lambda (first lame coef.) = ",tg_elastic_material(imat)%lam1," Pa"
324  write(ig_lst_unit,'(a,e15.7e3,a)') "s-wave modulus = ",tg_elastic_material(imat)%lam2," Pa"
325  write(ig_lst_unit,'(a,e15.7e3,a)') "p-wave modulus = ",tg_elastic_material(imat)%pwmo," Pa"
326  write(ig_lst_unit,'(a,e15.7e3,a)') "bulk-modulus = ",tg_elastic_material(imat)%bulk," Pa"
327  endif
328 !
329 !
330 !------->viscoelastic
331  if (lg_visco) then
332 
333  allocate(tg_viscoelastic_material(imat)%wkqs(ig_nrelax),tg_viscoelastic_material(imat)%wkqp(ig_nrelax))
334 
335  dtmpqs = (3.071+1.433*tg_viscoelastic_material(imat)%qs**(-1.158)*log(tg_viscoelastic_material(imat)%qs/5.0))/(1.0+0.415*tg_viscoelastic_material(imat)%qs)
336  dtmpqp = (3.071+1.433*tg_viscoelastic_material(imat)%qp**(-1.158)*log(tg_viscoelastic_material(imat)%qp/5.0))/(1.0+0.415*tg_viscoelastic_material(imat)%qp)
337  do k = 1,ig_nrelax
338  tg_viscoelastic_material(imat)%wkqs(k) = dtmpqs*(dtmpqs*rg_relax_coeff(k,2)+rg_relax_coeff(k,3))
339  tg_viscoelastic_material(imat)%wkqp(k) = dtmpqp*(dtmpqp*rg_relax_coeff(k,2)+rg_relax_coeff(k,3))
340  enddo
341 
342  ctmpqs = cmplx(0.0,0.0)
343  ctmpqp = cmplx(0.0,0.0)
344  do k = 1,ig_nrelax
345  ctmpqs = ctmpqs + cmplx(tg_viscoelastic_material(imat)%wkqs(k),0.0)/(1.0 + cmplx(0.0,1.0)*2.0*rg_pi*tg_viscoelastic_material(imat)%freq*rg_relax_coeff(k,1))
346  ctmpqp = ctmpqp + cmplx(tg_viscoelastic_material(imat)%wkqp(k),0.0)/(1.0 + cmplx(0.0,1.0)*2.0*rg_pi*tg_viscoelastic_material(imat)%freq*rg_relax_coeff(k,1))
347  enddo
348 
349  tg_viscoelastic_material(imat)%uswm = tg_elastic_material(imat)%lam2/abs(1.0-ctmpqs)
350  tg_viscoelastic_material(imat)%upwm = tg_elastic_material(imat)%pwmo/abs(1.0-ctmpqp)
351  if (ig_myrank == 0) then
352  write(ig_lst_unit,'(a,f15.8 ,a)') "frequency for viscosity = ",tg_viscoelastic_material(imat)%freq," hz"
353  write(ig_lst_unit,'(a,f15.8 ,a)') "s-wave quality factor = ",tg_viscoelastic_material(imat)%qs ," "
354  write(ig_lst_unit,'(a,f15.8 ,a)') "p-wave quality factor = ",tg_viscoelastic_material(imat)%qp ," "
355  write(ig_lst_unit,'(a,e15.7e3,a)') "unrelaxed s-wave modulus = ",tg_viscoelastic_material(imat)%uswm," pa"
356  write(ig_lst_unit,'(a,e15.7e3,a)') "unrelaxed p-wave modulus = ",tg_viscoelastic_material(imat)%upwm," pa"
357  endif
358 
359  endif !endif of material type viscoelastic
360 
361  enddo !enddo of loop on material
362 
363 
364  !
365  !
366  !************************************************************************************************
367  !Distribute material properties to gll nodes
368  !************************************************************************************************
369 !
370 !
371 !------->elastic
372  if (.not.lg_visco) then
373  do iel = 1,ig_nhexa
374  do k = 1,ig_ngll !zeta
375  do l = 1,ig_ngll !eta
376  do m = 1,ig_ngll !xi
377 
378  imat = ig_hexa_material_number(iel)
379  rg_hexa_gll_rho(m,l,k,iel) = tg_elastic_material(imat)%dens
380  rg_hexa_gll_rhovs2(m,l,k,iel) = tg_elastic_material(imat)%dens*tg_elastic_material(imat)%svel**2
381  rg_hexa_gll_rhovp2(m,l,k,iel) = tg_elastic_material(imat)%dens*tg_elastic_material(imat)%pvel**2
382 
383  enddo
384  enddo
385  enddo
386  enddo
387 !
388 !
389 !------->viscoelastic
390  else
391  do iel = 1,ig_nhexa
392  do k = 1,ig_ngll !zeta
393  do l = 1,ig_ngll !eta
394  do m = 1,ig_ngll !xi
395 
396  imat = ig_hexa_material_number(iel)
397  rg_hexa_gll_rho(m,l,k,iel) = tg_elastic_material(imat)%dens
398  rg_hexa_gll_rhovs2(m,l,k,iel) = tg_viscoelastic_material(imat)%uswm
399  rg_hexa_gll_rhovp2(m,l,k,iel) = tg_viscoelastic_material(imat)%upwm
400 
401  do n = 1,ig_nrelax
402  rg_hexa_gll_wkqs(n,m,l,k,iel) = tg_viscoelastic_material(imat)%wkqs(n)
403  rg_hexa_gll_wkqp(n,m,l,k,iel) = tg_viscoelastic_material(imat)%wkqp(n)
404  enddo
405 
406  enddo
407  enddo
408  enddo
409  enddo
410  endif
411 
412  !
413  !
414  !************************************************************************************************
415  !write material in VTK xml format
416  !************************************************************************************************
417 
418 ! if (LG_OUTPUT_MEDIUM_VTK) call write_medium_vtk_cell_xml()
419  if (lg_output_medium_vtk) call write_medium_vtk_node_xml(var1,var2)
420 
421  return
422 !***********************************************************************************************************************************************************************************
423  end subroutine init_hexa_medium
424 !***********************************************************************************************************************************************************************************
425 
426 !
427 !
432 !***********************************************************************************************************************************************************************************
433  subroutine init_quadp_medium()
434 !***********************************************************************************************************************************************************************************
435  use mpi
436  use mod_global_variables, only :&
437  ig_ngll&
438  ,ig_myrank&
439  ,ig_nquad_parax&
440  ,ig_quadp_neighbor_hexa&
441  ,ig_quadp_neighbor_hexaface&
442  ,rg_quadp_gll_rhovs&
443  ,rg_quadp_gll_rhovp&
444  ,rg_hexa_gll_rho&
445  ,rg_hexa_gll_rhovs2&
446  ,rg_hexa_gll_rhovp2&
447  ,error_stop&
448  ,lg_boundary_absorption
449 
450  use mod_init_memory
451 
452  implicit none
453 
454  integer :: ios
455  integer :: i
456  integer :: j
457  integer :: k
458  integer :: iquad
459  integer :: ihexa
460  integer :: iface
461 
462  character(len=255) :: info
463 
464  !
465  !See numbering convention in subroutine propagate_gll_nodes_quad of module_mesh_gll_num.f90
466  !
467 
468  if (ig_nquad_parax > 0 .and. lg_boundary_absorption) then
469 
470  ios = init_array_real(rg_quadp_gll_rhovs,ig_nquad_parax,ig_ngll,ig_ngll,"rg_quadp_gll_rhovs")
471 
472  ios = init_array_real(rg_quadp_gll_rhovp,ig_nquad_parax,ig_ngll,ig_ngll,"rg_quadp_gll_rhovp")
473 
474  do iquad = 1,ig_nquad_parax
475 
476  ihexa = ig_quadp_neighbor_hexa(iquad)
477  iface = ig_quadp_neighbor_hexaface(iquad)
478 
479  select case(iface)
480  case(1)
481  do j = 1,ig_ngll
482  do i = 1,ig_ngll
483  rg_quadp_gll_rhovs(i,j,iquad) = rg_hexa_gll_rho(i,j,1,ihexa)*sqrt(rg_hexa_gll_rhovs2(i,j,1,ihexa)/rg_hexa_gll_rho(i,j,1,ihexa))
484  rg_quadp_gll_rhovp(i,j,iquad) = rg_hexa_gll_rho(i,j,1,ihexa)*sqrt(rg_hexa_gll_rhovp2(i,j,1,ihexa)/rg_hexa_gll_rho(i,j,1,ihexa))
485  enddo
486  enddo
487 
488  case(2)
489  do k = 1,ig_ngll
490  do i = 1,ig_ngll
491  rg_quadp_gll_rhovs(k,i,iquad) = rg_hexa_gll_rho(i,1,k,ihexa)*sqrt(rg_hexa_gll_rhovs2(i,1,k,ihexa)/rg_hexa_gll_rho(i,1,k,ihexa))
492  rg_quadp_gll_rhovp(k,i,iquad) = rg_hexa_gll_rho(i,1,k,ihexa)*sqrt(rg_hexa_gll_rhovp2(i,1,k,ihexa)/rg_hexa_gll_rho(i,1,k,ihexa))
493  enddo
494  enddo
495 
496  case(3)
497  do k = 1,ig_ngll
498  do j = 1,ig_ngll
499  rg_quadp_gll_rhovs(k,j,iquad) = rg_hexa_gll_rho(ig_ngll,j,k,ihexa)*sqrt(rg_hexa_gll_rhovs2(ig_ngll,j,k,ihexa)/rg_hexa_gll_rho(ig_ngll,j,k,ihexa))
500  rg_quadp_gll_rhovp(k,j,iquad) = rg_hexa_gll_rho(ig_ngll,j,k,ihexa)*sqrt(rg_hexa_gll_rhovp2(ig_ngll,j,k,ihexa)/rg_hexa_gll_rho(ig_ngll,j,k,ihexa))
501  enddo
502  enddo
503 
504  case(4)
505  do k = 1,ig_ngll
506  do i = 1,ig_ngll
507  rg_quadp_gll_rhovs(i,k,iquad) = rg_hexa_gll_rho(i,ig_ngll,k,ihexa)*sqrt(rg_hexa_gll_rhovs2(i,ig_ngll,k,ihexa)/rg_hexa_gll_rho(i,ig_ngll,k,ihexa))
508  rg_quadp_gll_rhovp(i,k,iquad) = rg_hexa_gll_rho(i,ig_ngll,k,ihexa)*sqrt(rg_hexa_gll_rhovp2(i,ig_ngll,k,ihexa)/rg_hexa_gll_rho(i,ig_ngll,k,ihexa))
509  enddo
510  enddo
511 
512  case(5)
513  do k = 1,ig_ngll
514  do j = 1,ig_ngll
515  rg_quadp_gll_rhovs(j,k,iquad) = rg_hexa_gll_rho(1,j,k,ihexa)*sqrt(rg_hexa_gll_rhovs2(1,j,k,ihexa)/rg_hexa_gll_rho(1,j,k,ihexa))
516  rg_quadp_gll_rhovp(j,k,iquad) = rg_hexa_gll_rho(1,j,k,ihexa)*sqrt(rg_hexa_gll_rhovp2(1,j,k,ihexa)/rg_hexa_gll_rho(1,j,k,ihexa))
517  enddo
518  enddo
519 
520  case(6)
521  do j = 1,ig_ngll
522  do i = 1,ig_ngll
523  rg_quadp_gll_rhovs(j,i,iquad) = rg_hexa_gll_rho(i,j,ig_ngll,ihexa)*sqrt(rg_hexa_gll_rhovs2(i,j,ig_ngll,ihexa)/rg_hexa_gll_rho(i,j,ig_ngll,ihexa))
524  rg_quadp_gll_rhovp(j,i,iquad) = rg_hexa_gll_rho(i,j,ig_ngll,ihexa)*sqrt(rg_hexa_gll_rhovp2(i,j,ig_ngll,ihexa)/rg_hexa_gll_rho(i,j,ig_ngll,ihexa))
525  enddo
526  enddo
527 
528  case default
529  write(info,'(a)') "error in subroutine init_quadp_medium: invalid face number"
530  call error_stop(info)
531 
532  end select
533 
534  enddo
535 
536  endif
537 
538  if (allocated(ig_quadp_neighbor_hexa )) deallocate(ig_quadp_neighbor_hexa)
539  if (allocated(ig_quadp_neighbor_hexaface)) deallocate(ig_quadp_neighbor_hexaface)
540 
541  return
542 !***********************************************************************************************************************************************************************************
543  end subroutine init_quadp_medium
544 !***********************************************************************************************************************************************************************************
545 
546 !
547 !
551 !***********************************************************************************************************************************************************************************
552  subroutine init_medium_from_mat_file(il_mat_num_of_hexa)
553 !***********************************************************************************************************************************************************************************
554 
555  use mpi
556 
557  use mod_global_variables, only :&
558  ig_nhexa&
559  ,cg_myrank&
560  ,cg_prefix
561 
562  implicit none
563 
564  integer, intent(out) :: il_mat_num_of_hexa(:)
565  integer :: ihexa
566  !
567  !->fill the array 'il_mat_num_of_hexa' using the file *.mat specific for myrank
568  open(unit=90,file=trim(cg_prefix)//".cpu."//trim(cg_myrank)//".mat")
569  do ihexa = 1,ig_nhexa
570  read(unit=90,fmt=*) il_mat_num_of_hexa(ihexa)
571  enddo
572  close(90)
573 
574  return
575 !***********************************************************************************************************************************************************************************
576  end subroutine init_medium_from_mat_file
577 !***********************************************************************************************************************************************************************************
578 
579 !
580 !
584 !***********************************************************************************************************************************************************************************
585  subroutine write_medium_vtk_cell_xml()
586 !***********************************************************************************************************************************************************************************
587 
588  use mpi
589 
590  use mod_global_variables, only :&
591  ig_hexa_material_number&
592  ,ig_nhexa&
593  ,ig_hexa_nnode&
594  ,ig_hexa_gnode_glonum&
595  ,ig_mesh_nnode&
596  ,rg_gnode_x&
597  ,rg_gnode_y&
598  ,rg_gnode_z&
599  ,tg_elastic_material&
600  ,cg_prefix&
601  ,ig_myrank&
602  ,cg_myrank
603 
604  use mod_vtk_io
605 
606  implicit none
607 
608  real , dimension(ig_nhexa) :: var
609 
610  integer(kind=1), dimension(ig_nhexa) :: cell_type
611  integer , dimension(ig_nhexa) :: offset
612  integer , dimension(ig_hexa_nnode*ig_nhexa) :: connectivity
613  integer :: ihexa
614  integer :: inode
615  integer :: imat
616  integer :: icon
617  integer :: ios
618 
619  character(len=255) :: fname
620 
621 !
622 !---->file name to store the medium of each cpu
623  fname = trim(cg_prefix)//'.medium.cpu.'//trim(cg_myrank)//".vtu"
624 
625 !
626 !---->fill the array cell_type with 12 --> hexahedron 8 nodes
627  do ihexa = 1,ig_nhexa
628  cell_type(ihexa) = 12
629  enddo
630 
631 !
632 !---->fill the array offset : sum of vertices over all hexahedra
633  do ihexa = 1,ig_nhexa
634  offset(ihexa) = ihexa*ig_hexa_nnode
635  enddo
636 
637 !
638 !---->fill the array connectivity
639  icon = 0
640  do ihexa = 1,ig_nhexa
641  do inode = 1,ig_hexa_nnode
642 
643  icon = icon + 1
644  connectivity(icon) = ig_hexa_gnode_glonum(inode,ihexa) - 1
645 
646  enddo
647  enddo
648 
649  ios = vtk_ini_xml( &
650  output_format = 'BINARY' &
651  ,filename = fname &
652  ,mesh_topology = 'UnstructuredGrid' )
653 
654  ios = vtk_geo_xml( &
655  nn = ig_mesh_nnode &
656  ,nc = ig_nhexa &
657  ,x = rg_gnode_x &
658  ,y = rg_gnode_y &
659  ,z = rg_gnode_z )
660 
661  ios = vtk_con_xml( &
662  nc = ig_nhexa &
663  ,connect = connectivity &
664  ,offset = offset &
665  ,cell_type = cell_type )
666 
667  ios = vtk_dat_xml( &
668  var_location = 'CELL' &
669  ,var_block_action = 'OPEN' )
670 
671 !
672 !---->fill the array var with P-wave velocity value
673  do ihexa = 1,ig_nhexa
674 
675  imat = ig_hexa_material_number(ihexa)
676  var(ihexa) = tg_elastic_material(imat)%pvel
677 
678  enddo
679 
680  ios = vtk_var_xml( &
681  nc_nn = ig_nhexa &
682  ,varname = 'Vp' &
683  ,var = var &
684  )
685 !
686 !---->fill the array var with S-wave velocity value
687  do ihexa = 1,ig_nhexa
688  imat = ig_hexa_material_number(ihexa)
689  var(ihexa) = tg_elastic_material(imat)%svel
690  enddo
691 
692  ios = vtk_var_xml( &
693  nc_nn = ig_nhexa &
694  ,varname = 'Vs' &
695  ,var = var &
696  )
697 
698  ios = vtk_dat_xml( &
699  var_location = 'CELL' &
700  ,var_block_action = 'CLOSE' &
701  )
702 
703  ios = vtk_geo_xml()
704 
705  ios = vtk_end_xml()
706 
707  if (ig_myrank == 0) call write_medium_collection()
708 
709  return
710 !***********************************************************************************************************************************************************************************
711  end subroutine write_medium_vtk_cell_xml
712 !***********************************************************************************************************************************************************************************
713 
714 !
715 !
719 !***********************************************************************************************************************************************************************************
720  subroutine write_medium_vtk_node_xml(var1,var2)
721 !***********************************************************************************************************************************************************************************
722 
723  use mpi
724 
725  use mod_global_variables, only :&
726  ig_hexa_material_number&
727  ,ig_ngll&
728  ,ig_nhexa&
729  ,ig_hexa_nnode&
730  ,ig_hexa_gnode_glonum&
731  ,ig_mesh_nnode&
732  ,ig_hexa_node2gll&
733  ,rg_gnode_x&
734  ,rg_gnode_y&
735  ,rg_gnode_z&
736  ,rg_hexa_gll_rho&
737  ,rg_hexa_gll_rhovs2&
738  ,cg_prefix&
739  ,ig_myrank&
740  ,cg_myrank
741 
742  use mod_vtk_io
743 
744  implicit none
745 
746  real, intent(in), dimension(IG_NGLL,IG_NGLL,IG_NGLL,ig_nhexa) :: var1
747  real, intent(in), dimension(IG_NGLL,IG_NGLL,IG_NGLL,ig_nhexa) :: var2
748 
749  real , dimension(ig_mesh_nnode) :: var
750  real :: myvar
751 
752  integer(kind=1), dimension(ig_nhexa) :: cell_type
753  integer , dimension(ig_nhexa) :: offset
754  integer , dimension(ig_hexa_nnode*ig_nhexa) :: connectivity
755  integer :: ihexa
756  integer :: inode
757  integer :: global_gnode
758  integer :: icon
759  integer :: ios
760  integer :: k
761  integer :: l
762  integer :: m
763 
764  character(len=255) :: fname
765 
766 !
767 !---->file name to store the medium of each cpu
768  fname = trim(cg_prefix)//'.medium.cpu.'//trim(cg_myrank)//".vtu"
769 
770 !
771 !---->fill the array cell_type with 12 --> hexahedron 8 nodes
772  do ihexa = 1,ig_nhexa
773  cell_type(ihexa) = 12
774  enddo
775 
776 !
777 !---->fill the array offset : sum of vertices over all hexahedra
778  do ihexa = 1,ig_nhexa
779  offset(ihexa) = ihexa*ig_hexa_nnode
780  enddo
781 
782 !
783 !---->fill the array connectivity
784  icon = 0
785  do ihexa = 1,ig_nhexa
786  do inode = 1,ig_hexa_nnode
787 
788  icon = icon + 1
789  connectivity(icon) = ig_hexa_gnode_glonum(inode,ihexa) - 1
790 
791  enddo
792  enddo
793 
794 !
795 !---->initialize VTK xml file
796  ios = vtk_ini_xml( &
797  output_format = 'BINARY' &
798  ,filename = fname &
799  ,mesh_topology = 'UnstructuredGrid' )
800 
801 !
802 !---->geometry information for VTK xml file
803  ios = vtk_geo_xml( &
804  nn = ig_mesh_nnode &
805  ,nc = ig_nhexa &
806  ,x = rg_gnode_x &
807  ,y = rg_gnode_y &
808  ,z = rg_gnode_z &
809  )
810 
811 !
812 !---->connectivity information for VTK xml file
813  ios = vtk_con_xml( &
814  nc = ig_nhexa &
815  ,connect = connectivity &
816  ,offset = offset &
817  ,cell_type = cell_type &
818  )
819 
820 !
821 !---->data storage type for VTK xml file : cell or node
822  ios = vtk_dat_xml( &
823  var_location = 'NODE' &
824  ,var_block_action = 'OPEN' &
825  )
826 
827 !
828 !---->fill the array var with gll_node value
829 
830  do ihexa = 1,ig_nhexa
831 
832  do inode = 1,ig_hexa_nnode
833 
834  k = ig_hexa_node2gll(1,inode)
835  l = ig_hexa_node2gll(2,inode)
836  m = ig_hexa_node2gll(3,inode)
837 
838  global_gnode = ig_hexa_gnode_glonum(inode,ihexa)
839 ! myvar = sqrt(rg_hexa_gll_rhovs2(m,l,k,ihexa)/rg_hexa_gll_rho(m,l,k,ihexa))
840  myvar = var1(m,l,k,ihexa)
841 
842  var(global_gnode) = myvar
843 
844  enddo
845 
846  enddo
847 
848 !
849 !---->write array var in the xml VTK file
850  ios = vtk_var_xml( &
851  nc_nn = ig_mesh_nnode &
852  ,varname = 'var1' &
853  ,var = var &
854  )
855 
856 !
857 !---->fill the array var with gll_node value
858 
859  do ihexa = 1,ig_nhexa
860 
861  do inode = 1,ig_hexa_nnode
862 
863  k = ig_hexa_node2gll(1,inode)
864  l = ig_hexa_node2gll(2,inode)
865  m = ig_hexa_node2gll(3,inode)
866 
867  global_gnode = ig_hexa_gnode_glonum(inode,ihexa)
868 ! myvar = sqrt(rg_hexa_gll_rhovs2(m,l,k,ihexa)/rg_hexa_gll_rho(m,l,k,ihexa))
869  myvar = var2(m,l,k,ihexa)
870 
871  var(global_gnode) = myvar
872 
873  enddo
874 
875  enddo
876 
877 !
878 !---->write array var in the xml VTK file
879  ios = vtk_var_xml( &
880  nc_nn = ig_mesh_nnode &
881  ,varname = 'var2' &
882  ,var = var &
883  )
884 
885 !
886 !---->finalize xml VTK file
887  ios = vtk_dat_xml( &
888  var_location = 'NODE' &
889  ,var_block_action = 'CLOSE' &
890  )
891  ios = vtk_geo_xml()
892 
893  ios = vtk_end_xml()
894 
895  if (ig_myrank == 0) call write_medium_collection()
896 
897  return
898 !***********************************************************************************************************************************************************************************
899  end subroutine write_medium_vtk_node_xml
900 !***********************************************************************************************************************************************************************************
901 
902 !
903 !
907 !***********************************************************************************************************************************************************************************
908  subroutine write_medium_collection()
909 !***********************************************************************************************************************************************************************************
910 
911  use mod_global_variables, only :&
912  get_newunit&
913  ,cg_prefix&
914  ,ig_ncpu
915 
916 
917  implicit none
918 
919  integer :: icpu
920  integer :: myunit
921 
922  character(len=255) :: fname
923  character(len=6 ) :: cpart
924 
925  open(unit=get_newunit(myunit),file=trim(cg_prefix)//".medium.collection.pvd")
926 
927  write(unit=myunit,fmt='(a)') "<?xml version=""1.0""?>"
928  write(unit=myunit,fmt='(a)') "<VTKFile type=""Collection"" version=""0.1"" byte_order=""BigEndian"">"
929  write(unit=myunit,fmt='(a)') " <Collection>"
930 
931  do icpu = 1,ig_ncpu
932 
933  write(cpart,'(I6.6)') icpu-1
934 
935  fname = trim(cg_prefix)//".medium.cpu."//trim(cpart)//".vtu"
936 
937  write(unit=myunit,fmt='(3a)') " <DataSet group="""" part=""0"" file=""",trim(fname),"""/>"
938 
939  enddo
940 
941  write(unit=myunit,fmt='(a)') " </Collection>"
942  write(unit=myunit,fmt='(a)') "</VTKFile>"
943 
944  close(myunit)
945 
946  return
947 !***********************************************************************************************************************************************************************************
948  end subroutine write_medium_collection
949 !***********************************************************************************************************************************************************************************
950 
951 end module mod_init_medium
mod_init_memory
This module contains subroutines to allocate arrays and to compute an approximation of the total RAM ...
Definition: module_init_memory.f90:126
mod_vtk_io::vtk_end_xml
integer(i4p) function, public vtk_end_xml()
The VTK_END_XML function finalizes opened files.
Definition: module_vtk_io.f90:960
mod_init_medium
This module contains subroutines to initialize medium for hexahedron and quadrangle elements...
Definition: module_init_medium.f90:126
mod_vtk_io::vtk_con_xml
integer(i4p) function, public vtk_con_xml(NC, connect, offset, cell_type)
The VTK_CON_XML function must be called when unstructured grid topology is used. It saves the connect...
Definition: module_vtk_io.f90:620
mod_init_medium::init_quadp_medium
subroutine, public init_quadp_medium()
This subroutine fills medium properties at GLL nodes of paraxial quadrangle elements (i...
Definition: module_init_medium.f90:433
mod_vtk_io::vtk_ini_xml
integer(i4p) function, public vtk_ini_xml(output_format, filename, mesh_topology, nx1, nx2, ny1, ny2, nz1, nz2)
The VTK_INI_XML function is used for initializing file. This function must be the first to be called...
Definition: module_vtk_io.f90:288
mod_init_medium::write_medium_vtk_node_xml
subroutine, private write_medium_vtk_node_xml(var1, var2)
This subroutine writes S and P-wave velocities of the physical domain of cpu myrank in UnstructuredGr...
Definition: module_init_medium.f90:720
mod_vtk_io::vtk_geo_xml
overloading of VTK_GEO_XML
Definition: module_vtk_io.f90:157
mod_global_variables
This module defines all global variables of EFISPEC3D. Scalar variables are initialized directly in t...
Definition: module_global_variables.f90:128
mod_init_memory::init_array_real
Interface init_array_real to redirect allocation to n-rank arrays.
Definition: module_init_memory.f90:148
mod_init_medium::write_medium_vtk_cell_xml
subroutine, private write_medium_vtk_cell_xml()
This subroutine writes S and P-wave velocities of the physical domain of cpu myrank in UnstructuredGr...
Definition: module_init_medium.f90:585
mod_init_medium::init_medium_from_mat_file
subroutine, private init_medium_from_mat_file(il_mat_num_of_hexa)
This subroutine fills medium properties at GLL nodes of hexahedron elements based on material number ...
Definition: module_init_medium.f90:552
mod_vtk_io::vtk_dat_xml
integer(i4p) function, public vtk_dat_xml(var_location, var_block_action)
The VTK_DAT_XML function opens or closes CellData/PointData structures.
Definition: module_vtk_io.f90:695
mod_init_medium::write_medium_collection
subroutine, private write_medium_collection()
This subroutine writes a VTK xml collection file *.pvd which contains all *.vtu files generated by wr...
Definition: module_init_medium.f90:908
mod_vtk_io
This module programmed by S. Zaghi (https://github.com/szaghi/Lib_VTK_IO) contains functions to write...
Definition: module_vtk_io.f90:130
mod_vtk_io::vtk_var_xml
overloading of VTK_VAR_XML
Definition: module_vtk_io.f90:166
mod_init_medium::init_hexa_medium
subroutine, public init_hexa_medium()
This subroutine fills medium properties (elastic or viscoelastic) at GLL nodes of hexahedron elements...
Definition: module_init_medium.f90:158