Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implements the anisotropy/isotropic model example from specfem2d #326

Closed
wants to merge 7 commits into from
Closed

Implements the anisotropy/isotropic model example from specfem2d #326

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
71d373a
added anisotropy crystal
lsawade Dec 17, 2024
0696797
Merge branch 'issue-228-anisotropy' into issue-306
lsawade Dec 18, 2024
adea203
The example runs but the traces are nans
lsawade Dec 18, 2024
9d2c1ce
The fixed a few issues: DT, NSTEPs, and some file writing. The seismo…
lsawade Dec 18, 2024
9067670
Merge branch 'issue-228-anisotropy' into issue-306
lsawade Dec 18, 2024
6c19280
Added basic structure for the anisotropic-isotropic-model example
lsawade Dec 19, 2024
ab0d526
Added the example for the anisotropic isotropic model. Updated the pl…
lsawade Dec 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions examples/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
cmake_minimum_required(VERSION 3.17.5)

add_subdirectory(anisotropic-crystal)
add_subdirectory(anisotropic-isotropic-model)
add_subdirectory(homogeneous-medium-flat-topography)
add_subdirectory(fluid-solid-interface)
add_subdirectory(Tromp_2005)
6 changes: 6 additions & 0 deletions examples/anisotropic-crystal/.gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
Snakefile
specfem_config.yaml
Par_file
output.log
OUTPUT_FILES
.snakemake
135 changes: 135 additions & 0 deletions examples/anisotropic-crystal/CMakeFiles/Par_File.in
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
215 changes: 215 additions & 0 deletions examples/anisotropic-crystal/CMakeFiles/Snakefile.in
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
SPECFEM_BIN = "specfem2d"
MESHFEM_BIN = "xmeshfem2D"


rule all:
input:
plot="OUTPUT_FILES/results/plot.png",
localrule: True


rule generate_mesh:
input:
"Par_File",
output:
database="@CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/database.bin",
stations="@CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/STATIONS",
localrule: True
shell:
"""
mkdir -p OUTPUT_FILES
{MESHFEM_BIN} -p {input}
"""


rule run_solver:
input:
database="@CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/database.bin",
stations="@CMAKE_SOURCE_DIR@/examples/anisotropic-crystal/OUTPUT_FILES/STATIONS",
source="sources.yaml",
config="specfem_config.yaml",
output:
seismograms=expand(
"OUTPUT_FILES/results/{station_name}{network_name}{component}.semd",
station_name=[
"S0001",
"S0002",
"S0003",
"S0004",
"S0005",
"S0006",
"S0007",
"S0008",
"S0009",
"S0010",
"S0011",
"S0012",
"S0013",
"S0014",
"S0015",
"S0016",
"S0017",
"S0018",
"S0019",
"S0020",
"S0021",
"S0022",
"S0023",
"S0024",
"S0025",
"S0026",
"S0027",
"S0028",
"S0029",
"S0030",
"S0031",
"S0032",
"S0033",
"S0034",
"S0035",
"S0036",
"S0037",
"S0038",
"S0039",
"S0040",
"S0041",
"S0042",
"S0043",
"S0044",
"S0045",
"S0046",
"S0047",
"S0048",
"S0049",
"S0050",
],
network_name=["AA"],
component=["BXX", "BXZ"],
),
resources:
nodes=1,
tasks=1,
cpus_per_task=1,
runtime=10,
shell:
"""
# 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
"""


rule plot_seismogram:
input:
trace_files=expand(
"OUTPUT_FILES/results/{station_name}{network_name}{component}.semd",
station_name=[
"S0001",
"S0002",
"S0003",
"S0004",
"S0005",
"S0006",
"S0007",
"S0008",
"S0009",
"S0010",
"S0011",
"S0012",
"S0013",
"S0014",
"S0015",
"S0016",
"S0017",
"S0018",
"S0019",
"S0020",
"S0021",
"S0022",
"S0023",
"S0024",
"S0025",
"S0026",
"S0027",
"S0028",
"S0029",
"S0030",
"S0031",
"S0032",
"S0033",
"S0034",
"S0035",
"S0036",
"S0037",
"S0038",
"S0039",
"S0040",
"S0041",
"S0042",
"S0043",
"S0044",
"S0045",
"S0046",
"S0047",
"S0048",
"S0049",
"S0050",
],
network_name=["AA"],
component=["BXX", "BXZ"],
),
output:
traces="OUTPUT_FILES/results/plot.png",
localrule: True
run:
import glob
import os
import numpy as np
import obspy
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("agg")

def get_traces(directory):
traces = []
files = glob.glob(directory + "/*.sem*")
## iterate over all seismograms
for filename in files:
station_name = os.path.splitext(filename)[0]
station_name = station_name.split("/")[-1]
trace = np.loadtxt(filename, delimiter=" ")
starttime = trace[0, 0]
dt = trace[1, 0] - trace[0, 0]
traces.append(
obspy.Trace(
trace[:, 1],
{"network": station_name, "starttime": starttime, "delta": dt},
)
)

stream = obspy.Stream(traces)

return stream


stream = get_traces("OUTPUT_FILES/results")

N_traces = len(stream)
Amax = np.max(stream.max())
plt.figure(figsize=(6, 8))

for i, trace in enumerate(stream):
plt.plot(trace.times('matplotlib'), trace.data + i * Amax,
label=f"{trace.stats.network}.{trace.stats.station}")

plt.savefig(output.traces, dpi=300)
plt.close('all')


rule clean:
shell:
"""
rm -rf OUTPUT_FILES
"""
Loading
Loading