All Classes Files Functions Variables Pages
module_source.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_source
127 
128  use mpi
129 
130  implicit none
131 
132  public :: compute_double_couple_source
133  public :: init_double_couple_source
134  public :: init_single_force_source
135  private :: compute_source_local_coordinate
136 
137  contains
138 
139 !
140 !
146 !***********************************************************************************************************************************************************************************
147  subroutine compute_double_couple_source()
148 !***********************************************************************************************************************************************************************************
149 
150  use mpi
151 
152  use mod_global_variables, only : ig_ndcsource&
153  ,ig_ngll&
154  ,ig_ndof&
155  ,tg_dcsource&
156  ,rg_dcsource_gll_force&
157  ,rg_hexa_gll_dxidx&
158  ,rg_hexa_gll_dxidy&
159  ,rg_hexa_gll_dxidz&
160  ,rg_hexa_gll_detdx&
161  ,rg_hexa_gll_detdy&
162  ,rg_hexa_gll_detdz&
163  ,rg_hexa_gll_dzedx&
164  ,rg_hexa_gll_dzedy&
165  ,rg_hexa_gll_dzedz&
166  ,ig_myrank&
167  ,error_stop&
168  ,ig_lst_unit
169 
170  use mod_init_memory
171 
172  use mod_lagrange, only : lagrad,lagrap
173 
174  implicit none
175 
176  real :: ftmp(3,3,ig_ngll,ig_ngll,ig_ngll)
177 
178  integer :: iso
179  integer :: iel
180  integer :: k
181  integer :: l
182  integer :: m
183  integer :: n
184  integer :: o
185  integer :: p
186  integer :: ios
187 
188 
189 
190  !**************************************************************************************************************!
191  !References: Komatitsch and Tromp, 1999 \cite Komatitsch1999a; Komatitsch and Tromp, 2002 \cite Komatitsch2002b!
192  !**************************************************************************************************************!
193 
194 
195  !
196  !
197  !***********************************************************
198  !implementation of the double couple point source
199  !***********************************************************
200  if (ig_myrank == 0) then
201  write(ig_lst_unit,'(" ",/,a)') "computing double-couple sources if any..."
202  call flush(ig_lst_unit)
203  endif
204 
205  if (ig_ndcsource > 0) then
206 
207  ios = init_array_real(rg_dcsource_gll_force,ig_ndcsource,ig_ngll,ig_ngll,ig_ngll,ig_ndof,"rg_dcsource_gll_force")
208 
209  endif
210 
211  !
212  !inject double couple points sources at any location (i.e., not necessarily at a gll node)
213  do iso = 1,ig_ndcsource
214 
215  iel = tg_dcsource(iso)%iel
216  !
217  !->compute first part at gll node first, then an interpolation (cf. below) is done at the xi, eta, zeta of the source
218  do k = 1,ig_ngll
219  do l = 1,ig_ngll
220  do m = 1,ig_ngll
221  ftmp(1,1,m,l,k) = tg_dcsource(iso)%mxx*rg_hexa_gll_dxidx(m,l,k,iel) &
222  + tg_dcsource(iso)%mxy*rg_hexa_gll_dxidy(m,l,k,iel) &
223  + tg_dcsource(iso)%mxz*rg_hexa_gll_dxidz(m,l,k,iel)
224  ftmp(2,1,m,l,k) = tg_dcsource(iso)%mxx*rg_hexa_gll_detdx(m,l,k,iel) &
225  + tg_dcsource(iso)%mxy*rg_hexa_gll_detdy(m,l,k,iel) &
226  + tg_dcsource(iso)%mxz*rg_hexa_gll_detdz(m,l,k,iel)
227  ftmp(3,1,m,l,k) = tg_dcsource(iso)%mxx*rg_hexa_gll_dzedx(m,l,k,iel) &
228  + tg_dcsource(iso)%mxy*rg_hexa_gll_dzedy(m,l,k,iel) &
229  + tg_dcsource(iso)%mxz*rg_hexa_gll_dzedz(m,l,k,iel)
230 
231  ftmp(1,2,m,l,k) = tg_dcsource(iso)%mxy*rg_hexa_gll_dxidx(m,l,k,iel) &
232  + tg_dcsource(iso)%myy*rg_hexa_gll_dxidy(m,l,k,iel) &
233  + tg_dcsource(iso)%myz*rg_hexa_gll_dxidz(m,l,k,iel)
234  ftmp(2,2,m,l,k) = tg_dcsource(iso)%mxy*rg_hexa_gll_detdx(m,l,k,iel) &
235  + tg_dcsource(iso)%myy*rg_hexa_gll_detdy(m,l,k,iel) &
236  + tg_dcsource(iso)%myz*rg_hexa_gll_detdz(m,l,k,iel)
237  ftmp(3,2,m,l,k) = tg_dcsource(iso)%mxy*rg_hexa_gll_dzedx(m,l,k,iel) &
238  + tg_dcsource(iso)%myy*rg_hexa_gll_dzedy(m,l,k,iel) &
239  + tg_dcsource(iso)%myz*rg_hexa_gll_dzedz(m,l,k,iel)
240 
241  ftmp(1,3,m,l,k) = tg_dcsource(iso)%mxz*rg_hexa_gll_dxidx(m,l,k,iel) &
242  + tg_dcsource(iso)%myz*rg_hexa_gll_dxidy(m,l,k,iel) &
243  + tg_dcsource(iso)%mzz*rg_hexa_gll_dxidz(m,l,k,iel)
244  ftmp(2,3,m,l,k) = tg_dcsource(iso)%mxz*rg_hexa_gll_detdx(m,l,k,iel) &
245  + tg_dcsource(iso)%myz*rg_hexa_gll_detdy(m,l,k,iel) &
246  + tg_dcsource(iso)%mzz*rg_hexa_gll_detdz(m,l,k,iel)
247  ftmp(3,3,m,l,k) = tg_dcsource(iso)%mxz*rg_hexa_gll_dzedx(m,l,k,iel) &
248  + tg_dcsource(iso)%myz*rg_hexa_gll_dzedy(m,l,k,iel) &
249  + tg_dcsource(iso)%mzz*rg_hexa_gll_dzedz(m,l,k,iel)
250  enddo
251  enddo
252  enddo
253  !
254  !->interpolation at the xi, eta, zeta of the source + multiplication by the test function
255  do k = 1,ig_ngll
256  do l = 1,ig_ngll
257  do m = 1,ig_ngll
258  rg_dcsource_gll_force(:,m,l,k,iso) = 0.0
259  do n = 1,ig_ngll
260  do o = 1,ig_ngll
261  do p = 1,ig_ngll
262  rg_dcsource_gll_force(1,m,l,k,iso) = rg_dcsource_gll_force(1,m,l,k,iso) &
263  + lagrap(p,tg_dcsource(iso)%xi,ig_ngll) &
264  * lagrap(o,tg_dcsource(iso)%et,ig_ngll) &
265  * lagrap(n,tg_dcsource(iso)%ze,ig_ngll) &
266  *(ftmp(1,1,p,o,n) * lagrad(m,tg_dcsource(iso)%xi,ig_ngll) &
267  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
268  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
269  + ftmp(2,1,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
270  * lagrad(l,tg_dcsource(iso)%et,ig_ngll) &
271  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
272  + ftmp(3,1,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
273  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
274  * lagrad(k,tg_dcsource(iso)%ze,ig_ngll))
275 
276  rg_dcsource_gll_force(2,m,l,k,iso) = rg_dcsource_gll_force(2,m,l,k,iso) &
277  + lagrap(p,tg_dcsource(iso)%xi,ig_ngll) &
278  * lagrap(o,tg_dcsource(iso)%et,ig_ngll) &
279  * lagrap(n,tg_dcsource(iso)%ze,ig_ngll) &
280  *(ftmp(1,2,p,o,n) * lagrad(m,tg_dcsource(iso)%xi,ig_ngll) &
281  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
282  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
283  + ftmp(2,2,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
284  * lagrad(l,tg_dcsource(iso)%et,ig_ngll) &
285  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
286  + ftmp(3,2,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
287  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
288  * lagrad(k,tg_dcsource(iso)%ze,ig_ngll))
289 
290  rg_dcsource_gll_force(3,m,l,k,iso) = rg_dcsource_gll_force(3,m,l,k,iso) &
291  + lagrap(p,tg_dcsource(iso)%xi,ig_ngll) &
292  * lagrap(o,tg_dcsource(iso)%et,ig_ngll) &
293  * lagrap(n,tg_dcsource(iso)%ze,ig_ngll) &
294  *(ftmp(1,3,p,o,n) * lagrad(m,tg_dcsource(iso)%xi,ig_ngll) &
295  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
296  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
297  + ftmp(2,3,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
298  * lagrad(l,tg_dcsource(iso)%et,ig_ngll) &
299  * lagrap(k,tg_dcsource(iso)%ze,ig_ngll) &
300  + ftmp(3,3,p,o,n) * lagrap(m,tg_dcsource(iso)%xi,ig_ngll) &
301  * lagrap(l,tg_dcsource(iso)%et,ig_ngll) &
302  * lagrad(k,tg_dcsource(iso)%ze,ig_ngll))
303  enddo
304  enddo
305  enddo
306  enddo
307  enddo
308  enddo
309  enddo
310 
311  if (ig_myrank == 0) then
312  write(ig_lst_unit,'(a)') "done"
313  call flush(ig_lst_unit)
314  endif
315 
316  return
317 !***********************************************************************************************************************************************************************************
318  end subroutine compute_double_couple_source
319 !***********************************************************************************************************************************************************************************
320 
321 !
322 !
329 !***********************************************************************************************************************************************************************************
330  subroutine init_double_couple_source()
331 !***********************************************************************************************************************************************************************************
332 
333  use mpi
334 
335  use mod_global_variables, only :&
336  ig_lst_unit&
337  ,tg_dcsource&
338  ,ig_ndcsource&
339  ,ig_ncpu&
340  ,ig_myrank&
341  ,type_double_couple_source&
342  ,cg_prefix&
343  ,rg_dt&
344  ,ig_ndt&
345  ,rg_pi&
346  ,rg_dcsource_user_func
347 
348  use mod_source_function , only : interp_linear
349 
350  use mod_receiver , only : search_closest_hexa_gll
351 
352  implicit none
353 
354  real, allocatable, dimension(:) :: dum1
355  real, allocatable, dimension(:) :: dum2
356  real, allocatable, dimension(:) :: dum3
357  real , dimension(ig_ncpu) :: dummy
358  real :: st
359  real :: di
360  real :: ra
361  real :: m0
362  real :: rldmin
363  real :: maxdmin
364 
365  integer , dimension(MPI_STATUS_SIZE) :: statut
366  integer :: min_loc(1)
367  integer , dimension(ig_ncpu) :: ilrcpu
368  integer :: i
369  integer :: j
370  integer :: iso
371  integer :: ios
372  integer :: ios1
373  integer :: ilndco
374  integer :: itempo
375 
376  character(len=1) :: ctempo
377 
378  type(type_double_couple_source), allocatable :: tldoco(:)
379 
380  !
381  !
382  !********************************************************************
383  !find element containing double couple point source in cpu 'myrank'
384  !********************************************************************
385 
386  ilndco = 0
387 
388  open(unit=100,file=trim(cg_prefix)//".dcs",status='old',iostat=ios)
389 
390  if (ios == 0) then
391 
392  if (ig_myrank == 0) write(ig_lst_unit,'(a)') " "
393 
394  read(unit=100,fmt='(i10)') ig_ndcsource
395 
396  allocate(tldoco(ig_ndcsource))
397 
398  do iso = 1,ig_ndcsource
399 
400  !initialize tldoco
401  tldoco(iso)%x = 0.0
402  tldoco(iso)%y = 0.0
403  tldoco(iso)%z = 0.0
404  tldoco(iso)%mxx = 0.0
405  tldoco(iso)%myy = 0.0
406  tldoco(iso)%mzz = 0.0
407  tldoco(iso)%mxy = 0.0
408  tldoco(iso)%mxz = 0.0
409  tldoco(iso)%myz = 0.0
410  tldoco(iso)%shift_time = 0.0
411  tldoco(iso)%rise_time = 0.0
412  tldoco(iso)%xi = 0.0
413  tldoco(iso)%et = 0.0
414  tldoco(iso)%ze = 0.0
415  tldoco(iso)%dmin = 0.0
416  tldoco(iso)%str = 0.0
417  tldoco(iso)%dip = 0.0
418  tldoco(iso)%rak = 0.0
419  tldoco(iso)%mw = 0.0
420  tldoco(iso)%icur = 0
421  tldoco(iso)%cpu = 0
422  tldoco(iso)%iel = 0
423  tldoco(iso)%kgll = 0
424  tldoco(iso)%lgll = 0
425  tldoco(iso)%mgll = 0
426 
427 
428  read(unit=100,fmt=*) tldoco(iso)%x,tldoco(iso)%y,tldoco(iso)%z,tldoco(iso)%mw,tldoco(iso)%str,tldoco(iso)%dip,tldoco(iso)%rak,tldoco(iso)%shift_time,tldoco(iso)%rise_time,tldoco(iso)%icur
429 
430  st = tldoco(iso)%str*rg_pi/180.0
431  di = tldoco(iso)%dip*rg_pi/180.0
432  ra = tldoco(iso)%rak*rg_pi/180.0
433  m0 = 10**(1.5*tldoco(iso)%mw + 9.1)
434  tldoco(iso)%mxx = +m0*(sin(di)*cos(ra)*sin(2.0*st) - sin(2.0*di)*sin(ra)*(cos(st))**2) !+myy for aki
435  tldoco(iso)%mxy = +m0*(sin(di)*cos(ra)*cos(2.0*st) + 0.5*sin(2.0*di)*sin(ra)* sin(2.0*st)) !+mxy
436  tldoco(iso)%mxz = +m0*(cos(di)*cos(ra)*sin( st) - cos(2.0*di)*sin(ra)* cos(st)) !-myz
437  tldoco(iso)%myy = -m0*(sin(di)*cos(ra)*sin(2.0*st) + sin(2.0*di)*sin(ra)*(sin(st))**2) !+mxx
438  tldoco(iso)%myz = +m0*(cos(di)*cos(ra)*cos( st) + cos(2.0*di)*sin(ra)* sin(st)) !-mxz
439  tldoco(iso)%mzz = +m0*sin(2.0*di)*sin(ra) !+mzz
440 
441  if ( (iso == 1) .and. (tldoco(iso)%icur == 10) ) then
442  itempo = 0
443  open(unit=99,file=trim(cg_prefix)//".dcf")
444  do while (.true.)
445  read(unit= 99,fmt=*,iostat=ios1) ctempo
446  if (ios1 /= 0) exit
447  itempo = itempo + 1
448  enddo
449  rewind(99)
450  allocate(rg_dcsource_user_func(ig_ndt+1),dum3(ig_ndt+1),dum1(itempo),dum2(itempo))
451  do i = 1,itempo
452  read(unit=99,fmt=*) dum1(i),dum2(i)
453  enddo
454  close(99)
455  do i = 1,ig_ndt+1
456  dum3(i) = (i-1)*rg_dt
457  enddo
458  call interp_linear(1,itempo,dum1,dum2,ig_ndt+1,dum3,rg_dcsource_user_func)
459  if (ig_myrank == 0) then
460  open(unit=90,file=trim(cg_prefix)//".dcf.interp")
461  do i = 1,ig_ndt+1
462  write(unit=90,fmt='(2(E14.7,1X))') dum3(i),rg_dcsource_user_func(i)
463  enddo
464  close(90)
465  endif
466  endif
467  !
468  !->find the closest element and gll point to the source iso
469  call search_closest_hexa_gll(tldoco(iso)%x,tldoco(iso)%y,tldoco(iso)%z &
470  ,tldoco(iso)%dmin,tldoco(iso)%kgll,tldoco(iso)%lgll,tldoco(iso)%mgll,tldoco(iso)%iel)
471  !
472  !->check which cpu has the smallest dmin. this cpu will compute this source and the other won't
473  call mpi_allgather(tldoco(iso)%dmin,1,mpi_real,dummy,1,mpi_real,mpi_comm_world,ios)
474  min_loc = minloc(dummy(1:ig_ncpu))
475  if ( (min_loc(1)-1) == ig_myrank) then
476  tldoco(iso)%cpu = ig_myrank
477  ilndco = ilndco + 1
478  else
479  tldoco(iso)%cpu = -1
480  endif
481  if (ig_myrank == 0) then
482  write(ig_lst_unit,'(a,i8,a,i0)') "double couple source ",iso," computed by cpu ",min_loc(1)-1
483  endif
484  enddo
485  close(100)
486  !
487  !->transfer array tldoco which knows all the sources to array tg_dcsource that needs only to know the sources of cpu 'ig_myrank'
488  if (ilndco > 0) then
489  allocate(tg_dcsource(ilndco))
490  j = 0
491  do iso = 1,ig_ndcsource
492  if (tldoco(iso)%cpu == ig_myrank) then
493  j = j + 1
494  tg_dcsource(j) = tldoco(iso)
495  endif
496  enddo
497  ig_ndcsource = ilndco
498  else
499  ig_ndcsource = 0
500  endif
501  deallocate(tldoco)
502  !
503  !->cpu 0 gather the number of source (ig_ndcsource) of the other cpus
504  call mpi_gather(ig_ndcsource,1,mpi_integer,ilrcpu,1,mpi_integer,0,mpi_comm_world,ios)
505  !
506  do iso = 1,ig_ndcsource
507  call compute_source_local_coordinate(tg_dcsource(iso))
508  if (ig_myrank /= 0) call mpi_send(tg_dcsource(iso)%dmin,1,mpi_real,0,100,mpi_comm_world,ios)
509  enddo
510  !
511  !->cpu 0 gathers dmin of all cpus and write the maximum dmin in *.lst
512  if (ig_myrank == 0) then
513  maxdmin = 0.0
514  do i = 1,ig_ncpu
515  if (i == 1) then
516  do iso = 1,ig_ndcsource
517  maxdmin = max(maxdmin,tg_dcsource(iso)%dmin)
518  enddo
519  else
520  do iso = 1,ilrcpu(i)
521  call mpi_recv(rldmin,1,mpi_real,i-1,100,mpi_comm_world,statut,ios)
522  maxdmin = max(maxdmin,rldmin)
523  enddo
524  endif
525  enddo
526  write(ig_lst_unit,'(a,e14.7)') "maximum localisation error of all double couple point sources = ",maxdmin
527  call flush(ig_lst_unit)
528  endif
529  else
530  if (ig_myrank == 0) write(ig_lst_unit,'(/,a)') "no double couple point source computed"
531  ig_ndcsource = 0
532  endif
533 
534  return
535 !***********************************************************************************************************************************************************************************
536  end subroutine init_double_couple_source
537 !***********************************************************************************************************************************************************************************
538 
539 !
540 !
547 !***********************************************************************************************************************************************************************************
548  subroutine init_single_force_source()
549 !***********************************************************************************************************************************************************************************
550 
551  use mpi
552 
553  use mod_global_variables, only :&
554  ig_lst_unit&
555  ,tg_sfsource&
556  ,ig_nsfsource&
557  ,ig_ncpu&
558  ,ig_myrank&
559  ,type_single_force_source&
560  ,cg_prefix&
561  ,rg_dt&
562  ,ig_ndt&
563  ,rg_sfsource_user_func&
564  ,ig_hexa_gll_glonum
565 
566  use mod_source_function , only : interp_linear
567 
568  use mod_receiver , only : search_closest_hexa_gll
569 
570  implicit none
571 
572  real, allocatable, dimension(:) :: dum1,dum2,dum3
573  real, dimension(ig_ncpu) :: dummy
574  real :: maxdmin
575  real :: rldmin
576 
577  integer, dimension(mpi_status_size) :: statut
578  integer, dimension(1) :: min_loc
579  integer, dimension(ig_ncpu) :: ilrcpu
580  integer :: i
581  integer :: j
582  integer :: ipf
583  integer :: ios
584  integer :: ios1
585  integer :: ilnpfo
586  integer :: itempo
587 
588  character(len=1) :: ctempo
589 
590  type(type_single_force_source), allocatable, dimension(:) :: tlpofo
591 
592  !
593  !
594  !*****************************************************************
595  !find element containing single force point source in cpu myrank
596  !*****************************************************************
597 
598  ilnpfo = 0
599 
600  open(unit=100,file=trim(cg_prefix)//".sfs",status='old',iostat=ios)
601 
602  if (ios == 0) then
603 
604  if (ig_myrank == 0) write(ig_lst_unit,'(a)') " "
605 
606  read(unit=100,fmt='(i10)') ig_nsfsource
607 
608  allocate(tlpofo(ig_nsfsource))
609 
610  do ipf = 1,ig_nsfsource
611 
612  !initialize tlpofo
613  tlpofo(ipf)%x = 0.0
614  tlpofo(ipf)%y = 0.0
615  tlpofo(ipf)%z = 0.0
616  tlpofo(ipf)%fac = 0.0
617  tlpofo(ipf)%rise_time = 0.0
618  tlpofo(ipf)%shift_time = 0.0
619  tlpofo(ipf)%var1 = 0.0
620  tlpofo(ipf)%dmin = 0.0
621  tlpofo(ipf)%icur = 0.0
622  tlpofo(ipf)%idir = 0
623  tlpofo(ipf)%cpu = 0
624  tlpofo(ipf)%iel = 0
625  tlpofo(ipf)%iequ = 0
626  tlpofo(ipf)%kgll = 0
627  tlpofo(ipf)%lgll = 0
628  tlpofo(ipf)%mgll = 0
629 
630  read(unit=100,fmt=*) tlpofo(ipf)%x,tlpofo(ipf)%y,tlpofo(ipf)%z,tlpofo(ipf)%fac,tlpofo(ipf)%idir,tlpofo(ipf)%shift_time,tlpofo(ipf)%rise_time,tlpofo(ipf)%icur
631 
632  if ( (ipf == 1) .and. (tlpofo(ipf)%icur == 10) ) then
633  itempo = 0
634  open(unit=99,file=trim(cg_prefix)//".sff")
635  do while (.true.)
636  read(unit= 99,fmt=*,iostat=ios1) ctempo
637  if (ios1 /= 0) exit
638  itempo = itempo + 1
639  enddo
640  rewind(99)
641  allocate(rg_sfsource_user_func(ig_ndt+1),dum3(ig_ndt+1),dum1(itempo),dum2(itempo))
642  do i = 1,itempo
643  read(unit=99,fmt=*) dum1(i),dum2(i)
644  enddo
645  close(99)
646  do i = 1,ig_ndt+1
647  dum3(i) = (i-1)*rg_dt
648  enddo
649  call interp_linear(1,itempo,dum1,dum2,ig_ndt+1,dum3,rg_sfsource_user_func)
650  if (ig_myrank == 0) then
651  open(unit=90,file=trim(cg_prefix)//".sff.interp")
652  do i = 1,ig_ndt+1
653  write(unit=90,fmt='(2(E14.7,1X))') dum3(i),rg_sfsource_user_func(i)
654  enddo
655  close(90)
656  endif
657  endif
658  !
659  !---->find the closest element and gll point to the source ipf
660  call search_closest_hexa_gll(tlpofo(ipf)%x,tlpofo(ipf)%y,tlpofo(ipf)%z &
661  ,tlpofo(ipf)%dmin,tlpofo(ipf)%kgll,tlpofo(ipf)%lgll,tlpofo(ipf)%mgll,tlpofo(ipf)%iel)
662  !
663  !->check which cpu has the smallest dmin. this cpu will compute this source and the other won't
664  call mpi_allgather(tlpofo(ipf)%dmin,1,mpi_real,dummy,1,mpi_real,mpi_comm_world,ios)
665  min_loc = minloc(dummy(1:ig_ncpu))
666  if ( (min_loc(1)-1) == ig_myrank) then
667  tlpofo(ipf)%cpu = ig_myrank
668  ilnpfo = ilnpfo + 1
669  else
670  tlpofo(ipf)%cpu = -1
671  endif
672  if (ig_myrank == 0) then
673  write(ig_lst_unit,'(a,i8,a,i0)') "single force point source ",ipf," computed by cpu ",min_loc(1)-1
674  endif
675  enddo
676  close(100)
677  !
678  !->transfer array tlpofo which knows all the sources to array tg_sfsource that needs only to know the sources of cpu 'ig_myrank'
679  if (ilnpfo > 0) then
680  allocate(tg_sfsource(ilnpfo))
681  j = 0
682  do ipf = 1,ig_nsfsource
683  if (tlpofo(ipf)%cpu == ig_myrank) then
684  j = j + 1
685  tg_sfsource(j) = tlpofo(ipf)
686  endif
687  enddo
688  ig_nsfsource = ilnpfo
689  else
690  ig_nsfsource = 0
691  endif
692  deallocate(tlpofo)
693  !
694  !->cpu 0 gather the number of single force point source (ig_nsfsource) of the other cpus
695  call mpi_gather(ig_nsfsource,1,mpi_integer,ilrcpu,1,mpi_integer,0,mpi_comm_world,ios)
696  !
697  do ipf = 1,ig_nsfsource
698  ! call pfoinf(tg_sfsource(ipf))
699  tg_sfsource(ipf)%iequ = ig_hexa_gll_glonum(tg_sfsource(ipf)%mgll,tg_sfsource(ipf)%lgll,tg_sfsource(ipf)%kgll,tg_sfsource(ipf)%iel)
700  if (ig_myrank /= 0) call mpi_send(tg_sfsource(ipf)%dmin,1,mpi_real,0,100,mpi_comm_world,ios)
701  enddo
702  !
703  !->cpu 0 gathers dmin of all cpus and write the maximum dmin in *.lst
704  if (ig_myrank == 0) then
705  maxdmin = 0.0
706  do i = 1,ig_ncpu
707  if (i == 1) then
708  do ipf = 1,ig_nsfsource
709  maxdmin = max(maxdmin,tg_sfsource(ipf)%dmin)
710  enddo
711  else
712  do ipf = 1,ilrcpu(i)
713  call mpi_recv(rldmin,1,mpi_real,i-1,100,mpi_comm_world,statut,ios)
714  maxdmin = max(maxdmin,rldmin)
715  enddo
716  endif
717  enddo
718  write(ig_lst_unit,'(a,e14.7)') "maximum localisation error of all single force point sources = ",maxdmin
719  call flush(ig_lst_unit)
720  endif
721  else
722  if (ig_myrank == 0) write(ig_lst_unit,'(/,a)') "no single force point source computed"
723  ig_nsfsource = 0
724  endif
725 
726  return
727 !***********************************************************************************************************************************************************************************
728  end subroutine init_single_force_source
729 !***********************************************************************************************************************************************************************************
730 
731 !
732 !
738 !***********************************************************************************************************************************************************************************
739  subroutine compute_source_local_coordinate(tlsrc)
740 !***********************************************************************************************************************************************************************************
741 
742  use mpi
743 
744  use mod_global_variables, only :&
745  rg_gll_abscissa&
746  ,ig_myrank&
747  ,type_double_couple_source
748 
749  use mod_jacobian , only : compute_hexa_jacobian
750 
751  use mod_coordinate , only : compute_hexa_point_coord
752 
753 
754  implicit none
755 
756  type(type_double_couple_source), intent(inout) :: tlsrc
757 
758  real :: eps
759  real :: dmin
760  real :: xisol
761  real :: etsol
762  real :: zesol
763  real :: newx
764  real :: newy
765  real :: newz
766  real :: dxidx
767  real :: dxidy
768  real :: dxidz
769  real :: detdx
770  real :: detdy
771  real :: detdz
772  real :: dzedx
773  real :: dzedy
774  real :: dzedz
775  real :: dx
776  real :: dy
777  real :: dz
778  real :: dxi
779  real :: det
780  real :: dze
781  real :: deter
782 
783  integer :: ihexa
784  integer :: iter
785 
786  integer, parameter :: iter_max=100
787 
788  !
789  !->initialize value
790  xisol = rg_gll_abscissa(tlsrc%mgll)
791  etsol = rg_gll_abscissa(tlsrc%lgll)
792  zesol = rg_gll_abscissa(tlsrc%kgll)
793  newx = tlsrc%x
794  newy = tlsrc%y
795  newz = tlsrc%z
796  dmin = tlsrc%dmin
797  ihexa = tlsrc%iel
798  eps = 1.0e-2
799  iter = 0
800 
801  !
802  !->solve the nonlinear system
803  do while(dmin > eps)
804 
805  call compute_hexa_jacobian(ihexa,xisol,etsol,zesol,dxidx,dxidy,dxidz,detdx,detdy,detdz,dzedx,dzedy,dzedz,deter)
806 
807  call compute_hexa_point_coord(ihexa,xisol,etsol,zesol,newx,newy,newz)
808 
809  dx = tlsrc%x - newx
810  dy = tlsrc%y - newy
811  dz = tlsrc%z - newz
812  dxi = dxidx*dx + dxidy*dy + dxidz*dz
813  det = detdx*dx + detdy*dy + detdz*dz
814  dze = dzedx*dx + dzedy*dy + dzedz*dz
815  xisol = xisol + dxi
816  etsol = etsol + det
817  zesol = zesol + dze
818  dmin = sqrt( (newx-tlsrc%x)**2 + (newy-tlsrc%y)**2 + (newz-tlsrc%z)**2 )
819 
820  iter = iter + 1
821 
822  if (mod(iter,iter_max) == 0) then
823  eps = eps*2.0
824  endif
825 
826  enddo
827 
828  tlsrc%xi = xisol
829  tlsrc%et = etsol
830  tlsrc%ze = zesol
831  tlsrc%dmin = dmin
832 
833  return
834 
835 !***********************************************************************************************************************************************************************************
836  end subroutine compute_source_local_coordinate
837 !***********************************************************************************************************************************************************************************
838 
839 
840 end module mod_source
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_receiver
This module contains subroutines to compute information about receivers and to write receivers' time ...
Definition: module_receiver.f90:126
mod_lagrange::lagrap
real function, public lagrap(i, x, n)
function to compute value of order Lagrange polynomial of the GLL node i at abscissa : ...
Definition: module_lagrange.f90:150
mod_receiver::search_closest_hexa_gll
subroutine, public search_closest_hexa_gll(x, y, z, dmin, kgll, lgll, mgll, iel)
This subroutine searches among all hexahedron elements in cpu myrank the closest GLL node to a point ...
Definition: module_receiver.f90:650
mod_source_function::interp_linear
subroutine, public interp_linear(dim_num, data_num, t_data, p_data, interp_num, t_interp, p_interp)
INTERP_LINEAR applies piecewise linear interpolation to data.
Definition: module_source_function.f90:394
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_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_source::compute_double_couple_source
subroutine, public compute_double_couple_source()
This subroutine pre-computes all double couple point sources defined by type mod_global_variables::ty...
Definition: module_source.f90:147
mod_source::init_double_couple_source
subroutine, public init_double_couple_source()
This subroutine reads all double couple point sources in file *.dcs; sets double couple point sources...
Definition: module_source.f90:330
mod_coordinate
This module contains subroutines to compute global -coordinates of a given point in the physical doma...
Definition: module_coordinate.f90:126
mod_coordinate::compute_hexa_point_coord
subroutine, public compute_hexa_point_coord(ihexa, xi, eta, zeta, x, y, z)
subroutine to compute -coordinates of a point knowing to which hexahedron element it belongs and its ...
Definition: module_coordinate.f90:151
mod_global_variables::type_single_force_source
type for single force point sources
Definition: module_global_variables.f90:661
mod_source
This module contains subroutines to compute information about sources.
Definition: module_source.f90:126
mod_source::init_single_force_source
subroutine, public init_single_force_source()
This subroutine reads all single force point sources in file *.sfs; determines to which cpu belong si...
Definition: module_source.f90:548
mod_source_function
This module contains functions and subroutines to compute tsource functions's time history...
Definition: module_source_function.f90:126
mod_global_variables::type_double_couple_source
type for double couple point sources
Definition: module_global_variables.f90:689
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_source::compute_source_local_coordinate
subroutine, private compute_source_local_coordinate(tlsrc)
This subroutine solves a nonlinear system by a gradient method to compute local coordinates of a poi...
Definition: module_source.f90:739