Merge pull request #1 from nqdu/devel

version 2.0 pull request

README.md

@@ -8,9 +8,7 @@
 Part of the code are adapted from [Instaseis](https://github.com/krischer/instaseis), so LGPL license is applied.
 
 
-## This is the test version, user manual and examples will be updated later ...
-
-## Download required packages
+## 1. Download required packages
 1. **Compilers:** c++/Fortran Compilers which support c++14 (tested on `GCC >=7.5`, `ICC >=19.2.0`), `cmake >= 3.12`, and MPI libraries.
 
 2. create a new environment with conda:
@@ -20,7 +18,7 @@ conda activate axisem_lib
 conda install numpy scipy numba pyproj
 pip install pybind11-global
 ```
-3. Install [parallel-hdf5](https://support.hdfgroup.org/HDF5/PHDF5/), [mpi4py](https://mpi4py.readthedocs.io/en/stable/install.html) and [h5py-mpi](https://docs.h5py.org/en/stable/mpi.html), and .
+3. Install [parallel-hdf5](https://support.hdfgroup.org/HDF5/PHDF5/), [mpi4py](https://mpi4py.readthedocs.io/en/stable/install.html) and [h5py-mpi](https://docs.h5py.org/en/stable/mpi.html), and [netcdf-fortran](https://docs.unidata.ucar.edu/netcdf-fortran/current/).
 
 4. build `AxiSEMLib` by using:
 ```bash
@@ -29,18 +27,39 @@ cmake .. -DCXX=g++ -DFC=gfortran  -DPYTHON_EXECUTABLE=`which python`
 make -j4; make install 
 ```
 
-## Download AxiSEM-1.4
-Download version 1.4 of [axisem](https://github.com/geodynamics/axisem). And install all required libraries from it's manual.
+## 2. Prepare AxiSEM Mesh
+- Go to directory `axisem` and change compiler options in `make_axisem.macros`. Remember to set `USE_NETCDF = true`, set `NETCDF_PATH`.
+
+- Go to `MESHER/`, set parameters including `DOMINANT_PERIOD` and number of slices in `inparam_mesh`. If you want to use a smoothed version of ak135 model, you can run scripts under `smooth_model/ak135.py`, and set the parameters like:
+```bash
+BACKGROUND_MODEL external
+EXT_MODEL ak135.smooth.bm
+```
+
+- Run mesh generation `./submit.csh` and `./movemesh.csh mesh_name`. Then the mesh files will be moved to `SOLVER/MESHES` as the `mesh_name` you set.
+
+## Prepare AxiSEM Solver files
+There are two files you should edit: `inparam_basic` and `inparam_advanced`.
 
-Ssubstitute `nc_routines.F90` in `axisem/SOLVER/` by the same file in `axisem_files.tar.gz`.
+In `param_basic` you should set `SEISMOGRAM_LENGTH` as you required, and you should set `ATTENUATION` to `false` because the current version only support isotropic model. And you can set some other parameters like `SIMULATION_TYPE`. If `SIMULATION_TYPE` is not `moment`, you should also edit `inparam_source`.
 
-## Prepare AxiSEM Files
-In `SOLVER/inparam_advanced`, you should set the several parameters:
+In `inparam_advanced`, you should set the part of several parameters as below:
 ```bash 
+# GLL points to save, starting and ending GLL point index 
+# (overwritten with 0 and npol for dumptype displ_only)
+KERNEL_IBEG         0
+KERNEL_IEND         4
+KERNEL_JBEG         0
+KERNEL_JEND         4
+
 KERNEL_WAVEFIELDS   true
 KERNEL_DUMPTYPE     displ_only
 KERNEL_SPP          8/16/32 (depend on your dominant frequency)
 
+# you should add this one 
+# KERNEL  dump  after KERNEL_T0
+KERNEL_T0       200. 
+
 # epicenter distance
 KERNEL_COLAT_MIN   25.
 KERNEL_COLAT_MAX   100.
@@ -50,11 +69,10 @@ KERNEL_COLAT_MAX   100.
 KERNEL_RMIN        5000.
 KERNEL_RMAX        6372.
 ```
-Then you can edit several files: `inparam_mesh`,`inparam_basic`, `inparam_source`,`CMTSOLUTION`,`STATIONS`.
-And if you want to use a smoothed-ak135 model, please check parameters in  `smooth_model/ak135_smooth.py`.
+Then you can prepare your `CMTSOLUTION` and `STATIONS`
 
 ## Run AxiSEM simulation
-Build **AxiSEM** with `USE_NETCDF` mode, and run it on your cluster. 
+Run it on your cluster.
 
 ## Transpose the output field.
 Set the variables in `submit_transpose.sh`, then:

bak/database.py

@@ -59,6 +59,7 @@ def read_basic(self,ncfile:str) -> None:
         self.shift = fio.attrs['source shift factor in sec'][0]
         self.nt = len(fio['snapshots'])
         self.nglob = len(fio['gllpoints_all'])
+        self.t0 = fio.attrs['dump_t0']
 
         # read source parameters
         self.evdp = fio.attrs['source depth in km'][0] * 1.

coupling.py

@@ -21,11 +21,8 @@ def get_wavefield_proc(args):
     db.set_iodata(basedir)
 
     # time vector
-    t0 = np.arange(db.nt) * db.dtsamp - db.shift
-    if tvec is None:
-        t1 = np.arange(20000) * 0.05 - db.shift 
-    else:
-        t1 = tvec - db.shift
+    t0 = np.arange(db.nt) * db.dtsamp + db.t0
+    t1 = tvec.copy()
     nt1 = len(t1)
     dt1 = t1[1] - t1[0]
 
@@ -124,11 +121,8 @@ def get_trac_proc(args):
     stel = -6371000 + r
 
     # create dataset
-    t0 = np.arange(db.nt) * db.dtsamp - db.shift
-    if tvec is None:
-        t1 = np.arange(20000) * 0.05 - db.shift 
-    else:
-        t1 = tvec - db.shift
+    t0 = np.arange(db.nt) * db.dtsamp + db.t0 
+    t1 = tvec.copy()
     nt1 = len(t1)
     dset = f.create_dataset("field",(nt1,len(r),3),dtype=np.float32,chunks=True)
     field = np.zeros((nt1,3),dtype=np.float32)
@@ -180,11 +174,8 @@ def get_displ_proc(args):
     stel = -6371000 + r
 
     # create dataset
-    t0 = np.arange(db.nt) * db.dtsamp - db.shift
-    if tvec is None:
-        t1 = np.arange(20000) * 0.05 - db.shift 
-    else:
-        t1 = tvec - db.shift
+    t0 = np.arange(db.nt) * db.dtsamp + db.t0 
+    t1 = tvec.copy()
     nt1 = len(t1)
     dt1 = t1[1] - t1[0]
     dset = f.create_dataset("field",(nt1,len(r),6),dtype=np.float32,chunks=True)

coupling_specfem.py

@@ -6,15 +6,7 @@
 from scipy.interpolate import interp1d 
 import sys 
 from pyproj import Proj
-
-def diff1(u,dt):
-    u1 = u * 0. 
-    nt = len(u)
-    u1[1:nt-1] = (u[2:] - u[0:nt-2]) / (2 * dt)
-    u1[0] = (u[1] - u[0]) / dt 
-    u1[-1] = (u[nt-1] - u[nt-2]) / dt 
-
-    return u1
+from utils import diff1
 
 
 def read_boundary_points(coordir:str,iproc:int):
@@ -39,7 +31,7 @@ def read_boundary_points(coordir:str,iproc:int):
 
 def get_field_proc(args):
     # unpack input paramters
-    iproc,basedir,coordir,outdir,tvec = args
+    iproc,basedir,coordir,outdir,tvec,UTM_ZONE = args
 
     # read database
     db = AxiBasicDB()
@@ -51,11 +43,8 @@ def get_field_proc(args):
     npts = len(xx)
 
     # create dataset
-    t0 = np.arange(db.nt) * db.dtsamp - db.shift
-    if tvec is None:
-        t1 = np.arange(20000) * 0.05 - db.shift 
-    else:
-        t1 = tvec - db.shift
+    t0 = np.arange(db.nt) * db.dtsamp + db.t0
+    t1 = tvec.copy()
     dt1 = t1[1] - t1[0]
     nt1 = len(t1)
 
@@ -73,7 +62,7 @@ def get_field_proc(args):
         return 0
 
     # convert to spherical coordinates
-    p = Proj(proj='utm',zone=10,ellps='WGS84')
+    p = Proj(proj='utm',zone=UTM_ZONE,ellps='WGS84')
     stlo,stla = p(xx,yy,inverse=True)
     r = zz + 6371000
     stel = -6371000 + r
@@ -124,14 +113,15 @@ def get_field_proc(args):
 
 def main():
     from utils import allocate_task
-    if len(sys.argv) !=7:
-        print("Usage: ./this h5dir coor_dir t0 dt nt outdir")
+    if len(sys.argv) !=8:
+        print("Usage: ./this h5dir coor_dir UTM_ZONE t0 dt nt outdir")
         exit(1)
     basedir = sys.argv[1]
     coordir = sys.argv[2]
-    t0,dt = map(lambda x: float(x),sys.argv[3:5])
-    nt = int(sys.argv[5])
-    outdir = sys.argv[6]
+    UTM_ZONE = int(sys.argv[3])
+    t0,dt = map(lambda x: float(x),sys.argv[4:6])
+    nt = int(sys.argv[6])
+    outdir = sys.argv[7]
     os.system(f'mkdir -p {outdir}')
 
     # mpi 
@@ -155,10 +145,8 @@ def main():
 
     # time window
     t1 = np.arange(nt) * dt + t0
-
-    #basedir = '/home/l/liuqy/nqdu/scratch/axisem/SOLVER/ak135'
     for i in range(startid,endid+1):
-        args = (i,basedir,coordir,outdir,t1)
+        args = (i,basedir,coordir,outdir,t1,UTM_ZONE)
         get_field_proc(args)
 
     MPI.Finalize()

database.py

@@ -24,6 +24,7 @@ def read_basic(self,ncfile:str) -> None:
         self.shift = fio.attrs['source shift factor in sec'][0]
         self.nt = len(fio['snapshots'])
         self.nglob = len(fio['gllpoints_all'])
+        self.t0 = fio.attrs['dump_t0']
 
         # read source parameters
         self.evdp = fio.attrs['source depth in km'][0] * 1.
@@ -92,6 +93,7 @@ def __copy__(self):
         db.shift = self.shift
         db.nt = self.nt
         db.nglob = self.nglob
+        db.t0 = self.t0
 
         # source parameters
         db.evdp = self.evdp

reciprocity.py

@@ -108,7 +108,7 @@ def reciprocity(db:AxiBasicDB,cmtfile):
     data_z = np.loadtxt(name + "Z.dat")
     data_e = np.loadtxt(name + "E.dat")
     data_n = np.loadtxt(name + "N.dat")
-    t = np.arange(db.nt) * db.dtsamp - db.shift
+    t = np.arange(db.nt) * db.dtsamp + db.t0
 
     plt.figure(1,figsize=(14,16))
     plt.subplot(3,1,1)