diff --git a/.DS_Store b/.DS_Store index 08afbf7..803c053 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/coupling.py b/coupling.py index d5905c4..2c6c52e 100644 --- a/coupling.py +++ b/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() diff --git a/coupling_specfem.py b/coupling_specfem.py index 8b94f5b..b39f92b 100644 --- a/coupling_specfem.py +++ b/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 diff --git a/submit_prepare_field.sh b/submit_prepare_field.sh deleted file mode 100644 index 7f7a014..0000000 --- a/submit_prepare_field.sh +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash -#SBATCH --nodes=3 -#SBATCH --ntasks-per-node=40 -#SBATCH --time=02:00:00 -#SBATCH --array=4-4 -#SBATCH --account=rrg-liuqy -#SBATCH --mem=0 -#SBATCH --partition=compute - -#module load gcc openmpi python/3.9.8 parallel-hdf5/gcc-8.3.0 -basedir=/home/l/liuqy/nqdu/scratch/axisem/SOLVER/ak135.$SLURM_ARRAY_TASK_ID/ -outdir=output.$SLURM_ARRAY_TASK_ID -coordir=/home/l/liuqy/nqdu/scratch/NEChina_mesh/DATABASES_MPI/ - -#basedir=/home/l/liuqy/nqdu/scratch/axisem/SOLVER/ak135_th -#outdir=output.th -#coordir=/scratch/l/liuqy/liutia97/cube2sph_injection_test/SPECFEM3D-with-Cube2sph-and-PML/utils/cube2sph/cube2sph_examples/northeast_china/DATABASES_MPI/ -mpirun -np 120 python coupling.py $basedir $coordir P 250 100 $outdir diff --git a/submit_transpose.sh b/submit_transpose.sh index 7747f0d..f6566f8 100644 --- a/submit_transpose.sh +++ b/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 \ No newline at end of file +mpirun -np 8 python ./transpose_fields.py $axisem_data_dir MZZ MXZ_MYZ MXY_MXX_M_MYY MXX_P_MYY \ No newline at end of file diff --git a/syn_seismogram.py b/syn_seismogram.py index 00d23d6..0ebd487 100644 --- a/syn_seismogram.py +++ b/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) diff --git a/test/specfem3d_cart/README.md b/test/specfem3d_cart/README.md index 10331cd..c4cf5fb 100644 --- a/test/specfem3d_cart/README.md +++ b/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`. \ No newline at end of file