new format

coupling.py

@@ -6,6 +6,90 @@
 from scipy.interpolate import interp1d 
 import sys 
 from utils import cart2sph,diff1,allocate_task
+from scipy.io import FortranFile
+
+def get_wavefield_proc(args):
+    iproc,basedir,coordir,outdir,tvec = args
+    datadir = coordir
+    file_trac = datadir + "proc%06d_wavefield_discontinuity_faces"%iproc
+    file_disp = datadir + "proc%06d_wavefield_discontinuity_points"%iproc
+    outbin = "%s/proc%06d_wavefield_discontinuity.bin"%(outdir,iproc)
+
+     # read database
+    db = AxiBasicDB()
+    db.read_basic(basedir + "/MZZ/Data/axisem_output.nc4")
+    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
+    nt1 = len(t1)
+    dt1 = t1[1] - t1[0]
+
+    # create datafile for displ/accel
+    data = np.loadtxt(file_disp,ndmin=2)
+    r,stla,stlo = cart2sph(data[:,0],data[:,1],data[:,2])
+    stel = -6371000 + r
+    f = h5py.File("%s/displ_proc%06d.h5"%(outdir,iproc),"w")
+    dset_d = f.create_dataset("displ",(nt1,len(r),3),dtype=np.float32,chunks=True)
+    dset_a = f.create_dataset("accel",(nt1,len(r),3),dtype=np.float32,chunks=True)
+
+    # compute displ/accel on the injection boundaries
+    if iproc == 0: print("synthetic displ/accel ...")
+    for ir in range(len(stla)):
+        #print(f"synthetic displ/accel for point {ir+1} in proc {iproc} ...")
+        ux1,uy1,uz1 = db.syn_seismo(stla[ir],stlo[ir],stel[ir],'xyz',basedir + '/CMTSOLUTION')
+
+        ux = interp1d(t0,ux1,bounds_error=False,fill_value=0.)(t1)
+        uy = interp1d(t0,uy1,bounds_error=False,fill_value=0.)(t1)
+        uz = interp1d(t0,uz1,bounds_error=False,fill_value=0.)(t1)
+        dset_d[:,ir,0] = np.float32(ux)
+        dset_d[:,ir,1] = np.float32(uy)
+        dset_d[:,ir,2] = np.float32(uz)
+
+        dset_a[:,ir,0] = diff1(diff1(ux,dt1),dt1)
+        dset_a[:,ir,1] = diff1(diff1(uy,dt1),dt1)
+        dset_a[:,ir,2] = diff1(diff1(uz,dt1),dt1)
+    f.close()
+
+    # compute traction on the injection boundaries
+    data = np.loadtxt(file_trac,ndmin=2)
+    r,stla,stlo = cart2sph(data[:,0],data[:,1],data[:,2])
+    stel = -6371000 + r
+    f1 = h5py.File("%s/traction_proc%06d.h5"%(outdir,iproc),"w")
+    dset_t = f1.create_dataset("trac",(nt1,len(r),3),dtype=np.float32,chunks=True)
+    if iproc == 0: print("synthetic traction ...")
+
+    for ir in range(len(stla)):
+        #print(f"synthetic traction for point {ir+1} in proc {iproc} ...")
+        sig_xyz = db.syn_stress(stla[ir],stlo[ir],stel[ir],basedir + '/CMTSOLUTION')
+        Tx = np.zeros((db.nt)); Ty = Tx *  1.; Tz = Tx * 1. 
+
+        nx = data[ir,3]; ny = data[ir,4]; nz = data[ir,5]
+        Tx = sig_xyz[0,:] * nx + sig_xyz[5,:] * ny + sig_xyz[4,:] * nz 
+        Ty = sig_xyz[5,:] * nx + sig_xyz[1,:] * ny + sig_xyz[3,:] * nz 
+        Tz = sig_xyz[4,:] * nx + sig_xyz[3,:] * ny + sig_xyz[2,:] * nz 
+
+        dset_t[:,ir,0] = interp1d(t0,Tx,bounds_error=False,fill_value=0.)(t1)
+        dset_t[:,ir,1] = interp1d(t0,Ty,bounds_error=False,fill_value=0.)(t1)
+        dset_t[:,ir,2] = interp1d(t0,Tz,bounds_error=False,fill_value=0.)(t1)
+
+    f1.close()
+
+    # write final binary for specfem_injection
+    f = h5py.File("%s/displ_proc%06d.h5"%(outdir,iproc),"r")
+    f1 = h5py.File("%s/traction_proc%06d.h5"%(outdir,iproc),"r")
+    fileio = FortranFile(outbin,"w")
+    for it in range(nt1):
+        fileio.write_record(f['displ'][it,:,:])
+        fileio.write_record(f['accel'][it,:,:])
+        fileio.write_record(f1['trac'][it,:,:])
+    fileio.close()
+    f.close()
+    f1.close()
 
 def get_trac_proc(args):
     iproc,basedir,coordir,outdir,tvec = args
@@ -32,7 +116,7 @@ def get_trac_proc(args):
     if tvec is None:
         t1 = np.arange(20000) * 0.05 - db.shift 
     else:
-        t1 = tvec
+        t1 = tvec - db.shift
     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)
@@ -88,7 +172,7 @@ def get_displ_proc(args):
     if tvec is None:
         t1 = np.arange(20000) * 0.05 - db.shift 
     else:
-        t1 = tvec
+        t1 = tvec - db.shift
     nt1 = len(t1)
     dt1 = t1[1] - t1[0]
     dset = f.create_dataset("field",(nt1,len(r),6),dtype=np.float32,chunks=True)
@@ -147,8 +231,9 @@ def main():
     #basedir = '/home/l/liuqy/nqdu/scratch/axisem/SOLVER/ak135'
     for i in range(startid,endid+1):
         args = (i,basedir,coordir,outdir,t1)
-        get_trac_proc(args)
-        get_displ_proc(args)
+        get_wavefield_proc(args)
+        #get_trac_proc(args)
+        #get_displ_proc(args)
 
     MPI.Finalize()
 

coupling_specfem.py

@@ -80,12 +80,12 @@ def get_field_proc(args):
 
     if iproc == 0: print("synthetic traction/velocity ...")
     for ir in range(npts):
-        print(f"synthetic traction for point {ir+1} of {npts} in proc {iproc} ...")
+        #print(f"synthetic traction for point {ir+1} of {npts} in proc {iproc} ...")
 
         # get rotation matrix from (xyz) to (enz)
         R = rotation_matrix(np.deg2rad(90-stla[ir]),np.deg2rad(stlo[ir]))
         tmp = R[:,1] * 1.
-        R[:,1] = -R[:,0] * 1. # \hat{n}_n is -\hat{\theta}
+        R[:,1] = -R[:,0] * 1. # \hat{e}_n is -\hat{\theta}
         R[:,0] = tmp * 1.
         R = R.T
 

submit_transpose.sh

@@ -10,4 +10,4 @@ module load gcc openmpi python/3.9.8 parallel-hdf5/gcc-8.3.0
 
 # solver dir
 axisem_data_dir=.. # like /path/to/axisem/SOLVER/ak135
-mpirun -np 8 python ../transpose_fields.py $axisem_data_dir MZZ MXZ_MYZ MXY_MXX_M_MYY MXX_P_MYY
+mpirun -np 8 python ./transpose_fields.py $axisem_data_dir MZZ MXZ_MYZ MXY_MXX_M_MYY MXX_P_MYY

syn_seismogram.py

@@ -1,21 +1,25 @@
 #from database1 import AxiBasicDB as AxiBasicDB1
 from bak.database import AxiBasicDB
 import numpy as np 
-from obspy import read_events,Trace
+from obspy import read_events
 import matplotlib.pyplot as plt 
+import os 
+
 basedir = '../SOLVER/ak135_th/'
-#basedir = '../SOLVER/ak135_th.old/'
+
+# create dir
+os.makedirs("SEISMOGRAMS",exist_ok=True)
 
 db:AxiBasicDB = AxiBasicDB()
 db.read_basic(basedir + "/MZZ/Data/axisem_output.nc4")
 db.set_iodata(basedir)
 
-cmtfile = "CMTSOLUTION"
-cat = read_events(cmtfile)[0]
+cmtfile = "/CMTSOLUTION"
+cat = read_events(basedir + cmtfile)[0]
 org = cat.origins[0]
 starttime = org.time
 xs = np.array([org.longitude,org.latitude,org.depth])
-mzz,mxx,myy,mxz,myz,mxy = db.read_cmt(cmtfile)
+mzz,mxx,myy,mxz,myz,mxy = db.read_cmt(basedir + cmtfile)
 mij = np.array([mxx,myy,mzz,2 * myz,2 * mxz,2 * mxy])
 
 stacords = np.loadtxt(basedir + "/MZZ/STATIONS",dtype = str)
@@ -27,13 +31,13 @@
     stla,stlo = np.float32(stacords[i,2:4])
     ue,un,uz = db.syn_seismo(stla,stlo,0.,'enz',basedir + cmtfile)
     name = stacords[i,0] + "_" + stacords[i,1]
-    newname = "SEISMOGRAMS/" + name + "_disp_post_mij_conv0000"
+    newname = "SEISMOGRAMS/" + name
     dataout = np.zeros((len(un),2)) 
     dataout = np.zeros((len(un),2)) 
     dataout[:,0] = t 
     dataout[:,1] = ue 
-    np.savetxt(f'{newname}_E.dat',dataout)
+    np.savetxt(f'{newname}.BXE.dat',dataout)
     dataout[:,1] = un
-    np.savetxt(f'{newname}_N.dat',dataout)
+    np.savetxt(f'{newname}.BXN.dat',dataout)
     dataout[:,1] = uz
-    np.savetxt(f'{newname}_Z.dat',dataout)
+    np.savetxt(f'{newname}.BXZ.dat',dataout)

test/specfem3d_cart/README.md

@@ -1,5 +1,5 @@
 # Prepare Boundary Points
-your should provide files `proc*_normal.txt`. This file can be obtained from `DATABASES_MPI` in **SPECFEM3D**'s `Par_file`. 
+your should provide files `proc*_normal.txt`. This file can be obtained from `LOCAL_PATH` in **SPECFEM3D**'s `Par_file`. Your can generate all these files by set `COUPLE_WITH_INJECTION_TECHNIQUE = .TRUE.` and `INJECTION_TYPE = 3` before you generate your mesh.  
 
 # OUTPUT
 The program will print `proc*_axisem_sol` files that contain velocity/tractions at each time step. You can set parameters in `create_interfaces.sh`.