adds Moon example

Author: danielpeter

EXAMPLES/real_world/Moon_Apollo_impact/DATA/FORCESOLUTION

@@ -0,0 +1,11 @@
+APOLLO 16 S-IVB BOOSTER IMPACT
+time shift:      0.0000
+hdurorf0:        0.1
+latorUTM:        1.9210 
+longorUTM:     335.3770
+depth:           0.0
+source time function:            0
+factor force source:             3.6d8
+component dir vect source E:     -0.1845634
+component dir vect source N:     -0.0484193
+component dir vect source Z_UP:  -0.9816272

EXAMPLES/real_world/Moon_Apollo_impact/DATA/Par_file

@@ -0,0 +1,2 @@
+S-IVB      AP16       1.9210       335.3770       0.0    0.0
+S12        XA        -3.01084     -23.42456       0.0    0.0

EXAMPLES/real_world/Moon_Apollo_impact/DATA/STATIONS

@@ -0,0 +1,112 @@
+#-----------------------------------------------------------
+#
+# Meshing input parameters
+#
+#-----------------------------------------------------------
+
+# coordinates of mesh block in latitude/longitude and depth in km
+LATITUDE_MIN                    = -5.0
+LATITUDE_MAX                    = 5.0
+LONGITUDE_MIN                   = -30.0
+LONGITUDE_MAX                   = -20.0
+DEPTH_BLOCK_KM                  = 100.d0
+UTM_PROJECTION_ZONE             = 20
+SUPPRESS_UTM_PROJECTION         = .false.
+
+# file that contains the interfaces of the model / mesh
+INTERFACES_FILE                 = interfaces.dat
+
+# file that contains the cavity
+CAVITY_FILE                     = dummy 
+
+# number of elements at the surface along edges of the mesh at the surface
+# (must be 8 * multiple of NPROC below if mesh is not regular and contains mesh doublings)
+# (must be multiple of NPROC below if mesh is regular)
+NEX_XI                          = 128 
+NEX_ETA                         = 128
+
+# number of MPI processors along xi and eta (can be different)
+NPROC_XI                        = 4 
+NPROC_ETA                       = 2
+
+#-----------------------------------------------------------
+#
+# Doubling layers
+#
+#-----------------------------------------------------------
+
+# Regular/irregular mesh
+USE_REGULAR_MESH                = .false.
+# Only for irregular meshes, number of doubling layers and their position
+NDOUBLINGS                      = 4 
+# NZ_DOUBLING_1 is the parameter to set up if there is only one doubling layer
+# (more doubling entries can be added if needed to match NDOUBLINGS value)
+NZ_DOUBLING_1                   = 4
+NZ_DOUBLING_2                   = 6 
+NZ_DOUBLING_3                   = 12
+NZ_DOUBLING_4                   = 14
+
+#-----------------------------------------------------------
+#
+# Visualization
+#
+#-----------------------------------------------------------
+
+# create mesh files for visualisation or further checking
+CREATE_ABAQUS_FILES             = .false.
+CREATE_DX_FILES                 = .false.
+CREATE_VTK_FILES                = .true.
+
+# stores mesh files as cubit-exported files into directory MESH/ (for single process run)
+SAVE_MESH_AS_CUBIT              = .false.
+
+# path to store the databases files
+LOCAL_PATH                      = ./DATABASES_MPI
+
+#-----------------------------------------------------------
+#
+# CPML
+#
+#-----------------------------------------------------------
+
+# CPML perfectly matched absorbing layers
+THICKNESS_OF_X_PML              = 0.05d0
+THICKNESS_OF_Y_PML              = 0.05d0
+THICKNESS_OF_Z_PML              = 0.1d0
+
+# add PML layers as extra outer mesh layers
+ADD_PML_AS_EXTRA_MESH_LAYERS    = .true.
+NUMBER_OF_PML_LAYERS_TO_ADD     = 2
+
+#-----------------------------------------------------------
+#
+# Domain materials
+#
+#-----------------------------------------------------------
+
+# number of materials
+NMATERIALS                      = 3
+# define the different materials in the model as:
+# #material_id  #rho  #vp  #vs  #Q_Kappa  #Q_mu  #anisotropy_flag  #domain_id
+#     Q_Kappa          : Q_Kappa attenuation quality factor
+#     Q_mu             : Q_mu attenuation quality factor
+#     anisotropy_flag  : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
+#     domain_id        : 1 = acoustic / 2 = elastic / 3 = poroelastic
+1  3312  7540 4340  9999. 6750.0  0  2      # moon mantle
+2  2762  3200 1800  9999. 6750.0  0  2      # moon upper crust
+3  2600  1000  500  9999. 6750.0  0  2      # moon regolith
+
+#-----------------------------------------------------------
+#
+# Domain regions
+#
+#-----------------------------------------------------------
+
+# number of regions
+NREGIONS                        = 3
+# define the different regions of the model as :
+#NEX_XI_BEGIN  #NEX_XI_END  #NEX_ETA_BEGIN  #NEX_ETA_END  #NZ_BEGIN #NZ_END  #material_id
+1              -1            1              -1            1         4        1
+1              -1            1              -1            5         14       2
+1              -1            1              -1            15        15       3
+

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/Mesh_Par_file

@@ -0,0 +1,4 @@
+-8000.0
+-8000.0
+-8000.0
+-8000.0

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/interface.inter.dat

@@ -0,0 +1,4 @@
+-30800.0
+-30800.0
+-30800.0
+-30800.0

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/interface.moho.dat

@@ -0,0 +1,34 @@
+# number of interfaces
+4  
+
+# We describe each interface below, structured as a 2D-grid, with several parameters : 
+# number of points along XI and ETA, minimal XI ETA coordinates 
+# and spacing between points which must be constant.
+# Then the records contain the Z coordinates of the NXI x NETA points.
+#
+# interface moho
+# SUPPRESS_UTM_PROJECTION  NXI  NETA LONG_MIN   LAT_MIN    SPACING_XI SPACING_ETA
+.true.                     2    2    68634.88   155340.91  0.0045     -0.0045
+interface.moho.dat
+
+# interface intermediate 
+.true.                     2    2    68634.88   155340.91  0.0045     -0.0045
+interface.inter.dat
+
+# interface at 1km depth
+.false. 203 205 -30.0 5.1 0.05000000000000071 -0.04999999999999982
+ptopo.xyz.2.dat
+
+# interface (topography, top of the mesh)
+.false. 203 205 -30.0 5.1 0.05000000000000071 -0.04999999999999982
+ptopo.xyz.1.dat
+
+# for each layer, we give the number of spectral elements in the vertical direction
+# layer number 1 (mantle layer)
+ 5 
+# layer number 2 (crust layer)
+ 7 
+# layer number 3 (crust layer)
+ 2 
+# layer number 4 (top layer)
+ 1 

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/interfaces.dat

@@ -0,0 +1 @@
+../../topo_data/ptopo.xyz.1.dat

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/ptopo.xyz.1.dat

@@ -0,0 +1 @@
+../../topo_data/ptopo.xyz.2.dat

EXAMPLES/real_world/Moon_Apollo_impact/DATA/meshfem3D_files/ptopo.xyz.2.dat

@@ -0,0 +1,127 @@
+# Apollo 16 booster impact
+
+---
+
+![screenshot of mesh with Vs velocities](./REF_SEIS/image.mesh-vs.jpg)
+
+This is a simple example for a Moon simulation. It takes the lunar topography and sets up a model
+around the Apollo 16 booster impact and the seismic station installed by Apollo 12.
+
+The Apollo 16 S-IVB booster impact (Onodera et al. 2022; Wagner et al. 2017) details are:
+* time of impact approximately on April 19, 1972, 21:02:04 (UTC)
+* location is at lat = 1.9210, lon = 335.3770  (deg)
+* impact angle specified by angle from horizon:  79, heading angle (NE): 255.3 (deg)
+* relatively low impact speed v ~ 2.6 km/s, and assumed booster mass (stage at impact) m ~ 13,973 kg
+* momentum p = m * v ~ 3.6 * 10^7 [kg m/s], with an estimated deceleration time dt ~ 0.1 s,
+  leads to an average force F = p / dt ~ 3.6 * 10^8 [N]
+
+
+Moon stations:
+  ```
+  S12        XA        -3.01084     -23.42456       0.0    0.0      # Apollo 12
+  S14        XA        -3.64450     -17.47753       0.0    0.0
+  S15        XA        26.13407       3.62981       0.0    0.0
+  S16        XA        -8.97577      15.49649       0.0    0.0
+  ```
+
+For meshing, we use a Moho interface at an estimated 30.8 km depth (at Apollo 12 landing site) and a 3-layer model.
+The top, middle and bottom layers have regolith, crustal and mantle rock properties respectively.
+The values are taken from Garcia et al. (2011) VPREMOON model.
+
+----
+
+## Setup
+
+The setup of the model is meant for a low-resolution movie simulation that runs on 8 CPU-cores.
+The mesh is build for the region around the booster impact and the Apollo 12 seismic station,
+an area from longitude -30 to -20 degrees, latitude from -5 to 5 degrees.
+
+The setup is done by the following steps:
+
+* **Step 1:** Mesh setup
+
+  We will use the in-house mesher `xmeshfem3D` to create our simulation mesh. The main files in `DATA/meshfem3D_files/Mesh_Par_file`
+  and `interfaces.dat` are setup such that we will have 3 layers, where the Moho surface at 30.8 km depth is added manually by
+  creating a flat interface file `interface.moho.dat`.
+
+  To facilitate the meshing with doubling-layers, we also add an additional intermediate flat interface at 8 km depth from where the mesh
+  will start the stretching to account for the topography. This avoids issues with elements having negative Jacobians due to the stretching
+  and doubling of the mesh element sizes.
+
+  In `interfaces.dat`, we also added two more entries for the surfaces `ptopo.xyz.*.dat` that will be created and modified accordingly in the next step.
+
+
+* **Step 2:** Topography
+
+  To create the mesh topographic surface of the Moon, we do:
+  ```
+  > cd EXAMPLES/real_world/Moon_Apollo_impact
+  > ln -s ../../../utils/scripts/run_get_simulation_topography.py
+  ```
+  and then
+  ```
+  > ./run_get_simulation_topography.py -30 -5 -20 5 --SRTM=moon --toposhift=1000.0 --toposcale=0.9
+  ```
+
+  This script creates the files needed in a subfolder `topo_data/`, in particular the topography surface (`ptopo.xyz.1.dat`) and
+  a down-shifted surface (`ptopo.xyz.2.dat`, by approximately 1 km and a slight down-scaling of the topography by factor `0.9`)
+  to defined our regolith elements within this layer.
+
+  The script also modified the entries in the `interfaces.dat` to read the correct number of entries and increments.
+  It also modifies `Mesh_Par_file` to set the correct mesh lat/lon-dimensions and lunar projection zone number (as the `UTM_PROJECTION_ZONE`).
+
+  > [!NOTE]
+  > For Moon simulations, instead of the Universal Transverse Mercator (UTM) projection valid for Earth, the mesher will use a Lunar Transverse Mercator (LTM) projection (for latitudes between [-82,82] degrees) or a Lunar Polar Stereographic (LPS) at North/South pole regions. To enable these lunar projections, use a model name starting with `moon_***` (see `moon_default` in `Par_file`; or use `moon_tomo` for combining it with a tomography model).
+  > The zone numbers are positive for the Northern hemisphere and negative for Southern hemisphere. LTM uses zone numbers in the range +/- [1,45], LPS uses 46 for North pole and -46 for South pole regions. The script `run_get_simulation_topography.py` will output the corresponding zone number for the specified area (using its midpoint to determine the zone).
+
+
+
+* **Step 3:** source and station
+
+  The station "S12" added in `DATA/STATIONS` is the Apollo 12 seismic station at location lat/lon = (-3.01084, -23.42456) (deg).
+  For the Apollo 16 S-IVB booster impact, we use a point force `DATA/FORCESOLUTION` at location lat/lon = (1.9210, 335.3770) (deg).
+  The source time function is a Gaussian with a force factor and direction estimated by the informations given above.
+
+
+For convenience, this example folder has all these required setup files to run directly the wave simulation.
+
+---
+
+## Scattering perturbations
+
+For this moon model, we will add random 3D scattering perturbations with a von Karman distribution in the regolith and crustal layers. The following entries can be added in the `DATA/Par_file` to add such perturbations:
+```
+## scattering
+# adds scattering perturbations to velocity model
+SCATTERING_PERTURBATIONS        = .true.
+# perturbation strength
+SCATTERING_STRENGTH             = 0.4d0
+# correlation factor
+SCATTERING_CORRELATION          = 5.d0
+# list of material ids to apply scattering perturbations (0 == for all materials)
+SCATTERING_MATERIAL_IDS         = 2,3
+```
+Here we chose some arbitrary values just to visualize its effect. The settings will create perturbations with a strength of 40% and a von Karman correlation factor 5.0 that leads to an approximated correlation length of ~664 m (see output in `OUTPUT_FILES/output_generate_databases.txt`). The material IDs "2,3" are depicting the crustal and regolith material layers defined in the mesh by `Mesh_Par_file`.
+
+
+----
+
+## Wave simulation:
+
+  To run our simulation, just type:
+  ```
+  > ./run_this_example.sh
+  ```
+
+  The simulation is meant as a movie simulation. To suppress numerical noise for a cleaner movie visualization, we enlarged the source half-duration
+  by `HDUR_MOVIE = 8.0` in the `Par_file`.
+
+  The solver will create binary movie files in `OUTPUT_FILES/` folder. To convert them to VTK-files that can be visualize e.g. with Paraview,
+  we run the script:
+  ```
+  > ./xcreate_movie_files.sh
+  ```
+
+
+
+