Merge branch 'devel' of https://github.com/geodynamics/specfem3d into…

… devel

CUBIT_GEOCUBIT/examples/homogeneous_halfspace/make_mesh.sh

@@ -16,4 +16,4 @@
 
 GEOCUBIT.py --build_volume --mesh --cfg=homogeneous_halfspace.cfg
 GEOCUBIT.py --collect --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH
-cp nummaterial_velocity_file.elastic_halfspace MESH/nummaterial_velocity_file
+cp nummaterial_velocity_file.elastic_halfspace MESH/nummaterial_velocity_file

CUBIT_GEOCUBIT/examples/homogeneous_halfspace_hex27/make_mesh_hex27.sh

@@ -16,4 +16,4 @@
 
 GEOCUBIT.py --build_volume --mesh --cfg=homogeneous_halfspace.cfg
 GEOCUBIT.py --collect --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH --hex27
-cp nummaterial_velocity_file.elastic_halfspace MESH/nummaterial_velocity_file
+cp nummaterial_velocity_file.elastic_halfspace MESH/nummaterial_velocity_file

CUBIT_GEOCUBIT/examples/layered_halfspace/make_mesh_notripling.sh

@@ -16,4 +16,4 @@
 
 GEOCUBIT.py --build_volume --mesh --cfg='./notripling/layered_halfspace_notripling.cfg'
 GEOCUBIT.py --collect --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH
-cp nummaterial_velocity_file.notripling MESH/nummaterial_velocity_file
+cp nummaterial_velocity_file.notripling MESH/nummaterial_velocity_file

CUBIT_GEOCUBIT/examples/layered_halfspace/make_mesh_tripling.sh

@@ -16,4 +16,4 @@
 
 GEOCUBIT.py --build_volume --mesh --cfg='./notripling/layered_halfspace_tripling.cfg'
 GEOCUBIT.py --collect --meshfiles=MESH_GEOCUBIT/mesh_vol_0.e --export2SPECFEM3D --SEMoutput=MESH
-cp nummaterial_velocity_file.tripling MESH/nummaterial_velocity_file
+cp nummaterial_velocity_file.tripling MESH/nummaterial_velocity_file

CUBIT_GEOCUBIT/examples/tomographic_model/tomography/create_tomography_model_file.sh

@@ -1,21 +1,21 @@
 #!/bin/sh
 
-# creates a simple, formatted tomography model with a constant velocity gradient 
+# creates a simple, formatted tomography model with a constant velocity gradient
 # for a block model with dimensions 134000 x 134000 x 60000
 
 # origin points
 ORIG_X=0.
-ORIG_Y=0. 
-ORIG_Z=0. 
+ORIG_Y=0.
+ORIG_Z=0.
 
 # end points
-END_X=134000. 
-END_Y=134000. 
+END_X=134000.
+END_Y=134000.
 END_Z=-60000.  # depth in negative z-direction
 
 # spacing of given tomography points
-SPACING_X=2000. 
-SPACING_Y=2000. 
+SPACING_X=2000.
+SPACING_Y=2000.
 SPACING_Z=-2000.
 
 # number of cell increments
@@ -24,11 +24,11 @@ NY=68
 NZ=31
 
 # min/max values
-VP_MIN=2500. 
-VP_MAX=8500. 
+VP_MIN=2500.
+VP_MAX=8500.
 VS_MIN=1500.
-VS_MAX=7500. 
-RHO_MIN=1500. 
+VS_MAX=7500.
+RHO_MIN=1500.
 RHO_MAX=1500.
 
 

CUBIT_GEOCUBIT/geocubitlib/save_fault_nodes_elements.py

@@ -0,0 +1,117 @@
+#!python
+# Opening fault cracks
+# Input : surface up and down.
+import cubit 
+
+class fault_input:
+   def __init__(self,id,surface_u,surface_d): 
+       self.id = id
+       self.surface_u = surface_u
+       self.surface_d = surface_d
+       self.name = 'MESH/fault_file_'+str(id)+'.dat'
+
+       quads_Aup,quads_Adp = save_cracks(self.name,self.surface_u,self.surface_d)
+       #Unpacking list.
+       quads_Au=unpack_list(quads_Aup)
+       quads_Ad=unpack_list(quads_Adp)
+
+       print 'len(Au):',len(quads_Au)
+       print 'len(Ad):',len(quads_Ad)
+
+       if not (len(quads_Au)==len(quads_Ad)):
+           print 'Number of elements for each fauld side up and down do not concide'
+           sys.exit('goodbye')
+
+       save_elements_nodes(self.name,quads_Au,quads_Ad)
+
+
+def save_cracks(name,list_surface_up,list_surface_down):
+   quads_fault_up = [] 
+   quads_fault_down = []
+   for surface in list_surface_up   :
+       quads_fault = cubit.get_surface_quads(surface)
+       quads_fault_up.append(quads_fault)
+   for surface in list_surface_down :
+       quads_fault = cubit.get_surface_quads(surface)
+       quads_fault_down.append(quads_fault)
+    # TO DO : stop python properly in case fault nodes at both sides
+    #         do not match.
+    #   if len(quads_fault_u) != len(quads_fault_d): stop
+    #
+    # SAVING FAULT ELEMENTS AND NODES
+   return quads_fault_up,quads_fault_down
+
+def unpack_list(fault_list):
+    list_fault = []
+    for i in range(0,len(fault_list)):
+        el=list(fault_list[i])
+        for j in el:
+            list_fault.append(j)
+    return list_fault 
+
+def save_elements_nodes(name,quads_fault_u,quads_fault_d):
+   fault_file = open(name,'w')
+   txt =''
+   list_hex=cubit.parse_cubit_list('hex','all')
+   txt='%10i %10i\n' % (len(quads_fault_u),len(quads_fault_d))
+   fault_file.write(txt)
+
+   dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u)) 
+   dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d)) 
+
+   # FAULT SIDE DOWN
+  # fault_file.write('upsurface')
+   for h in list_hex: 
+       faces = cubit.get_sub_elements('hex',h,2)  
+       for f in faces:
+           if dic_quads_fault_d.has_key(f):
+	      cubit.silent_cmd('group "nf" add Node in face '+str(f))
+	      group1 = cubit.get_id_from_name("nf")
+	      nodes = cubit.get_group_nodes(group1)
+	      cubit.silent_cmd('del group '+ str(group1))
+#              nodes=cubit.get_connectivity('Face',f)
+#              print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+#              txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+#                                                nodes[1],nodes[2],nodes[3])
+              ngnod2d = len(nodes)
+	      if ngnod2d == 9:
+#kangchen added               txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+#                                                nodes[1],nodes[2],nodes[3])
+                      txt='%10i %10i %10i %10i %10i %10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+                                                nodes[1],nodes[2],nodes[3],nodes[4],nodes[5],nodes[6],nodes[7],nodes[8])
+	      else:
+		      txt='%10i %10i %10i %10i %10i \n' % (h,nodes[0],\
+				                                                      nodes[1],nodes[2],nodes[3])
+
+
+
+
+              fault_file.write(txt)
+
+   # FAULT SIDE UP
+ #  fault_file.write('downsurface')
+   for h in list_hex: 
+       faces = cubit.get_sub_elements('hex',h,2)  
+       for f in faces:
+           if dic_quads_fault_u.has_key(f): 
+       	      cubit.silent_cmd('group "nf" add Node in face '+str(f))
+	      group1 = cubit.get_id_from_name("nf")
+	      nodes = cubit.get_group_nodes(group1)
+	      cubit.cmd('del group '+ str(group1))
+	      ngnod2d=len(nodes)
+	      if ngnod2d == 9:
+#	      nodes=cubit.get_connectivity('Face',f)
+#             print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
+#kangchen added               txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+#                                                nodes[1],nodes[2],nodes[3])
+                      txt='%10i %10i %10i %10i %10i %10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+                                                nodes[1],nodes[2],nodes[3],nodes[4],nodes[5],nodes[6],nodes[7],nodes[8])
+	      else:
+		      txt='%10i %10i %10i %10i %10i \n' % (h,nodes[0],\
+				                                                      nodes[1],nodes[2],nodes[3])
+
+
+
+              fault_file.write(txt)
+
+   fault_file.close()

EXAMPLES/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/DATA/Par_file

@@ -119,10 +119,12 @@ USE_FORCE_POINT_SOURCE          = .false.
 # to represent a (tilted) FORCESOLUTION force point source or a CMTSOLUTION moment-tensor source.
 USE_RICKER_TIME_FUNCTION        = .false.
 
-# decide if we save displacement, velocity and/or acceleration in forward runs (they can be set to true simultaneously)
+# decide if we save displacement, velocity, acceleration and/or pressure in forward runs (they can be set to true simultaneously)
+# currently pressure seismograms are implemented in acoustic (i.e. fluid) elements only
 SAVE_SEISMOGRAMS_DISPLACEMENT   = .true.
 SAVE_SEISMOGRAMS_VELOCITY       = .true.
 SAVE_SEISMOGRAMS_ACCELERATION   = .true.
+SAVE_SEISMOGRAMS_PRESSURE       = .false.  # currently implemented in acoustic (i.e. fluid) elements only
 
 # save seismograms in binary or ASCII format (binary is smaller but may not be portable between machines)
 USE_BINARY_FOR_SEISMOGRAMS      = .false.
@@ -137,6 +139,49 @@ WRITE_SEISMOGRAMS_BY_MASTER     = .false.
 # to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance
 SAVE_ALL_SEISMOS_IN_ONE_FILE    = .false.
 
+# use a trick to increase accuracy of pressure seismograms in fluid (acoustic) elements:
+# use the second derivative of the source for the source time function instead of the source itself,
+# and then record -potential_acoustic() as pressure seismograms instead of -potential_dot_dot_acoustic();
+# this is mathematically equivalent, but numerically significantly more accurate because in the explicit
+# Newmark time scheme acceleration is accurate at zeroth order while displacement is accurate at second order,
+# thus in fluid elements potential_dot_dot_acoustic() is accurate at zeroth order while potential_acoustic()
+# is accurate at second order and thus contains significantly less numerical noise.
+USE_TRICK_FOR_BETTER_PRESSURE   = .false.
+
+#
+# source encoding
+#
+# determines source encoding factor +/-1 depending on sign of moment tensor
+# (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (6), WCC177-WCC188.)
+USE_SOURCE_ENCODING             = .false.
+
+#
+# total energy calculation
+#
+# to plot total energy curves, for instance to monitor how CPML absorbing layers behave;
+# should be turned OFF in most cases because expensive
+OUTPUT_ENERGY                   = .false.
+# every how many time steps we compute energy (which is expensive to compute)
+NTSTEP_BETWEEN_OUTPUT_ENERGY    = 10
+
+#
+# adjoint kernel outputs
+#
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL                  = .false.
+
+# compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
+SAVE_TRANSVERSE_KL              = .false.
+
+# outputs approximate Hessian for preconditioning
+APPROXIMATE_HESS_KL             = .false.
+
+# save Moho mesh and compute Moho boundary kernels
+SAVE_MOHO_MESH                  = .false.
+
 # print source time function
 PRINT_SOURCE_TIME_FUNCTION      = .false.
 

EXAMPLES/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/DATA/Par_file

@@ -119,10 +119,12 @@ USE_FORCE_POINT_SOURCE          = .false.
 # to represent a (tilted) FORCESOLUTION force point source or a CMTSOLUTION moment-tensor source.
 USE_RICKER_TIME_FUNCTION        = .false.
 
-# decide if we save displacement, velocity and/or acceleration in forward runs (they can be set to true simultaneously)
+# decide if we save displacement, velocity, acceleration and/or pressure in forward runs (they can be set to true simultaneously)
+# currently pressure seismograms are implemented in acoustic (i.e. fluid) elements only
 SAVE_SEISMOGRAMS_DISPLACEMENT   = .true.
 SAVE_SEISMOGRAMS_VELOCITY       = .true.
 SAVE_SEISMOGRAMS_ACCELERATION   = .true.
+SAVE_SEISMOGRAMS_PRESSURE       = .false.  # currently implemented in acoustic (i.e. fluid) elements only
 
 # save seismograms in binary or ASCII format (binary is smaller but may not be portable between machines)
 USE_BINARY_FOR_SEISMOGRAMS      = .false.
@@ -137,6 +139,49 @@ WRITE_SEISMOGRAMS_BY_MASTER     = .false.
 # to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance
 SAVE_ALL_SEISMOS_IN_ONE_FILE    = .false.
 
+# use a trick to increase accuracy of pressure seismograms in fluid (acoustic) elements:
+# use the second derivative of the source for the source time function instead of the source itself,
+# and then record -potential_acoustic() as pressure seismograms instead of -potential_dot_dot_acoustic();
+# this is mathematically equivalent, but numerically significantly more accurate because in the explicit
+# Newmark time scheme acceleration is accurate at zeroth order while displacement is accurate at second order,
+# thus in fluid elements potential_dot_dot_acoustic() is accurate at zeroth order while potential_acoustic()
+# is accurate at second order and thus contains significantly less numerical noise.
+USE_TRICK_FOR_BETTER_PRESSURE   = .false.
+
+#
+# source encoding
+#
+# determines source encoding factor +/-1 depending on sign of moment tensor
+# (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (6), WCC177-WCC188.)
+USE_SOURCE_ENCODING             = .false.
+
+#
+# total energy calculation
+#
+# to plot total energy curves, for instance to monitor how CPML absorbing layers behave;
+# should be turned OFF in most cases because expensive
+OUTPUT_ENERGY                   = .false.
+# every how many time steps we compute energy (which is expensive to compute)
+NTSTEP_BETWEEN_OUTPUT_ENERGY    = 10
+
+#
+# adjoint kernel outputs
+#
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL                  = .false.
+
+# compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
+SAVE_TRANSVERSE_KL              = .false.
+
+# outputs approximate Hessian for preconditioning
+APPROXIMATE_HESS_KL             = .false.
+
+# save Moho mesh and compute Moho boundary kernels
+SAVE_MOHO_MESH                  = .false.
+
 # print source time function
 PRINT_SOURCE_TIME_FUNCTION      = .false.
 

EXAMPLES/CPML_examples/homogeneous_halfspace_HEX8_acoustic_absorbing_CPML_5sides/DATA/Par_file

@@ -119,10 +119,12 @@ USE_FORCE_POINT_SOURCE          = .false.
 # to represent a (tilted) FORCESOLUTION force point source or a CMTSOLUTION moment-tensor source.
 USE_RICKER_TIME_FUNCTION        = .false.
 
-# decide if we save displacement, velocity and/or acceleration in forward runs (they can be set to true simultaneously)
+# decide if we save displacement, velocity, acceleration and/or pressure in forward runs (they can be set to true simultaneously)
+# currently pressure seismograms are implemented in acoustic (i.e. fluid) elements only
 SAVE_SEISMOGRAMS_DISPLACEMENT   = .true.
 SAVE_SEISMOGRAMS_VELOCITY       = .true.
 SAVE_SEISMOGRAMS_ACCELERATION   = .true.
+SAVE_SEISMOGRAMS_PRESSURE       = .false.  # currently implemented in acoustic (i.e. fluid) elements only
 
 # save seismograms in binary or ASCII format (binary is smaller but may not be portable between machines)
 USE_BINARY_FOR_SEISMOGRAMS      = .false.
@@ -137,6 +139,49 @@ WRITE_SEISMOGRAMS_BY_MASTER     = .false.
 # to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance
 SAVE_ALL_SEISMOS_IN_ONE_FILE    = .false.
 
+# use a trick to increase accuracy of pressure seismograms in fluid (acoustic) elements:
+# use the second derivative of the source for the source time function instead of the source itself,
+# and then record -potential_acoustic() as pressure seismograms instead of -potential_dot_dot_acoustic();
+# this is mathematically equivalent, but numerically significantly more accurate because in the explicit
+# Newmark time scheme acceleration is accurate at zeroth order while displacement is accurate at second order,
+# thus in fluid elements potential_dot_dot_acoustic() is accurate at zeroth order while potential_acoustic()
+# is accurate at second order and thus contains significantly less numerical noise.
+USE_TRICK_FOR_BETTER_PRESSURE   = .false.
+
+#
+# source encoding
+#
+# determines source encoding factor +/-1 depending on sign of moment tensor
+# (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (6), WCC177-WCC188.)
+USE_SOURCE_ENCODING             = .false.
+
+#
+# total energy calculation
+#
+# to plot total energy curves, for instance to monitor how CPML absorbing layers behave;
+# should be turned OFF in most cases because expensive
+OUTPUT_ENERGY                   = .false.
+# every how many time steps we compute energy (which is expensive to compute)
+NTSTEP_BETWEEN_OUTPUT_ENERGY    = 10
+
+#
+# adjoint kernel outputs
+#
+# this parameter must be set to .true. to compute anisotropic kernels
+# in crust and mantle (related to the 21 Cij in geographical coordinates)
+# default is .false. to compute isotropic kernels (related to alpha and beta)
+ANISOTROPIC_KL                  = .false.
+
+# compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+# rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
+SAVE_TRANSVERSE_KL              = .false.
+
+# outputs approximate Hessian for preconditioning
+APPROXIMATE_HESS_KL             = .false.
+
+# save Moho mesh and compute Moho boundary kernels
+SAVE_MOHO_MESH                  = .false.
+
 # print source time function
 PRINT_SOURCE_TIME_FUNCTION      = .false.