All Classes Files Functions Variables Pages
module_init_efi.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_efi
127 
128  implicit none
129 
130  private
131 
132  public :: init_mpi
133  public :: init_input_variables
134  public :: init_temporal_domain
135  public :: init_temporal_saving
136  public :: init_gll_nodes
137  public :: init_gll_nodes_coordinates
138  public :: init_mass_matrix
139  public :: init_jacobian_matrix_hexa
140  public :: init_jacobian_matrix_quad
141 
142  contains
143 
150 !***********************************************************************************************************************************************************************************
151  subroutine init_mpi()
152 !***********************************************************************************************************************************************************************************
153 
154  use mpi
155 
156  use mod_global_variables, only :&
157  error_stop&
158  ,ig_ncpu&
159  ,cg_ncpu&
160  ,ig_myrank&
161  ,cg_myrank&
162  ,cg_cpu_name&
163  ,ig_cpu_name_len
164 
165  implicit none
166 
167  integer :: ios
168  character(len=255) :: info
169 
170  call mpi_init(ios)
171  if (ios /= 0) then
172  write(info,'(a)') "error while initializing mpi"
173  call error_stop(info)
174  endif
175 
176  call mpi_comm_size(mpi_comm_world,ig_ncpu,ios)
177  write(cg_ncpu,'(i6.6)') ig_ncpu
178 
179  call mpi_comm_rank(mpi_comm_world,ig_myrank,ios)
180  write(cg_myrank,'(i6.6)') ig_myrank
181 
182  call mpi_get_processor_name(cg_cpu_name,ig_cpu_name_len,ios)
183 
184  return
185 
186 !***********************************************************************************************************************************************************************************
187  end subroutine init_mpi
188 !***********************************************************************************************************************************************************************************
189 
190 !
191 !
194 !***********************************************************************************************************************************************************************************
195  subroutine init_input_variables()
196 !***********************************************************************************************************************************************************************************
197 
198  use mpi
199 
200  use mod_global_variables
201 
202  use mod_write_listing
203 
204  implicit none
205 
206  integer :: unit_input_file
207  integer :: ios
208  integer :: pos
209  integer :: imat
210  logical :: is_material_defined = .false.
211  logical :: is_file_exist = .false.
212  logical :: is_activate_dcs = .false.
213  logical :: is_activate_sfs = .false.
214 
215  character(len=100) :: buffer_line
216  character(len=100) :: buffer_value
217  character(len=100) :: keyword
218  character(len=255) :: info
219 
220 !
221 !---->Check if order of polynomial is correct
222  if (ig_lagrange_order < 4 .or. ig_lagrange_order > 6) then
223  write(info,'(a)') "error while checking order of lagrange polynomial"
224  call error_stop(info)
225  endif
226 
227 !
228 !---->read prefix
229  open(unit=100,file="prefix",iostat=ios)
230  if(ios /= 0) then
231  write(info,'(a)') "error while opening file prefix"
232  call error_stop(info)
233  endif
234  read(unit=100,fmt=*),cg_prefix
235  cg_prefix = trim(adjustl(cg_prefix))
236  close(100)
237 
238 !
239 !---->write header of simulation in listing file *.lst
240  call write_header()
241 
242 !
243 !
244 !*************************************************************************************************************
245 !---->process configuration file *.cfg
246 !*************************************************************************************************************
247  open(unit=get_newunit(unit_input_file),file=trim(cg_prefix)//".cfg",status='old',iostat=ios)
248  if (ios /= 0) then
249  write(info,'(a)') "error in subroutine init_input_variables while opening file *.cfg"
250  call error_stop(info)
251  endif
252 
253 !
254 !---->search for keyword in file *.cfg
255  do
256 
257  read(unit=unit_input_file,fmt='(a)',iostat=ios) buffer_line
258  if (ios /= 0) exit
259 
260  buffer_line = trim(adjustl(buffer_line))
261  pos = scan(buffer_line,'=')
262  if (pos >= 2) then
263  keyword = strlowcase(buffer_line(1:pos-1))
264  buffer_value = buffer_line(pos+1:100)
265  else
266  keyword = "na"
267  endif
268 
269  selectcase(trim(sweep_blanks(keyword)))
270 
271  case("durationofsimulation")
272  read(buffer_value,*) rg_simu_total_time
273 
274  case("timestep")
275  read(buffer_value,*) rg_dt
276  rg_dt2 = rg_dt**2
277 
278  case("receiversavingincrement")
279  read(buffer_value,*) ig_receiver_saving_incr
280 
281  case("snapshotsavingincrement", "surfacesnapshotsavingincrement")
282  read(buffer_value,*) ig_snapshot_saving_incr
283 
284  case("snapshotspaceincrement", "surfacesnapshotspaceincrement")
285  read(buffer_value,*) rg_receiver_snapshot_dxdy
286 
287  case("snapshotdisplacement", "surfacesnapshotdisplacement")
288  read(buffer_value,*) lg_snapshot_displacement
289 
290  case("snapshotvelocity", "surfacesnapshotvelocity")
291  read(buffer_value,*) lg_snapshot_velocity
292 
293  case("snapshotacceleration", "surfacesnapshotacceleration")
294  read(buffer_value,*) lg_snapshot_acceleration
295 
296  case("volumesnapshotsavingincrement")
297  read(buffer_value,*) ig_snapshot_volume_saving_incr
298 
299  case("volumesnapshotdisplacement")
300  read(buffer_value,*) lg_snapshot_volume_displacement
301 
302  case("volumesnapshotvelocity")
303  read(buffer_value,*) lg_snapshot_volume_velocity
304 
305  case("volumesnapshotacceleration")
306  read(buffer_value,*) lg_snapshot_volume_acceleration
307 
308  case("volumesnapshotxmin")
309  read(buffer_value,*) rg_snapshot_volume_xmin
310 
311  case("volumesnapshotxmax")
312  read(buffer_value,*) rg_snapshot_volume_xmax
313 
314  case("volumesnapshotymin")
315  read(buffer_value,*) rg_snapshot_volume_ymin
316 
317  case("volumesnapshotymax")
318  read(buffer_value,*) rg_snapshot_volume_ymax
319 
320  case("volumesnapshotzmin")
321  read(buffer_value,*) rg_snapshot_volume_zmin
322 
323  case("volumesnapshotzmax")
324  read(buffer_value,*) rg_snapshot_volume_zmax
325 
326  case("boundaryabsorption")
327  read(buffer_value,*) lg_boundary_absorption
328 
329  case("numberofmaterial")
330  read(buffer_value,*) ig_nmaterial
331  if (.not.allocated(tg_elastic_material)) allocate(tg_elastic_material(ig_nmaterial))
332  if (.not.allocated(ig_material_type)) allocate(ig_material_type(ig_nmaterial))
333  if (.not.allocated(tg_viscoelastic_material) .and. lg_visco) allocate(tg_viscoelastic_material(ig_nmaterial))
334 
335  case("material")
336  read(buffer_value,*) imat
337  is_material_defined = .true.
338  if (.not.lg_visco) then
339  ig_material_type(imat) = 1
340  else
341  ig_material_type(imat) = 2
342  endif
343 
344  case("vs")
345  read(buffer_value,*) tg_elastic_material(imat)%svel
346 
347  case("vp")
348  read(buffer_value,*) tg_elastic_material(imat)%pvel
349 
350  case("rho")
351  read(buffer_value,*) tg_elastic_material(imat)%dens
352 
353  case("qf")
354  if (lg_visco) then
355  read(buffer_value,*) tg_viscoelastic_material(imat)%freq
356  endif
357 
358  case("qs")
359  if (lg_visco) then
360  read(buffer_value,*) tg_viscoelastic_material(imat)%qs
361  endif
362 
363  case("qp")
364  if (lg_visco) then
365  read(buffer_value,*) tg_viscoelastic_material(imat)%qp
366  endif
367 
368  case default
369 
370  endselect
371 
372  enddo
373 
374  close(unit_input_file)
375 
376 
377 !
378 !
379 !*************************************************************************************************************
380 !check for illegal value of input variables
381 !*************************************************************************************************************
382  if (rg_simu_total_time <= tiny_real) then
383  write(info,'(a)') "error in subroutine init_input_variables: duration of simulation must be greater than TINY_REAL"
384  call error_stop(info)
385  endif
386 
387  do imat = 1,ig_nmaterial
388 
389  if (tg_elastic_material(imat)%svel <= tiny_real) then
390  write(info,'(a)') "error in subroutine init_input_variables: S-wave velocity must be greater than TINY_REAL"
391  call error_stop(info)
392  endif
393 
394  if (tg_elastic_material(imat)%pvel <= tiny_real) then
395  write(info,'(a)') "error in subroutine init_input_variables: P-wave velocity must be greater than TINY_REAL"
396  call error_stop(info)
397  endif
398 
399  if (tg_elastic_material(imat)%dens <= tiny_real) then
400  write(info,'(a)') "error in subroutine init_input_variables: density must be greater than TINY_REAL"
401  call error_stop(info)
402  endif
403 
404  if (lg_visco) then
405  if (tg_viscoelastic_material(imat)%freq <= tiny_real) then
406  write(info,'(a)') "error in subroutine init_input_variables: frequency for viscoelastic material must be greater than TINY_REAL"
407  call error_stop(info)
408  endif
409  if (tg_viscoelastic_material(imat)%qs <= tiny_real) then
410  write(info,'(a)') "error in subroutine init_input_variables: shear-wave quality factor for viscoelastic material must be greater than TINY_REAL"
411  call error_stop(info)
412  endif
413  if (tg_viscoelastic_material(imat)%qp <= tiny_real) then
414  write(info,'(a)') "error in subroutine init_input_variables: pressure-wave quality factor for viscoelastic material must be greater than TINY_REAL"
415  call error_stop(info)
416  endif
417  endif
418 
419  enddo
420 
421  if (.not.(is_material_defined)) then
422  write(info,'(a)') "error in subroutine init_input_variables: no material defined"
423  call error_stop(info)
424  endif
425 
426 
427 !
428 !
429 !*************************************************************************************************************
430 !set default value of input variables
431 !*************************************************************************************************************
432  if (rg_dt > tiny_real) then
433 
434  call init_temporal_domain(rg_simu_total_time,rg_dt)
435 
436  call init_temporal_saving(ig_ndt)
437 
438  endif
439 
440  call write_temporal_domain_info()
441 
442  inquire(file=trim(cg_prefix)//".dcs", exist=is_file_exist)
443  if(is_file_exist) then
444  is_activate_dcs = .true.
445  endif
446 
447  inquire(file=trim(cg_prefix)//".sfs", exist=is_file_exist)
448  if(is_file_exist) then
449  is_activate_sfs = .true.
450  endif
451 
452  if (.not.(is_activate_dcs) .and. .not.(is_activate_sfs)) then
453  write(info,'(a)') "error in subroutine init_input_variables: no source defined"
454  call error_stop(info)
455  endif
456 
457 
458  return
459 !***********************************************************************************************************************************************************************************
460  end subroutine init_input_variables
461 !***********************************************************************************************************************************************************************************
462 
463 
464 !
465 !
470 !***********************************************************************************************************************************************************************************
471  subroutine init_temporal_domain(t,dt)
472 !***********************************************************************************************************************************************************************************
473 
474  use mod_global_variables, only :&
475  ig_ndt&
476  ,rg_dt2
477 
478  implicit none
479 
480  real, intent(in) :: t
481  real, intent(in) :: dt
482 
483 
484  ig_ndt = int(t/dt) + 1
485 
486  rg_dt2 = dt**2
487 
488  return
489 !***********************************************************************************************************************************************************************************
490  end subroutine init_temporal_domain
491 !***********************************************************************************************************************************************************************************
492 
493 
494 !
495 !
503 !***********************************************************************************************************************************************************************************
504  subroutine init_temporal_saving(ndt)
505 !***********************************************************************************************************************************************************************************
506 
507  use mod_global_variables, only :&
508  ig_receiver_saving_incr&
509  ,ig_snapshot_saving_incr&
510  ,lg_snapshot&
511  ,lg_snapshot_displacement&
512  ,lg_snapshot_velocity&
513  ,lg_snapshot_acceleration&
514  ,ig_snapshot_volume_saving_incr&
515  ,lg_snapshot_volume&
516  ,lg_snapshot_volume_displacement&
517  ,lg_snapshot_volume_velocity&
518  ,lg_snapshot_volume_acceleration
519 
520  implicit none
521 
522  integer, intent(in) :: ndt
523 
524  if (ig_receiver_saving_incr == 0) then
525 
526  ig_receiver_saving_incr = ndt + 10
527 
528  endif
529 
530  if ( (ig_snapshot_saving_incr == 0) .or. (.not.lg_snapshot_displacement .and. .not.lg_snapshot_velocity .and. .not.lg_snapshot_acceleration) ) then
531 
532  ig_snapshot_saving_incr = ndt + 10
533  lg_snapshot = .false.
534  lg_snapshot_displacement = .false.
535  lg_snapshot_velocity = .false.
536  lg_snapshot_acceleration = .false.
537 
538  else
539 
540  lg_snapshot = .true.
541 
542  endif
543 
544  if ( (ig_snapshot_volume_saving_incr == 0) .or. (.not.lg_snapshot_volume_displacement .and. .not.lg_snapshot_volume_velocity .and. .not.lg_snapshot_volume_acceleration) ) then
545 
546  ig_snapshot_volume_saving_incr = ndt + 10
547  lg_snapshot_volume = .false.
548  lg_snapshot_volume_displacement = .false.
549  lg_snapshot_volume_velocity = .false.
550  lg_snapshot_volume_acceleration = .false.
551 
552  else
553 
554  lg_snapshot_volume = .true.
555 
556  endif
557 
558  return
559 !***********************************************************************************************************************************************************************************
560  end subroutine init_temporal_saving
561 !***********************************************************************************************************************************************************************************
562 
563 
564 !
565 !
577 !***********************************************************************************************************************************************************************************
578  subroutine init_gll_nodes()
579 !***********************************************************************************************************************************************************************************
580 
581  use mpi
582 
583  use mod_global_variables, only : &
584  cg_prefix&
585  ,rg_gll_abscissa_dist&
586  ,rg_gll_abscissa&
587  ,rg_gll_weight&
588  ,rg_gll_lagrange_deriv&
589  ,ig_ngll&
590  ,ig_lst_unit&
591  ,rg_pi&
592  ,ig_myrank&
593  ,ig_lagrange_order
594 
595  use mod_lagrange , only : lagrad
596 
597  implicit none
598 
599  real, parameter :: eps = 1.0e-7
600  real, allocatable ,dimension(:) :: xold
601  real, dimension(IG_NGLL,IG_NGLL) :: vandermonde_matrix
602  integer :: ngll
603  integer :: i
604  integer :: j
605  integer :: k
606  integer :: n
607 
608  !
609  !
610  !higher integration order among different group is taken as the integration order
611  n = ig_lagrange_order
612 
613 
614  !
615  !
616  !ngll is the number of integration node in 1d along the segment [-1;1]
617  ngll = n + 1
618 
619 
620  !
621  !
622  !*******************************
623  !find gll abscissa and weight
624  !*******************************
625  !
626  !use the chebyshev-gauss-lobatto nodes as the first guess
627  do i = 1,ngll
628  rg_gll_abscissa(i) = cos(rg_pi*real(i-1)/real(n))
629  enddo
630  !
631  !initialize the legendre vandermonde matrix
632  vandermonde_matrix(:,:) = 0.0
633  !
634  !compute p_{n} using the recursion relation. compute its first and second derivatives and update rg_gll_abscissa using the newton-raphson method
635  allocate(xold(ngll))
636  xold(:) = 2.0
637  do while (maxval(abs(rg_gll_abscissa-xold)) > eps)
638 
639  xold(:) = rg_gll_abscissa(:)
640  vandermonde_matrix(:,1) = 1.0
641  vandermonde_matrix(:,2) = rg_gll_abscissa(:)
642 
643  do k=2,n
644  vandermonde_matrix(:,k+1) = ( (2.0*real(k)-1.0)*rg_gll_abscissa(:)*vandermonde_matrix(:,k)-(real(k)-1.0)*vandermonde_matrix(:,k-1) )/real(k)
645  enddo
646 
647  rg_gll_abscissa(:) = xold(:)-( rg_gll_abscissa(:)*vandermonde_matrix(:,ngll)-vandermonde_matrix(:,n) )/( ngll*vandermonde_matrix(:,ngll) )
648 
649  enddo
650  !
651  !compute corresponding weight
652  rg_gll_weight(:) = 2.0/(n*ngll*vandermonde_matrix(:,ngll)**2)
653 
654  !
655  !
656  !*******************************************************************************
657  !write results rg_gll_abscissa,rg_gll_weight in file *.lst
658  !*******************************************************************************
659  if (ig_myrank == 0) then
660  write(ig_lst_unit,'(" ",/,a)') " --> abscissa of GLLs - weight of GLLs"
661  do i = 1,ngll
662  write(ig_lst_unit,fmt='(2(7X,f10.7))') rg_gll_abscissa(i),rg_gll_weight(i)
663  enddo
664  endif
665 
666  !
667  !
668  !*******************************************************************************
669  !compute subtraction between gll nodes (used in subroutines lagrap and lagrad)
670  !*******************************************************************************
671  do i = 1,ig_ngll
672  do j = 1,ig_ngll
673  rg_gll_abscissa_dist(j,i) = 0.0
674  if (i /= j) rg_gll_abscissa_dist(j,i) =1.0/(rg_gll_abscissa(i) - rg_gll_abscissa(j))
675  enddo
676  enddo
677 
678  !
679  !
680  !***************************************************
681  !compute derivative of lagrange poly at GLL nodes
682  !***************************************************
683  !
684  !notation convention: rg_gll_lagrange_deriv(j,i) = h'_j(xgll_i)
685  do i = 1,ig_ngll
686  do j = 1,ig_ngll
687  rg_gll_lagrange_deriv(j,i) = lagrad(j,rg_gll_abscissa(i),ig_ngll)
688  enddo
689  enddo
690 
691  return
692 
693 !***********************************************************************************************************************************************************************************
694  end subroutine init_gll_nodes
695 !***********************************************************************************************************************************************************************************
696 
697 !
698 !
702 !***********************************************************************************************************************************************************************************
703  subroutine init_gll_nodes_coordinates()
704 !***********************************************************************************************************************************************************************************
705 
706  use mpi
707 
708  use mod_global_variables, only :&
709  rg_gll_coordinate&
710  ,ig_ngll_total&
711  ,ig_nhexa&
712  ,ig_ngll&
713  ,ig_ndof&
714  ,ig_hexa_gll_glonum&
715  ,ig_hexa_nnode&
716  ,ig_hexa_gnode_xiloc&
717  ,ig_hexa_gnode_etloc&
718  ,ig_hexa_gnode_zeloc&
719  ,rg_gll_abscissa&
720  ,ig_line_nnode&
721  ,ig_hexa_gnode_glonum&
722  ,rg_gnode_x&
723  ,rg_gnode_y&
724  ,rg_gnode_z&
725  ,cg_prefix&
726  ,cg_myrank&
727  ,lg_output_debug_file&
728  ,error_stop
729 
730  use mod_init_memory
731 
732  use mod_lagrange, only : lagrap_geom
733 
734  implicit none
735 
736  real :: lag_xi
737  real :: lag_eta
738  real :: lag_zeta
739 
740  real :: xgnode
741  real :: ygnode
742  real :: zgnode
743 
744  integer :: ihexa
745  integer :: k
746  integer :: l
747  integer :: m
748  integer :: n
749  integer :: igll
750  integer :: ios
751 
752  integer(kind=1), dimension(ig_ngll_total) :: maskeq
753 
754  !
755  !
756  !*************************************************************************************
757  !compute global x,y,z-coordinates of gll nodes once and for all
758  !*************************************************************************************
759 
760  ios = init_array_real(rg_gll_coordinate,ig_ngll_total,ig_ndof,"rg_gll_coordinate")
761 
762  maskeq(:) = 0 !-->coordinate not computed yet
763 
764  do ihexa = 1,ig_nhexa
765 
766  do k = 1,ig_ngll !zeta
767  do l = 1,ig_ngll !eta
768  do m = 1,ig_ngll !xi
769 
770  igll = ig_hexa_gll_glonum(m,l,k,ihexa)
771 
772  if (maskeq(igll) == 0) then
773 
774  do n = 1,ig_hexa_nnode
775 
776  lag_xi = lagrap_geom(ig_hexa_gnode_xiloc(n),rg_gll_abscissa(m),ig_line_nnode)
777  lag_eta = lagrap_geom(ig_hexa_gnode_etloc(n),rg_gll_abscissa(l),ig_line_nnode)
778  lag_zeta = lagrap_geom(ig_hexa_gnode_zeloc(n),rg_gll_abscissa(k),ig_line_nnode)
779 
780  xgnode = rg_gnode_x(ig_hexa_gnode_glonum(n,ihexa))
781  ygnode = rg_gnode_y(ig_hexa_gnode_glonum(n,ihexa))
782  zgnode = rg_gnode_z(ig_hexa_gnode_glonum(n,ihexa))
783 
784  rg_gll_coordinate(1,igll) = rg_gll_coordinate(1,igll) + lag_xi*lag_eta*lag_zeta*xgnode
785 
786  rg_gll_coordinate(2,igll) = rg_gll_coordinate(2,igll) + lag_xi*lag_eta*lag_zeta*ygnode
787 
788  rg_gll_coordinate(3,igll) = rg_gll_coordinate(3,igll) + lag_xi*lag_eta*lag_zeta*zgnode
789 
790  enddo
791 
792  endif
793 
794  maskeq(igll) = 1
795 
796  enddo
797  enddo
798  enddo
799 
800  enddo
801 
802  if (lg_output_debug_file) then
803  open(unit=100,file="debug."//trim(cg_prefix)//".global.element.gll.coordinates.cpu."//trim(cg_myrank))
804 
805  do ihexa = 1,ig_nhexa
806  do k = 1,ig_ngll
807  do l = 1,ig_ngll
808  do m = 1,ig_ngll
809 
810  igll = ig_hexa_gll_glonum(m,l,k,ihexa)
811 
812  write(100,'(3(f15.3,1x),i10)') rg_gll_coordinate(1,igll)&
813  ,rg_gll_coordinate(2,igll)&
814  ,rg_gll_coordinate(3,igll)&
815  ,ig_hexa_gll_glonum(m,l,k,ihexa)
816  enddo
817  enddo
818  enddo
819  enddo
820 
821  close(100)
822  endif
823 
824  return
825 
826 !***********************************************************************************************************************************************************************************
827  end subroutine init_gll_nodes_coordinates
828 !***********************************************************************************************************************************************************************************
829 
830 !
831 !
834 !***********************************************************************************************************************************************************************************
835  subroutine init_mass_matrix()
836 !***********************************************************************************************************************************************************************************
837 
838  use mpi
839 
840  use mod_global_variables, only : &
841  ig_ngll&
842  ,ig_ngll_total&
843  ,rg_gll_mass_matrix&
844  ,rg_gll_weight&
845  ,ig_hexa_gll_glonum&
846  ,ig_nhexa&
847  ,ig_hexa_material_number&
848  ,rg_hexa_gll_jacobian_det&
849  ,ig_ncpu&
850  ,tg_cpu_neighbor&
851  ,ig_myrank&
852  ,rg_hexa_gll_rho&
853  ,error_stop&
854  ,ig_mpi_buffer_sizemax&
855  ,ig_ncpu_neighbor&
856  ,ig_lst_unit
857 
858  use mod_init_memory
859 
860  implicit none
861 
862  real , dimension(ig_mpi_buffer_sizemax) :: dummyr
863  real , dimension(ig_mpi_buffer_sizemax,ig_ncpu_neighbor) :: dlxmas
864 
865  integer, dimension(MPI_STATUS_SIZE) :: statut
866  integer :: ios
867  integer :: ngll
868  integer :: icon
869  integer :: ihexa
870  integer :: igll
871  integer :: kgll
872  integer :: lgll
873  integer :: mgll
874 
875  if (ig_myrank == 0) then
876  write(ig_lst_unit,'(" ",/,a)') "computing mass matrix..."
877  call flush(ig_lst_unit)
878  endif
879 
880  !
881  !
882  !*****************************************************************************************************************************************
883  !allocate mass matrix array
884  !*****************************************************************************************************************************************
885  ios = init_array_real(rg_gll_mass_matrix,ig_ngll_total,"rg_gll_mass_matrix")
886 
887  !
888  !
889  !*****************************************************************************************************************************************
890  !compute mass matrix for cpu 'myrank'
891  !*****************************************************************************************************************************************
892  do ihexa = 1,ig_nhexa
893  do kgll = 1,ig_ngll !zeta
894  do lgll = 1,ig_ngll !eta
895  do mgll = 1,ig_ngll !xi
896 
897  igll = ig_hexa_gll_glonum(mgll,lgll,kgll,ihexa)
898 
899  rg_gll_mass_matrix(igll) = rg_gll_mass_matrix(igll)&
900  + rg_hexa_gll_jacobian_det(mgll,lgll,kgll,ihexa)*rg_hexa_gll_rho(mgll,lgll,kgll,ihexa)&
901  * rg_gll_weight(mgll)*rg_gll_weight(lgll)*rg_gll_weight(kgll)
902  enddo
903  enddo
904  enddo
905  enddo
906  !
907  !
908  !*****************************************************************************************************************************************
909  !fill buffer dlxmas with gll to be send
910  !*****************************************************************************************************************************************
911  do icon = 1,ig_ncpu_neighbor
912  do igll = 1,tg_cpu_neighbor(icon)%ngll
913 
914  kgll = tg_cpu_neighbor(icon)%gll_send(igll)
915  dlxmas(igll,icon) = rg_gll_mass_matrix(kgll)
916 
917  enddo
918  enddo
919  !
920  !
921  !*****************************************************************************************************************************************
922  !assemble mass matrix on all cpus once and for all
923  !*****************************************************************************************************************************************
924  do icon = 1,ig_ncpu_neighbor
925 
926  ngll = tg_cpu_neighbor(icon)%ngll
927 
928  call mpi_sendrecv(dlxmas(1,icon),ngll,mpi_real,tg_cpu_neighbor(icon)%icpu,100,dummyr,ngll,mpi_real,tg_cpu_neighbor(icon)%icpu,100,mpi_comm_world,statut,ios)
929 
930  do igll = 1,ngll
931 
932  kgll = tg_cpu_neighbor(icon)%gll_recv(igll)
933  rg_gll_mass_matrix(kgll) = rg_gll_mass_matrix(kgll) + dummyr(igll)
934 
935  enddo
936 
937  enddo
938  !
939  !
940  !
941  !******************************************************************************************************************
942  !Compute the inverse of mass matrix
943  !******************************************************************************************************************
944  rg_gll_mass_matrix(:) = 1.0/rg_gll_mass_matrix(:)
945 
946  if (ig_myrank == 0) then
947  write(ig_lst_unit,'(a)') "done"
948  call flush(ig_lst_unit)
949  endif
950 
951 ! deallocate(rg_hexa_gll_rho)
952 
953  return
954 
955 !***********************************************************************************************************************************************************************************
956  end subroutine init_mass_matrix
957 !***********************************************************************************************************************************************************************************
958 
959 !
960 !
973 !***********************************************************************************************************************************************************************************
974  subroutine init_jacobian_matrix_hexa()
975 !***********************************************************************************************************************************************************************************
976 
977  use mpi
978 
979  use mod_global_variables, only : ig_ngll&
980  ,rg_gll_abscissa&
981  ,rg_gll_weight&
982  ,rg_gnode_x&
983  ,rg_gnode_y&
984  ,rg_gnode_z&
985  ,ig_hexa_gnode_xiloc&
986  ,ig_hexa_gnode_etloc&
987  ,ig_hexa_gnode_zeloc&
988  ,ig_nhexa&
989  ,ig_line_nnode&
990  ,rg_hexa_gll_jacobian_det&
991  ,rg_hexa_gll_dxidx&
992  ,rg_hexa_gll_dxidy&
993  ,rg_hexa_gll_dxidz&
994  ,rg_hexa_gll_detdx&
995  ,rg_hexa_gll_detdy&
996  ,rg_hexa_gll_detdz&
997  ,rg_hexa_gll_dzedx&
998  ,rg_hexa_gll_dzedy&
999  ,rg_hexa_gll_dzedz&
1000  ,ig_hexa_nnode&
1001  ,ig_myrank&
1002  ,error_stop&
1003  ,ig_lst_unit
1004 
1005  use mod_init_memory
1006 
1007  use mod_jacobian
1008 
1009  use mod_lagrange , only :&
1010  lagrap_geom&
1011  ,lagrad_geom
1012 
1013  implicit none
1014 
1015  real ,parameter :: eps = 1.0e-8
1016  real :: xi
1017  real :: eta
1018  real :: zeta
1019  real :: dxidx
1020  real :: dxidy
1021  real :: dxidz
1022  real :: detdx
1023  real :: detdy
1024  real :: detdz
1025  real :: dzedx
1026  real :: dzedy
1027  real :: dzedz
1028  real :: det
1029 
1030  integer :: ihexa
1031  integer :: k
1032  integer :: l
1033  integer :: m
1034  integer :: ios
1035 
1036  if (ig_myrank == 0) then
1037  write(ig_lst_unit,'(" ",/,a)') "computing jacobian matrix in hexa..."
1038  call flush(ig_lst_unit)
1039  endif
1040 
1041  !
1042  !
1043  !********************************************************************************************************************
1044  !initialize memory
1045  !********************************************************************************************************************
1046  ios = init_array_real(rg_hexa_gll_jacobian_det,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_jacobian_det")
1047 
1048  ios = init_array_real(rg_hexa_gll_dxidx,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dxidx")
1049 
1050  ios = init_array_real(rg_hexa_gll_dxidy,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dxidy")
1051 
1052  ios = init_array_real(rg_hexa_gll_dxidz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dxidz")
1053 
1054  ios = init_array_real(rg_hexa_gll_detdx,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_detdx")
1055 
1056  ios = init_array_real(rg_hexa_gll_detdy,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_detdy")
1057 
1058  ios = init_array_real(rg_hexa_gll_detdz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_detdz")
1059 
1060  ios = init_array_real(rg_hexa_gll_dzedx,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dzedx")
1061 
1062  ios = init_array_real(rg_hexa_gll_dzedy,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dzedy")
1063 
1064  ios = init_array_real(rg_hexa_gll_dzedz,ig_nhexa,ig_ngll,ig_ngll,ig_ngll,"rg_hexa_gll_dzedz")
1065 
1066 
1067  !
1068  !
1069  !********************************************************************************************************************
1070  !initialize jacobian matrix at each GLL nodes of each hexahedron elements
1071  !********************************************************************************************************************
1072  do ihexa = 1,ig_nhexa
1073 
1074  do k = 1,ig_ngll
1075  do l = 1,ig_ngll
1076  do m = 1,ig_ngll
1077 
1078  xi = rg_gll_abscissa(m)
1079  eta = rg_gll_abscissa(l)
1080  zeta = rg_gll_abscissa(k)
1081 
1082  call compute_hexa_jacobian(ihexa,xi,eta,zeta,dxidx,dxidy,dxidz,detdx,detdy,detdz,dzedx,dzedy,dzedz,det)
1083 
1084  rg_hexa_gll_jacobian_det(m,l,k,ihexa) = det
1085 
1086  rg_hexa_gll_dxidx(m,l,k,ihexa) = dxidx
1087  rg_hexa_gll_dxidy(m,l,k,ihexa) = dxidy
1088  rg_hexa_gll_dxidz(m,l,k,ihexa) = dxidz
1089 
1090  rg_hexa_gll_detdx(m,l,k,ihexa) = detdx
1091  rg_hexa_gll_detdy(m,l,k,ihexa) = detdy
1092  rg_hexa_gll_detdz(m,l,k,ihexa) = detdz
1093 
1094  rg_hexa_gll_dzedx(m,l,k,ihexa) = dzedx
1095  rg_hexa_gll_dzedy(m,l,k,ihexa) = dzedy
1096  rg_hexa_gll_dzedz(m,l,k,ihexa) = dzedz
1097 
1098  enddo
1099  enddo
1100  enddo
1101 
1102  enddo
1103 
1104  if (ig_myrank == 0) then
1105  write(ig_lst_unit,'(a)') "done"
1106  call flush(ig_lst_unit)
1107  endif
1108 
1109  return
1110 
1111 !***********************************************************************************************************************************************************************************
1112  end subroutine init_jacobian_matrix_hexa
1113 !***********************************************************************************************************************************************************************************
1114 
1115 !
1116 !
1121 !***********************************************************************************************************************************************************************************
1122  subroutine init_jacobian_matrix_quad()
1123 !***********************************************************************************************************************************************************************************
1124 
1125  use mpi
1126 
1127  use mod_global_variables, only : ig_ngll&
1128  ,rg_gll_abscissa&
1129  ,rg_gll_weight&
1130  ,rg_gnode_x&
1131  ,rg_gnode_y&
1132  ,rg_gnode_z&
1133  ,ig_line_nnode&
1134  ,ig_quad_gnode_xiloc&
1135  ,ig_quad_gnode_etloc&
1136  ,ig_nquad_parax&
1137  ,rg_quadp_gll_jaco_det&
1138  ,rg_quadp_gll_normal&
1139  ,ig_quad_nnode&
1140  ,ig_myrank&
1141  ,error_stop&
1142  ,ig_quadp_gnode_glonum&
1143  ,ig_lst_unit&
1144  ,lg_boundary_absorption
1145 
1146  use mod_init_memory
1147 
1148  use mod_lagrange , only :&
1149  lagrap_geom&
1150  ,lagrad_geom
1151 
1152  implicit none
1153 
1154  real ,parameter :: eps = 1.0e-8
1155  real :: xxi
1156  real :: xet
1157  real :: yxi
1158  real :: yet
1159  real :: zxi
1160  real :: zet
1161  real :: det
1162  real :: inv_det
1163  real :: crprx
1164  real :: crpry
1165  real :: crprz
1166  real :: vn(3)
1167  real :: lag_deriv_xi
1168  real :: lag_deriv_et
1169  real :: lag_xi
1170  real :: lag_et
1171  real :: xgnode
1172  real :: ygnode
1173  real :: zgnode
1174 
1175  integer :: iquad
1176  integer :: k
1177  integer :: l
1178  integer :: m
1179  integer :: ios
1180 
1181  !
1182  !
1183  !*********************************************************************************************************************
1184  !if there are quad (eg, paraxial elt): compute jacobian and derivative of local coordinate with respect to global ones
1185  !*********************************************************************************************************************
1186 
1187  if (lg_boundary_absorption .and. ig_nquad_parax > 0) then
1188 
1189  if (ig_myrank == 0) then
1190  write(ig_lst_unit,'(" ",/,a)') "computing jacobian matrix in quad..."
1191  call flush(ig_lst_unit)
1192  endif
1193 
1194  !
1195  !
1196  !********************************************************************************************************************
1197  !initialize memory
1198  !********************************************************************************************************************
1199  ios = init_array_real(rg_quadp_gll_jaco_det,ig_nquad_parax,ig_ngll,ig_ngll,"rg_quadp_gll_jaco_det")
1200 
1201  ios = init_array_real(rg_quadp_gll_normal,ig_nquad_parax,ig_ngll,ig_ngll,3,"rg_quadp_gll_normal")
1202 
1203 
1204  !
1205  !
1206  !********************************************************************************************************************
1207  !initialize jacobian matrix and normal vector at each GLL nodes of each quadrangle elements
1208  !********************************************************************************************************************
1209  do iquad = 1,ig_nquad_parax
1210  do k = 1,ig_ngll
1211  do l = 1,ig_ngll
1212 
1213  xxi = 0.0
1214  xet = 0.0
1215  yxi = 0.0
1216  yet = 0.0
1217  zxi = 0.0
1218  zet = 0.0
1219 
1220  do m = 1,ig_quad_nnode
1221 
1222  lag_deriv_xi = lagrad_geom(ig_quad_gnode_xiloc(m),rg_gll_abscissa(l),ig_line_nnode)
1223  lag_deriv_et = lagrad_geom(ig_quad_gnode_etloc(m),rg_gll_abscissa(k),ig_line_nnode)
1224 
1225  lag_xi = lagrap_geom(ig_quad_gnode_xiloc(m),rg_gll_abscissa(l),ig_line_nnode)
1226  lag_et = lagrap_geom(ig_quad_gnode_etloc(m),rg_gll_abscissa(k),ig_line_nnode)
1227 
1228  xgnode = rg_gnode_x(ig_quadp_gnode_glonum(m,iquad))
1229  ygnode = rg_gnode_y(ig_quadp_gnode_glonum(m,iquad))
1230  zgnode = rg_gnode_z(ig_quadp_gnode_glonum(m,iquad))
1231 
1232  xxi = xxi + lag_deriv_xi*lag_et*xgnode
1233  yxi = yxi + lag_deriv_xi*lag_et*ygnode
1234  zxi = zxi + lag_deriv_xi*lag_et*zgnode
1235 
1236  xet = xet + lag_xi*lag_deriv_et*xgnode
1237  yet = yet + lag_xi*lag_deriv_et*ygnode
1238  zet = zet + lag_xi*lag_deriv_et*zgnode
1239 
1240  enddo
1241 
1242  crprx = yxi*zet - zxi*yet
1243  crpry = zxi*xet - xxi*zet
1244  crprz = xxi*yet - yxi*xet
1245 
1246  det = sqrt(crprx**2 + crpry**2 + crprz**2) !the determinant of the quad element is the magnitude of the cross-product
1247  inv_det = 1.0/det
1248 
1249  if (det < eps) then
1250  write(*,*) "***error in jacobi quad: det null or negative, det = ",det
1251  write(*,*) "***element ",iquad, "cpu ",ig_myrank
1252  call mpi_abort(mpi_comm_world,100,ios)
1253  stop
1254  endif
1255  !!
1256  !!------------>determinant of the quad element
1257  rg_quadp_gll_jaco_det(l,k,iquad) = det
1258  !!
1259  !!------------>unit vector normal to the quad at the node k,l
1260  vn(1) = crprx*inv_det
1261  vn(2) = crpry*inv_det
1262  vn(3) = crprz*inv_det
1263 
1264  rg_quadp_gll_normal(1,l,k,iquad) = vn(1) !x
1265  rg_quadp_gll_normal(2,l,k,iquad) = vn(2) !y
1266  rg_quadp_gll_normal(3,l,k,iquad) = vn(3) !z
1267 
1268  enddo
1269  enddo
1270  enddo
1271 
1272  if (ig_myrank == 0) then
1273  write(ig_lst_unit,'(a)') "done"
1274  call flush(ig_lst_unit)
1275  endif
1276 
1277  endif
1278 
1279  return
1280 
1281 !***********************************************************************************************************************************************************************************
1282  end subroutine init_jacobian_matrix_quad
1283 !***********************************************************************************************************************************************************************************
1284 
1285 end module mod_init_efi
mod_init_efi::init_input_variables
subroutine, public init_input_variables()
This subroutine initializes the simulation by writing header of listing file *.lst and reading variab...
Definition: module_init_efi.f90:195
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_init_efi::init_gll_nodes
subroutine, public init_gll_nodes()
This subroutine computes GLL nodes abscissa and weight in the reference domain [-1:1] as well as deri...
Definition: module_init_efi.f90:578
mod_init_efi::init_mpi
subroutine, public init_mpi()
This subroutine initializes MPI.
Definition: module_init_efi.f90:151
mod_jacobian::compute_hexa_jacobian
subroutine, public compute_hexa_jacobian(ihexa, xisol, etsol, zesol, dxidx, dxidy, dxidz, detdx, detdy, detdz, dzedx, dzedy, dzedz, det)
This subroutine computes jacobian matrix and its determinant at location in hexahedron element ihexa...
Definition: module_jacobian.f90:158
mod_lagrange::lagrad_geom
real function, public lagrad_geom(i, x, n)
function to compute value of the derivative of order Lagrange polynomial of geometric node i at absc...
Definition: module_lagrange.f90:262
mod_write_listing::write_header
subroutine, public write_header()
subroutine that writes the header of the listing file (*.lst).
Definition: module_write_listing.f90:145
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_efi::init_jacobian_matrix_quad
subroutine, public init_jacobian_matrix_quad()
This subroutine computes the determinant of Jacobian matrix and normal unit vector of quadrangle elem...
Definition: module_init_efi.f90:1122
mod_init_memory::init_array_real
Interface init_array_real to redirect allocation to n-rank arrays.
Definition: module_init_memory.f90:148
mod_write_listing::write_temporal_domain_info
subroutine, public write_temporal_domain_info()
subroutine that writes temporal information about the simulation in the listing file (*...
Definition: module_write_listing.f90:245
mod_init_efi::init_temporal_saving
subroutine, public init_temporal_saving(ndt)
This subroutine set saving information for receivers and snapshots.
Definition: module_init_efi.f90:504
mod_init_efi
This module contains subroutines to initialize some global variable vectors and matrices.
Definition: module_init_efi.f90:126
mod_lagrange::lagrap_geom
real function, public lagrap_geom(i, x, n)
function to compute value of order Lagrange polynomial of geometric node i at abscissa : ...
Definition: module_lagrange.f90:181
mod_init_efi::init_gll_nodes_coordinates
subroutine, public init_gll_nodes_coordinates()
This subroutine computes GLL nodes -coordinates in the physical domain, cartesian right-handed coordi...
Definition: module_init_efi.f90:703
mod_write_listing
This module contains subroutines to write information in the listing file *.lst.
Definition: module_write_listing.f90:126
mod_init_efi::init_temporal_domain
subroutine, public init_temporal_domain(t, dt)
This subroutine initializes the number of time step of the simulation and the squared time step...
Definition: module_init_efi.f90:471
mod_lagrange
This module contains functions to compute Lagrange polynomials and its derivatives; and to interpolat...
Definition: module_lagrange.f90:126
mod_jacobian
This module contains subroutines to compute jacobian matrix.
Definition: module_jacobian.f90:126
mod_lagrange::lagrad
real function, public lagrad(i, x, n)
function to compute value of the derivative of order Lagrange polynomial of GLL node i at abscissa ...
Definition: module_lagrange.f90:213
mod_init_efi::init_mass_matrix
subroutine, public init_mass_matrix()
This subroutine computes and assembles the mass matrix of the system .
Definition: module_init_efi.f90:835
mod_init_efi::init_jacobian_matrix_hexa
subroutine, public init_jacobian_matrix_hexa()
This subroutine computes Jacobian matrix and its determinant for hexahedron elements.
Definition: module_init_efi.f90:974