All Classes Files Functions Variables Pages
module_init_time_step.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_time_step
127 
128  use mpi
129 
130  implicit none
131 
132  private
133 
134  public :: init_time_step
135  private :: round_time_step
136  private :: manexp
137 
138  contains
139 
140 
141 !
142 !
146 !***********************************************************************************************************************************************************************************
147  subroutine init_time_step()
148 !***********************************************************************************************************************************************************************************
149 
150  use mpi
151 
152  use mod_global_variables, only :&
153  ig_ngll&
154  ,tiny_real&
155  ,rg_simu_total_time&
156  ,rg_dt&
157  ,ig_ndt&
158  ,ig_nhexa&
159  ,rg_gll_coordinate&
160  ,rg_hexa_gll_rho&
161  ,rg_hexa_gll_rhovp2&
162  ,ig_hexa_gll_glonum&
163  ,lg_boundary_absorption&
164  ,error_stop
165 
166  use mod_init_efi, only :&
167  init_temporal_domain&
168  ,init_temporal_saving
169 
170  use mod_write_listing
171 
172  real :: dt_stable
173  real :: dt_stable_all_cpu
174 
175  real, dimension(ig_nhexa) :: hexa_dt_stable
176  real :: igll_dist_k
177  real :: igll_dist_l
178  real :: igll_dist_m
179  real :: igll_corner_x
180  real :: igll_corner_y
181  real :: igll_corner_z
182  real :: igll_k_x
183  real :: igll_k_y
184  real :: igll_k_z
185  real :: igll_l_x
186  real :: igll_l_y
187  real :: igll_l_z
188  real :: igll_m_x
189  real :: igll_m_y
190  real :: igll_m_z
191 
192  real :: igll_corner_vp
193  real :: igll_k_vp
194  real :: igll_l_vp
195  real :: igll_m_vp
196  real :: m_vp_max
197  real :: l_vp_max
198  real :: k_vp_max
199 
200  real :: c_max
201 
202  integer :: ihexa
203  integer :: kgll
204  integer :: lgll
205  integer :: mgll
206  integer :: igll_corner
207  integer :: igll_k
208  integer :: igll_l
209  integer :: igll_m
210  integer :: kshift
211  integer :: lshift
212  integer :: mshift
213  integer :: ios
214 
215 
216  if(lg_boundary_absorption) then
217 
218  c_max = 0.29
219 
220  else
221 
222  c_max = 0.59
223 
224  endif
225 
226 
227  do ihexa = 1,ig_nhexa
228 
229  hexa_dt_stable(ihexa) = huge(c_max)
230 
231  do kgll = 1,ig_ngll,ig_ngll-1
232 
233  kshift = kgll - 1 + 2*mod(kgll,ig_ngll)
234 
235  do lgll = 1,ig_ngll,ig_ngll-1
236 
237  lshift = lgll - 1 + 2*mod(lgll,ig_ngll)
238 
239  do mgll = 1,ig_ngll,ig_ngll-1
240 
241  mshift = mgll - 1 + 2*mod(mgll,ig_ngll)
242 
243  igll_corner = ig_hexa_gll_glonum(mgll ,lgll ,kgll ,ihexa)
244  igll_k = ig_hexa_gll_glonum(mgll ,lgll ,kshift,ihexa)
245  igll_l = ig_hexa_gll_glonum(mgll ,lshift,kgll ,ihexa)
246  igll_m = ig_hexa_gll_glonum(mshift,lgll ,kgll ,ihexa)
247 
248 !
249 !------------------->get corner GLL coordinates
250  igll_corner_x = rg_gll_coordinate(1,igll_corner)
251  igll_corner_y = rg_gll_coordinate(2,igll_corner)
252  igll_corner_z = rg_gll_coordinate(3,igll_corner)
253 
254 !
255 !------------------->compute the distance between the corner GLL nodes and the closest GLL node towards m-direction
256  igll_m_x = rg_gll_coordinate(1,igll_m)
257  igll_m_y = rg_gll_coordinate(2,igll_m)
258  igll_m_z = rg_gll_coordinate(3,igll_m)
259 
260  igll_dist_m = sqrt((igll_corner_x-igll_m_x)**2 + (igll_corner_y-igll_m_y)**2 + (igll_corner_z-igll_m_z)**2)
261 
262 !
263 !------------------->compute the distance between the corner GLL nodes and the closest GLL node towards l-direction
264  igll_l_x = rg_gll_coordinate(1,igll_l)
265  igll_l_y = rg_gll_coordinate(2,igll_l)
266  igll_l_z = rg_gll_coordinate(3,igll_l)
267 
268  igll_dist_l = sqrt((igll_corner_x-igll_l_x)**2 +(igll_corner_y-igll_l_y)**2 + (igll_corner_z-igll_l_z)**2)
269 
270 !
271 !------------------->compute the distance between the corner GLL nodes and the closest GLL node towards k-direction
272  igll_k_x = rg_gll_coordinate(1,igll_k)
273  igll_k_y = rg_gll_coordinate(2,igll_k)
274  igll_k_z = rg_gll_coordinate(3,igll_k)
275 
276  igll_dist_k = sqrt((igll_corner_x-igll_k_x)**2 +(igll_corner_y-igll_k_y)**2 + (igll_corner_z-igll_k_z)**2)
277 
278 !
279 !------------------->get P-wave velocity around corner kgll,lgll,mgll
280  igll_corner_vp = sqrt(rg_hexa_gll_rhovp2(mgll ,lgll ,kgll ,ihexa)/rg_hexa_gll_rho(mgll ,lgll ,kgll ,ihexa))
281  igll_k_vp = sqrt(rg_hexa_gll_rhovp2(mgll ,lgll ,kshift,ihexa)/rg_hexa_gll_rho(mgll ,lgll ,kshift,ihexa))
282  igll_l_vp = sqrt(rg_hexa_gll_rhovp2(mgll ,lshift,kgll ,ihexa)/rg_hexa_gll_rho(mgll ,lshift,kgll ,ihexa))
283  igll_m_vp = sqrt(rg_hexa_gll_rhovp2(mshift,lgll ,kgll ,ihexa)/rg_hexa_gll_rho(mshift,lgll ,kgll ,ihexa))
284 
285 !
286 !------------------->compute maximal P-wave velocity around corner kgll,lgll,mgll
287  m_vp_max = max(igll_corner_vp,igll_m_vp)
288  l_vp_max = max(igll_corner_vp,igll_l_vp)
289  k_vp_max = max(igll_corner_vp,igll_k_vp)
290 
291  hexa_dt_stable(ihexa) = min(hexa_dt_stable(ihexa),c_max * igll_dist_k / igll_k_vp , c_max * igll_dist_l / igll_l_vp , c_max * igll_dist_m / igll_m_vp)
292 
293  enddo
294  enddo
295  enddo
296 
297  enddo
298 
299  dt_stable = minval(hexa_dt_stable)
300 
301  call round_time_step(dt_stable)
302 
303  call mpi_allreduce(dt_stable,dt_stable_all_cpu,1,mpi_real,mpi_min,mpi_comm_world,ios)
304 
305  if ( (rg_dt <= tiny_real) .or. (rg_dt > dt_stable_all_cpu) ) then
306 
307  rg_dt = dt_stable_all_cpu
308 
309  call init_temporal_domain(rg_simu_total_time,rg_dt)
310 
311  call init_temporal_saving(ig_ndt)
312 
313  call write_cfl_condition_ko(rg_dt)
314 
315  else
316 
317  call write_cfl_condition_ok(rg_dt)
318 
319  endif
320 
321  return
322 !***********************************************************************************************************************************************************************************
323  end subroutine init_time_step
324 !***********************************************************************************************************************************************************************************
325 
326 
327 !
328 !
332 !***********************************************************************************************************************************************************************************
333  subroutine round_time_step(x)
334 !***********************************************************************************************************************************************************************************
335 
336  real, intent(inout) :: x
337 
338  integer :: e
339  real :: m
340 
341 !
342 !------->get the mantissa and exponent of a real number
343  call manexp(x,m,e)
344 
345 !
346 !------->round the 10*mantissa to .0 or .5
347  if ( (10.0*m - floor(10.0*m)) >= 0.5 ) then
348 
349  x = ((floor(10.0*m)+0.5)/10.0)*10.0**(e)
350 
351  else
352 
353  x = ((floor(10.0*m))/10.0)*10.0**(e)
354 
355  endif
356 
357  return
358 
359 !***********************************************************************************************************************************************************************************
360  end subroutine round_time_step
361 !***********************************************************************************************************************************************************************************
362 
363 
364 !
365 !
370 !***********************************************************************************************************************************************************************************
371  subroutine manexp(x,m,e)
372 !***********************************************************************************************************************************************************************************
373 
374  real , intent(in ) :: x
375  integer, intent(out) :: e
376  real , intent(out) :: m
377 
378  real :: xx
379 
380  if (x < 0.0) then
381 
382  xx = -x
383  e = int(log10(xx))
384  if (e > 0) e = e + 1
385  m = - xx * 1.0e1 ** (-e)
386 
387  elseif (x > 0.0) then
388 
389  e = int(log10(x))
390  if (e > 0) e = e + 1
391  m = x * 1.0e1 ** (-e)
392 
393  else
394 
395  e = 0.0
396  m = 0.0
397 
398  endif
399 
400  return
401 
402 !***********************************************************************************************************************************************************************************
403  end subroutine manexp
404 !***********************************************************************************************************************************************************************************
405 
406 end module mod_init_time_step
mod_init_time_step::init_time_step
subroutine, public init_time_step()
This subroutine returns the time step mod_global_variables::rg_dt that makes the simulation stable...
Definition: module_init_time_step.f90:147
mod_init_time_step::round_time_step
subroutine, private round_time_step(x)
This subroutine rounds down the mantissa of the scientific notation of a real number to ...
Definition: module_init_time_step.f90:333
mod_write_listing::write_cfl_condition_ok
subroutine, public write_cfl_condition_ok(dt)
subroutine that writes valid CFL condition in the listing file (*.lst).
Definition: module_write_listing.f90:406
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_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_write_listing::write_cfl_condition_ko
subroutine, public write_cfl_condition_ko(dt)
subroutine that writes invalid CFL condition in the listing file (*.lst).
Definition: module_write_listing.f90:433
mod_write_listing
This module contains subroutines to write information in the listing file *.lst.
Definition: module_write_listing.f90:126
mod_init_time_step::manexp
subroutine, private manexp(x, m, e)
This subroutine returns the mantissa and the exponent of the exponential notation of a real number...
Definition: module_init_time_step.f90:371
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_init_time_step
This module contains subroutines to initialize the time step of the simulation.
Definition: module_init_time_step.f90:126