# SEM cartesian workflow

This note summarizes the general workflow to run SEM cartesian. To start, I recommend consulting the manual as it's comprehensively written. This note is intended to highlight the most crucial aspects to consider. And there are some notes on creating a cartesian mesh (with topography or tomography) using SPECFEM's internal mesher.

## 0. Download and install
Download from github repo and compile

git clone https://github.com/SPECFEM/specfem3d.git
cd specfem3d
./configure (./configure --with-cuda=10  if you have NVIDIA CUDA available)
make

## 1. Check DATA/Par_file

- a. if you want to work with geographical WGS84 coordinates for your CMTSOLUTION and STATIONS, set `SUPPRESS_UTM_PROJECTION` to `.false.` and set the corresponding `UTM_PROJECTION_ZONE` you want to set (otherwise if working in cartesian coordinates set `SUPPRESS_UTM_PROJECTION` to `.true.`)
- b. set `NPROC` and `NSTEP`, `DT` (can be revisited later)
- c. if you want to consider your own 3D tomography model, set `MODEL = tomo`
- d. if you have NVIDIA GPU with cuda, turn `GPU_MODE = .true.` so that you can benefit from it.

## 2. Define the mesh

### Check DATA/meshfem3D_files/Mesh_Par_file

- a. check also `SUPPRESS_UTM_PROJECTION` to be the same as how it set in `DATA/Par_file`
- b. set minimum and maximum of longitude and latitude, as well as `DEPTH_BLOCK_KM`
- c. `NEX_XI` and `NEX_ETA` define the number of elements along surface edges. For number of elements along the depth direction set this in the last line of `DATA/meshfem3D_files/interfaces.dat`
- d. `NPROC_XI*NPROC_ETA` would be the MPI processes to be used.
- e. (optional, creating simple model)
- If you want to create a simple model like a multiple homogeneous layers with some rectangular heterogeneity blocks. Check the ""Domain materials"" part and for different domain assign the corresponding value. Then define the horizontal and depth extension in Domain regions.
- For the simplest case of a homogeneous block, just set `1 NEX_XI 1 NEX_ETA 1 NEX_NZ 1`
- f. (optional, including topography)
- check `DATA/meshfem3D_files/interfaces.dat`. Its concept is that you assign a reference depth, and make perturbation in vertical direction at that depth with information provided by another file. So if we want include topography we want to set the reference depth to be at the surface and we provide a file in following format to describe the topography.

Now we focus only the surface of the mesh and suppose it is split into 3x2 elements, it would look like (assuming east for right and north for upward)

    x x x
    x x x

if I want to assign a topography for this I provide a file like this in one column

    -2000
    -3000
    -4000
    1000
    2000
    4000

it's defined along XI first and along ETA starting from the **south-west corner**. The direction is defined as positive upward, so a positive value represents bathymetry above sea level, and a negative value represents depths below sea level. The arrangement would look like this:

    1000  2000  4000
    -2000 -3000 -4000

Note that the unit for topography is meter. Make sure that the loop in your code would first loop over XI (or x) direction and then ETA (or y). Something like
do y=1,ny
do x=1,nx
write(*,*)topo(i,j)
enddo
enddo

Now if the file is created and its name is topo.dat. We back to interfaces.dat to specify the information of this file. For example, to include topography of the above case

    .false. 3 2 119.00 21.00 0.004 0.004
    topo.dat

The meaning of each argument is

    SUPPRESS_UTM_PROJECTION  NXI  NETA LONG_MIN   LAT_MIN    SPACING_XI SPACING_ETA 

Then at the end of interfaces.dat set the value to the number of grids along depth direction. for example, set to 4 if you want to have 4 elements along the depth direction so the entire mesh will have 3\*2*4 elements.

- z. if it's done then run (change `NPROC` to number of processors )

            mpirun -np ${NPROC} ./bin/xmeshfem3D 

            ./bin/xdecompose_mesh ${NPROC} ./MESH ./DATABASES_MPI

## 3. Assign the model parameters

1. If there is no tomography to be included, run directly:

        mpirun -np ${NPROC} ./bin/xgenerate_databases

2. (optional, including tomography)
- If you want to consider 3D tomography, you have to create a folder `DATA/tomo_files` and the name of the input will be `tomography_model.xyz`.
- I recommend you to have a look at the source code `src/generate_databases/model_tomography.f90`. It's actually self-explained (see subroutine `read_model_tomography()`). Following the code, the 3D tomography for elastic simulation will look like:

The first four lines define the header, then the tomography value from the first point to the last point:
1. `x0 y0 z0 xE yE zE`
2. `dx dy dz`
3. `nx ny nz`
4. `vp_min vp_max vs_min vs_max density_min density_max`
5. `x0 y0 z0 Vp Vs Rh Qp Qs`
6. `x1 y1 z1 Vp Vs Rh Qp Qs`
...

The first line defines the south-west-bottom corner of the box (`x0,y0,z0`) and the north-east-top corner of the box (`xE,yE,zE`). Note that we have to use UTM coordinates here. The value shall increment first along the XI (x) direction, then the ETA (y) direction, and finally the Z direction.

An example would look like this:

      -19783.2232 2341691.2263 -300000.0000 601788.6581 2931280.9724 6000
      3000 3000 3000
      201 201 103
      2000. 9000. 1000. 5000. 1100 3000
      -19783.2232 2341691.2263 -300000.0000    8189.3000    4786.4000    3515.0875  800 560
      -16647.5228 2341592.7422 -300000.0000    8189.1000    4786.4000    3515.0125  800 560
      -13511.9291 2341494.8593 -300000.0000    8188.9000    4786.4000    3514.9375  800 560
      -10376.4413 2341397.5775 -300000.0000    8189.1000    4786.4000    3515.0125  800 560
      -7241.0589 2341300.8966 -300000.0000    8189.3000    4786.4000    3515.0875  800 560
      -4105.7811 2341204.8165 -300000.0000    8189.8000    4786.4000    3515.2750  800 560
      -970.6073 2341109.3371 -300000.0000    8190.4000    4786.4000    3515.5000  800 560
      2164.4630 2341014.4583 -300000.0000    8190.3000    4786.4000    3515.4625  800 560
      5299.4307 2340920.1800 -300000.0000    8190.2000    4786.4000    3515.4250  800 560
      8434.2963 2340826.5020 -300000.0000    8189.3000    4786.4000    3515.0875  800 560

Note that the unit for velocity is m/s and for density is kg/m^3.

## 4. STATION and CMTSOLUTION
Check the specfem cartesian manual it's already well documented.

## 5. Run the simulation
Now the mesh and its property is defined, we can launch simulation with

    mpirun -np ${NPROC} ./bin/xspecfem3D

And the output files will be in ``OUTPUT_FILES``