Implement Anisotropic Crystal example

CMakeLists.txt

@@ -459,6 +459,7 @@ target_link_libraries(
 add_library(
         frechet_derivatives
         src/frechet_derivatives/impl/element_kernel/elastic_isotropic.cpp
+        src/frechet_derivatives/impl/element_kernel/elastic_anisotropic.cpp
         src/frechet_derivatives/impl/element_kernel/acoustic_isotropic.cpp
         src/frechet_derivatives/impl/frechet_element.cpp
         src/frechet_derivatives/frechet_derivatives.cpp

examples/CMakeLists.txt

@@ -1,5 +1,6 @@
 cmake_minimum_required(VERSION 3.17.5)
 
+add_subdirectory(anisotropic-crystal)
 add_subdirectory(homogeneous-medium-flat-topography)
 add_subdirectory(fluid-solid-interface)
 add_subdirectory(Tromp_2005)

examples/Tromp_2005/CMakeFiles/Snakefile.in

@@ -49,14 +49,16 @@ rule forward_configuration:
                             }
                         }
                     }
-                }
+        }
 
         with open(input.config, "r") as f:
             config = yaml.safe_load(f)
 
             ## Add forward node to the simulation setup
         config["parameters"]["simulation-setup"].update(forward)
         config["parameters"]["databases"]["source-file"] = "@CMAKE_SOURCE_DIR@/examples/Tromp_2005/forward_source.yaml"
+        config["parameters"]["receivers"]["stations-file"] = "@CMAKE_SOURCE_DIR@/examples/Tromp_2005/OUTPUT_FILES/STATIONS"
+        config["parameters"]["databases"]["mesh-database"] = "@CMAKE_SOURCE_DIR@/examples/Tromp_2005/OUTPUT_FILES/database.bin"
 
         with open(output.config, "w") as f:
             yaml.safe_dump(config, f)
@@ -90,8 +92,8 @@ rule forward_simulation:
         runtime=10,
     shell:
         """
-            module purge
-            module load boost/1.73.0
+            # module purge
+            # module load boost/1.73.0
             mkdir -p OUTPUT_FILES/results
             echo "Hostname: $(hostname)" > output.log
             {SPECFEM_BIN} -p {input.config} >> output.log
@@ -157,6 +159,7 @@ rule adjoint_configuration:
             ## Add adjoint node to the simulation setup
         config["parameters"]["simulation-setup"].update(adjoint)
         config["parameters"]["databases"]["source-file"] = "@CMAKE_SOURCE_DIR@/examples/Tromp_2005/adjoint_source.yaml"
+        config["parameters"]["databases"]["mesh-database"] = "@CMAKE_SOURCE_DIR@/examples/Tromp_2005/OUTPUT_FILES/database.bin"
 
         with open(output.config, "w") as f:
             yaml.safe_dump(config, f)
@@ -212,8 +215,8 @@ rule adjoint_simulation:
         runtime=10,
     shell:
         """
-            module purge
-            module load boost/1.73.0
+            # module purge
+            # module load boost/1.73.0
             mkdir -p OUTPUT_FILES/results
             echo "Hostname: $(hostname)" > output.log
             {SPECFEM_BIN} -p {input.config} >> output.log
@@ -228,7 +231,8 @@ rule plot_kernels:
         plot="OUTPUT_FILES/Kernels.png",
     run:
         from plot import plot_kernels
-
+        import matplotlib
+        matplotlib.use('agg')
         plot_kernels(input.kernels, output.plot)
 
 

examples/Tromp_2005/templates/specfem_config.yaml

@@ -25,7 +25,7 @@ parameters:
           t0: 8.0
 
   receivers:
-    stations-file: "/scratch/gpfs/rk9481/specfem2d_kokkos/examples/Tromp_2005/OUTPUT_FILES/STATIONS"
+    stations-file: placeholder/for/path/to/STATIONS
     angle: 0.0
     seismogram-type:
       - displacement
@@ -36,5 +36,5 @@ parameters:
     number-of-runs: 1
 
   databases:
-    mesh-database: "/scratch/gpfs/rk9481/specfem2d_kokkos/examples/Tromp_2005/OUTPUT_FILES/database.bin"
-    source-file: "/scratch/gpfs/rk9481/specfem2d_kokkos/examples/Tromp_2005/source.yaml"
+    mesh-database: placeholder/for/path/to/database.bin
+    source-file: placeholder/for/path/to/source.yaml

examples/anisotropic-crystal/.gitignore

@@ -0,0 +1,2 @@
+Snakefile
+specfem_config.yaml

examples/anisotropic-crystal/CMakeFiles/Par_File.in

@@ -0,0 +1,135 @@
+#-----------------------------------------------------------
+#
+# Simulation input parameters
+#
+#-----------------------------------------------------------
+
+# title of job
+title                           = Anisotropic Crystal
+# parameters concerning partitioning
+NPROC                           = 1              # number of processes
+
+# Output folder to store mesh related files
+OUTPUT_FILES                   = @CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES
+
+#-----------------------------------------------------------
+#
+# Mesh
+#
+#-----------------------------------------------------------
+
+# Partitioning algorithm for decompose_mesh
+PARTITIONING_TYPE               = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+
+# number of control nodes per element (4 or 9)
+NGNOD                           = 9
+
+# location to store the mesh
+database_filename               = @CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/database.bin
+
+#-----------------------------------------------------------
+#
+# Receivers
+#
+#-----------------------------------------------------------
+
+
+
+# use an existing STATION file found in ./DATA or create a new one from the receiver positions below in this Par_file
+use_existing_STATIONS           = .false.
+
+# number of receiver sets (i.e. number of receiver lines to create below)
+nreceiversets                   = 1
+
+# orientation
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
+nrec                            = 50             # number of receivers
+xdeb                            = 0.05           # first receiver x in meters
+zdeb                            = 0.2640         # first receiver z in meters
+xfin                            = 0.28           # last receiver x in meters (ignored if only one receiver)
+zfin                            = 0.2640         # last receiver z in meters (ignored if only one receiver)
+record_at_surface_same_vertical = .false.        # receivers inside the medium or at the surface
+
+# filename to store stations file
+stations_filename              = @CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/STATIONS
+
+#-----------------------------------------------------------
+#
+# Velocity and density models
+#
+#-----------------------------------------------------------
+
+# number of model materials
+nbmodels                        = 1
+# available material types (see user manual for more information)
+#   acoustic:    model_number 1 rho Vp 0  0 0 QKappa Qmu 0 0 0 0 0 0
+#   elastic:     model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0
+#   anistoropic: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 0 0 0
+#   poroelastic: model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu
+#   tomo:        model_number -1 0 9999 9999 A 0 0 9999 9999 0 0 0 0 0
+1 2 7100. 16.5d10 5.d10 0 6.2d10 0 3.96d10 0 0 0 0 0 0
+
+# external tomography file
+TOMOGRAPHY_FILE                 = ./DATA/tomo_file.xyz
+
+# use an external mesh created by an external meshing tool or use the internal mesher
+read_external_mesh              = .false.
+
+#-----------------------------------------------------------
+#
+# PARAMETERS FOR EXTERNAL MESHING
+#
+#-----------------------------------------------------------
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+axial_elements_file             = ./DATA/axial_elements_file   # file containing the axial elements if AXISYM is true
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+acoustic_forcing_surface_file   = ./DATA/MSH/Surf_acforcing_Bottom_enforcing_mesh   # file containing the acoustic forcing surface
+absorbing_cpml_file             = ./DATA/absorbing_cpml_file   # file containing the CPML element numbers
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes  # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------
+#
+# PARAMETERS FOR INTERNAL MESHING
+#
+#-----------------------------------------------------------
+
+# file containing interfaces for internal mesh
+interfacesfile                  = @CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/topoaniso.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 0.33           # abscissa of right side of the model
+nx                              = 60             # number of elements along X
+
+# Stacey ABC
+STACEY_ABSORBING_CONDITIONS     = .false.
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .false.
+absorbright                     = .false.
+absorbtop                       = .false.
+absorbleft                      = .false.
+
+# define the different regions of the model in the (nx,nz) spectral-element mesh
+nbregions                       = 1              # then set below the different regions and model number for each region
+# format of each line: nxmin nxmax nzmin nzmax material_number
+1 60 1   60 1
+
+#-----------------------------------------------------------
+#
+# Display parameters
+#
+#-----------------------------------------------------------
+
+# meshing output
+output_grid_Gnuplot             = .false.        # generate a GNUPLOT file containing the grid, and a script to plot it
+output_grid_ASCII               = .false.        # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not