Last update: 2020


  Specfem3D     1. configure: ./configure FC=gfortran CC=gcc MPIFC=mpif90 --with-mpi --with-scotch-dir=/home/yujiang/specfem3d/external_libs/scotch_5.1.12b  MPI_INC=/home/yujiang/mpi/mpich3/include   FCLIBS=" "    2. makefile problem: /usr/bin/ld: cannot find -lz        solution: sudo yum install zlib-devel  #redhad linux
Iridis4: ./configure FC=gfortran CC=gcc MPIFC=mpif90 --with-mpi --with-scotch-dir=/scratch/yx1r19/specfem3d_lucas/external_libs/scotch_5.1.12b  MPI_INC=/local/software/mpich/3.0.4/gcc-4.8.1/include  FCLIBS=" "

Iridis5: ./configure FC=gfortran CC=gcc MPIFC=mpif90 --with-mpi --with-scotch-dir=/scratch/yx1r19/specfem3d_lucas/external_libs/scotch_5.1.12b  MPI_INC=/local/software/mpich/3.2.1/gcc/include FCLIBS=" "

sdsc: cd "



 ParaView    1. installation problem: /home/yujiang/Downloads/ParaView-5.5.2-Qt5-MPI-Linux-64bit/lib/paraview: error while loading shared libraries: libGLU.so.1: cannot open shared object file: No such file or directory        solution: yum install /usr/lib*/libGLU*    to install  "mesa-libGLU.x86_64 0:9.0.0-4.el7"  and  "mesa-libGLU-devel.x86_64 0:9.0.0-4.el7"    Add software path:    export PATH=/home/yujiang/Downloads/ParaView-5.5.2-Qt5-MPI-Linux-64bit/bin:$PATH     source ~/.bashrc
  Specfem3D (in xgenerate_databases)    Error negative or null 3D Jacobian found    Error detected, aborting MPI... proc      solution: dense the element layers in depth, careful the dh of the grid.   Specfem3D (in xmeshfem3D)    the order saved for topography.dat in interfaces.dat is x, y, and z, going x direction first, then y,  and z, based on the mesh coordinates. 
  Specfem3D (in xmeshfem3D and xgenerate_databases )  1. interpolation for topography:  to compute  elevations at twice-dense grids of the element, by using x times y grids of elevations provided by the user (see interfaces.dat), where the number of x and y 'should' be dense than the twice-dense grids.  2. read external tomography files into mesh: to compute tomography values (e.g., vp,vs, rho, Qk and Qu) at the GLL points of each element by using trilinear interpolation, e.g., vp_final = interpolate_trilinear(vp1,vp2,vp3,vp4,vp5,vp6,vp7,vp8), where vp1,vp2,vp3,vp4,vp5,vp6,vp7,vp8 are the eight points of a cubic of the "tomography_model.xyz" provided by the user (external tomography files). i.e., read external tomography files into a mesh and store the tomo values in GLL points of each element of the mesh.   3. read tomography file with variable discretization intervals: the number of the "tomography_model_i.xyz" files used in subroutine read_model_tomography() (i.e., the nundefMat_ext_mesh) is read in subroutine read_partition_files, which means that the lines (rows) of tomography files defined in nummaterial_velocity_file should equal to NMATERIALS  in Mesh_Par_file when use internal mesh and when read an external tomography on to the internal mesh. In the sources, I consider the IMODEL_DEFAULT and IMODEL_TOMO separately although the tomography can be read in case of IMODEL_DEFAULT, which requires to set imaterial_def=2 (see get_mode.f90 for imaterial_def = mat_ext_mesh(2,ispec), and  read_partition_file.f90 for read mat_ext_mesh(2,ispec), and save_databases.f90 for write material_index(2,ispec) , where 2 is set there).#lucas: in Specfem3D and Multi-SEM3D, I use 'MODEL=tomo' for reading external models; this can be also done using the MODEL=default (have tested for Specfem3D but not tested yet for MultiSEM3D). # lucas: the tomo only can do all the model is from external tomo.
   updated(04/04/2019): two ways to read tomography file with variable discretization intervals:                 a)  interpolate the external tomograpy files into constant intervals. NX, NY, and NZ may be different, but constant each of them. In this case, we only need one imaterial in the mech, i.e., nundefMat_ext_mesh=1, which is read in subroutine read_partition_files that opens the "Database" files produced/resulted in the xmeshfem3D .               b)  seperate the tomography files into several tomography files and each has constant intervals. NX, NY, and NZ may be different, but constant. 
The order for external tomo file: x first and then y, z is started from the bottom and gradually upper to the model top surface. more details see, e.g., the tomography_model_01.xyz in my DATA folder.
Specfem3D_Cartesian (in xmeshfem3D)
fix the UTM_PROJECTION_ZONE in xmeshfem3D for several zones since in the source package, only one utm projection zone is adopted. This may be unsuitable when a research area covers several utm zones. Solution: set the left-bottom corner as a reference point (x,y)=(0,0), and other coordinates are relevant to this reference point. -------------------------------------% 1. reference point, to be 0 and 0 in x-y coordinate systemlon_min=-20.0167;lat_min=-3.3526;
% 2. max lon and latlon_max=-6.75;lat_max=2.6025;
% 3. range of lon and latdLon=lon_max - lon_min;dLat=lat_max - lat_min;
% 4. max lon and lat in mx_end_m=dLon*111000;y_end_m=dLat*111000;
LAT_min=0;LAT_max=(lat_max - lat_min)*111000;LON_min=0;LON_max=(lon_max - lon_min)*111000;MESH_range=cat(1,LAT_min,LAT_max,LON_min,LON_max);dlmwrite('MESH_range',MESH_range,'delimiter','\t','precision','%.4f');
% source coodinate source_lat_m=(0.13-lat_min)*111000;source_lon_m=(-17.72-lon_min)*111000;source_coor=cat(1,source_lat_m,source_lon_m);dlmwrite('source_latUp_longBot',source_coor,'delimiter','\t','precision','%.4f');
% receiver coodinate% import STATIONSreceiver_lat=STATIONS(:,3);receiver_lon=STATIONS(:,4);receiver_ele=STATIONS(:,6);receiver_0=  STATIONS(:,5);for i=1:1:23receiver_lat_m(i)=(receiver_lat(i)-lat_min)*111000;receiver_lon_m(i)=(receiver_lon(i)-lon_min)*111000;endreceiver_lat_m=receiver_lat_m';receiver_lon_m=receiver_lon_m';
receiver_par=cat(2,receiver_0,receiver_0,receiver_lat_m,receiver_lon_m,receiver_0,receiver_ele);dlmwrite('STATIONS_newCoor',receiver_par,'delimiter','\t','precision','%.4f');
% for the CMPL boundary in mesh_par_file, if SUPPRESS_UTM_PROJECTION=true, we can use "meter" directly in mesh_par_file, otherwise (SUPPRESS_UTM_PROJECTION=false), use degree (lat, lon)x=0.65*111000;y=0.45*111000;z=0.0459*111000;----------------------------------------
Specfem3D (in xmeshfem3D and xgenerate_databases )1. read tomography file with variable discretization intervals:   (a) in mesh_par_file, just changed the mat_id = -1, -2, .... until the last material.  (b) in codes, changed the tomography_model.xyz into tomography_model_i.xyz in "save_databases.F90" for each specific material with mat_id.  see here:  ===========================  ! writes out undefined materials  do i = 1,NMATERIALS    domain_id = material_properties(i,7)    mat_id = material_properties(i,8)    if (mat_id < 0) then      ! format:      ! #material_id #type-keyword #domain-name #tomo-filename #tomo_id #domain_id      ! format example tomography: -1 tomography elastic tomography_model.xyz 0 2      undef_mat_prop(:,:) = ''      ! material id      write(undef_mat_prop(1,1),*) mat_id      ! name      undef_mat_prop(2,1) = 'tomography'      select case (domain_id)      case (IDOMAIN_ACOUSTIC)        undef_mat_prop(3,1) = 'acoustic'      case (IDOMAIN_ELASTIC) ! lucas IDOMAIN_ELASTIC=2 set in constant.h        undef_mat_prop(3,1) = 'elastic'      end select      ! default name     ! debug by lucas      write(file_str,"(i2.2)") i      !print *,'lucas defined file_i: ',file_str      undef_mat_prop(4,1) = 'tomography_model_'//trim(file_str)//'.xyz' ! how to write vary file name 21-(i-1)      ! debug by lucas      !print *,'lucas defined tomo file name: ',undef_mat_prop(4,1)      ! default tomo-id (unused)      write(undef_mat_prop(5,1),*) 0      ! domain-id      write(undef_mat_prop(6,1),*) domain_id      ! debug by lucas      !print *,'lucas undef mat: ',undef_mat_prop      !writes out properties      write(IIN_database) undef_mat_prop    endif  enddo   ===============================
Specfem3D (specfem3D)//1. careful the order in case of computing the temp1_new, temp2_new, and temp3_new when used NGLLX=5 in subroutine compute_forces_acoustic_fast_Deville.//2. careful the computation of tempx1_att_new,tempx2_att_new,tempx3_att_new,                                                          tempy1_att_new,tempy2_att_new,tempy3_att_new,                                                        tempz1_att_new,tempz2_att_new,tempz3_att_new, in compute_strain_in_element() if(is_CPML(ispec)) in subroutine compute_forces_viscoelastic().

//3. careful the order of (kappa_z,d_z,alpha_z,kappa_y,d_y,alpha_y,kappa_x,d_x,alpha_x) in lijk_parameter_computation() when computed for A6, A7, A8, and A9 in case of pml_compute_memory_variables_acoustic() and pml_compute_memory_variables_elastic() in file pml_compute_memory_variables.f90.//4. careful the READ_ADJSRC_ASDF in compute_add_sources_viscoelastic(), see codes in details before use it.

code details of the specfem3D, done! 


Specfem3D (specfem3D)
In FORCESOLUTION, use degree for latorUTM and longorUTM, instead of using UTM coordinates.In STATIONS, the order saved is: station_name, network_name, stlat, stlon, stele, stbur. The stlat and stlon also indicate degree here, not the UTM coordinates.
Specfem3D (specfem3D)
 Frechet Kernels: when compute frechet kernels:1) mkdir SEM folder under OUTPUT _FILES directory 2) copy .adj into /SEM, here need to make adjoint sources using utils.3 changed codes: setup_sources_receivers.f90 and compute_array_sources.f90, e.g., /..SEM/ to SEM/ to make the .adj files can be read.

--------------------------

Specfem3D globe

Tips: use different models or change a new model needs to recompile the codes.

----------------------------------





use obspy that installs via anaconda should: in the shebang (#!) script, the python should be the one located in the obspy direction of the anaconda package, e.g., #!/home/yujiang/anaconda3/envs/obspy/bin/pythonfrom obspy.core   import   read...--------------------------------------------------------------that's why use "conda activate obspy" and then type "python" works. 

######################################################
Test1 : Specfem3D (specfem3D)
x_range:1448.6 km, y_range:781.538 km
------------------------------------------------------------------------------------                  x                                                              delta_x                                         number of element: 160                        9.053 km/element   (~2.26km grid size)                 y                                                              delta_ynumber of element:  96                         8.14 km/element     (~2km grid size)
------------------------------------------------------------------------------------                  x                                                              delta_x                                         number of element: 96                          15.089 km/element   (~3.77km grid size)                 y                                                              delta_ynumber of element:  48                         16.28 km/element    (~4.07km grid size)
------------------------------------------------------------------------------------                  x                                                              delta_x                                         number of element: 64                          22.63km/element   (~5.6km grid size)                 y                                                              delta_ynumber of element:  32                         24.42 km/element    (~6.1km grid size)
------------------------------------------------------------------------------------                 x                                                              delta_x                                         number of element: 32                         45.27km/element   (~11.32km grid size)                 y                                                              delta_ynumber of element:  16                        48.85 km/element    (~12.21km grid size)
------------------------------------------------------------------------------------
######################################################

Learn tips:


FFT spectral example:
1. matlabcodes:https://dadorran.wordpress.com/2014/02/17/plot_freq_spectrum/data: https://www.dropbox.com/s/n4kpd4pp2u3v9nf/tremor_analysis.txtfr_axis=(fs/N)*[0: N-1], where N = length(signal) and fs= the number of samples per second.
Learn-seismic-tomographyhttps://computational-seismology.github.io/learn-seismic-tomography/building-blocks/pre-processing/

obspy.signal - Signal processing routines for ObsPyhttps://docs.obspy.org/packages/obspy.signal.html

################# addressed #################1. why there are bout 0.004 seconds difference: this is because the delta or sampling rate for different station differs, some about 0.016 seconds and some 0.02 seconds.2. why the specfem2sac.sh always mention that, cannot open the files stored in directory, e.g., OUTPUT_FILES/e20160401_2km/? solution: just change the path, e.g., move the folder e20160401_2km out of the folder OUTPUT_FILES.3. check the accuracy of the mesh tomography? it is corrected based on the visiable model check.   Done! 24 May 2019.

###################### to be addressed (addressed) #########1. Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.Backtrace for this error:#0  0x7FB951430697#1  0x7FB951430CDE#2  0x7FB950D7E27F#3  0x41F32E in mxm5_single_two_arrays_at_a_time.2379 at compute_forces_acoustic.f90:992#4  0x422E59 in compute_forces_acoustic_fast_deville_ at compute_forces_acoustic.f90:505#5  0x41CE28 in compute_forces_acoustic_calling_ at compute_forces_acoustic_calling_routine.f90:94#6  0x496FB2 in iterate_time_ at iterate_time.F90:239 (discriminator 1)#7  0x402F89 in xspecfem3d at specfem3D.F90:387#8  0x7FB950D6A3D4
====================================================================================   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES=   PID 7202 RUNNING AT yjx.noc.soton.ac.uk=   EXIT CODE: 136=   CLEANING UP REMAINING PROCESSES=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES===================================================================================YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Floating point exception (signal 8)This typically refers to a problem with your application.Please see the FAQ page for debugging suggestions
one temporary solution: remove the water layer, but how the true issue is?Final solution: use compute_forces_acoustic_NGLL5_fast() from specfem3D_original ( the stable version)
2. why there are still about 3 seconds difference between the observed and tomography or 1D in terms of waveforms, e.g., for the main peak.3. check how to use detrend and taper? 4. check the delta for different station?   for e20140401-142554 station        delta  I04D         0.016  I34D         0.016  L13D         0.01  L18D         0.01  L21D         0.01  L33D        0.01  L37D        0.01  L39D        0.01  S03D       0.02  S06D       0.02  S10D        0.02  S11D         0.02  S19D        0.02  S22D       0.02  S26D       0.02  S29D       0.02  S31D        0.02  S38D       0.025. check the station location that may face with phase change6. sign convention for each component, e.g., RTZ or NEZ7. Do a shallow modeling test for surface wave 

set source time function can be found in the subroutine of  get_sft_viscoelastic() of compute_add_sources_viscoelastic() in folder of specfem3D.  The hdur_gaussian is set at setup_sources_receivers.f90.
e20160829
S26D. sign convention:
method    Z        N             E1D               +         +              +tomo         -         -               +Obs            +         -(BH1)    -(BH2)
Some notes used:1. Coping data from iridis5 to local machine needs to be done at local by: scp yx1r19@iridis5_a.soton.ac.uk:/home/yx1r19/specfem3d_lucas/EXAMPLES/meshfem3D_examples/many_interfaces_PiLab_1D/OUTPUT_FILES/XS.*  ./
2. copying data from local to iridis5 also needs to be done at local by: scp -r ./SOFI3D-Release  yx1r19@iridis5_a.soton.ac.uk:/home/yx1r19/
 Harvard CMT to Nm: CMT_Nm_moment_tensor(:,:) = _Harvard_CMT_moment_tensor(:,:) * 1.d-7.  since 1 dyne.cm = 1e-7 Newton.m

Relation: Harvard CMT moment tensor (coordinate: x(S)-y(E)-z(Up)) to coordinate: X(E)-Y(N)-Z(U) ((used in Specfem3D) is in locate_source.F90.  In SOFI3D, they use x(E)-z(N)-y(Up) , where y indicate vertical direction.

Some useful debug in HPC: https://www.sphenisc.com/doku.php/software/development/debug_and_optim


Use_external_source_file:
In CMTSOLUTION, at the last line of each source, add  "./DATA/external_stf", where the file external_stf located in DATA directory and in the external_sft file, the first line is dt, following by NSTEP of the samples, where NSTEP set in Par_file. That is the first line is dt, the other lines are the stf signal. =============================================================================================================================

Obspy 

basemap:

>>> from mpl_toolkits.basemap import BasemapTraceback (most recent call last):  File "<stdin>", line 1, in <module>  File "/home/yujiang/anaconda2/envs/obspy/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py", line 156, in <module>    epsgf = open(os.path.join(pyproj_datadir,'epsg'))IOError: [Errno 2] No such file or directory: '/home/yujiang/anaconda2/envs/obspy/share/proj/epsg'
solution:conda config --add channels conda-forge conda config --set channel_priority strict conda install basemap so when $ conda list ,  the basemap will be in conda-forge channel. 

=====================================================================================================

SES3D

1. compute gradient (kernel)

[yx1r19@cyan51 TOOLS]$ ./project_kernel.exe  project visco-elastic kernels (0=no, 1=yes): 0 directory for sensitivity densities: "../DATA/OUTPUT/1.8s/"    confirm directory: ../DATA/OUTPUT/1.8s/                                                                                                                         directory for output: "../MODELS/MODELS_3D/" confirm directory: ../MODELS/MODELS_3D/  
2. plot kernel, use model.py 3. Python version to be used    solution: a) directly use the python path at the top of each python file, e.g., #!/home/yujiang/anaconda3/envs/obspy/bin/python                       b) use "#!/usr/bin/env  python" at the top of each python file, but in this way, the command should be run under the python env, e.g.,                       (obspy)[usename@pcname directory]$                 or                       (lasif)[usename@pcname directory]$                 where obspy and lasif can be set as the env via Anaconda.              if we use the direct path as a), the python file (name.py) can be run either under bash, obspy, or lasif env.  Tips: one should install the dependencies in each specific env, e.g., obspy for SES3D or lasif for FWI . Always remember to check the python path and version before using it.             ====================================================================================================================================

LASIF Installation:

1.py.test 

error: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObjectsolution:  pip install numpy==1.15.4 for my local computer
error: SyntaxError: Non-ASCII character '\xc2' in file /mainfs/scratch/yx1r19/LASIF/lasif/domain.py on line 318, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for detailssolution: -*-  coding:utf-8 -*- one the top line, but below #!/../../ 
E   DistributionNotFound: The 'mpi4py' distribution was not found and is required by the applicationsolution: pip install mpi4py

LASIF USAGE

TclError: no display name and no $DISPLAY environment variablesolution: -Y when use ssh to  login in.

Specfem2D


tomography file input: In the tomography.xyz file, in the header parameters, only the NX and NZ are used in the code, other header parameters are just for the user to read and do not actually used in the code. In this case, in the Par_file, set:#-----------------------------------------------------------------------------## velocity and density models##-----------------------------------------------------------------------------nbmodels                        = 11 -1 2000.d0 2000.d0 1176.0d0  0 0 9999.d0  9999.d0 0 0 0 0 0 0
---i.e., the 15 parameters:n, indic, val0 (rho), val1 (vp), val2(vs), val3, val4, val5(Qkappa), val6 (Qmu), val7, val8, val9, val10, val11, val12. In acoustic case, val2 should be set to zero, and the Qkappa, and Qmu do not actually used here, and we can set in the tomography.xyz file.  In each layer composed of several elements in vertical direction, the order for storing the element (indicated by ispec) starts from bottom and continue along x direction, e.g.:
6 7 8 9 101 2  3 4  5
while in each element, the order for each GLL points is:21 22 23 24 25 ...6  7  8  9  101  2  3  4  5
This means that, for preparing the tomography.xyz, the vp, vs, rho should be started from the bottom of each material layer. And the tomography only charge for that material layer (indicated by -1 ) defined in Par_file. For the whole model (may composed of several materials), we also construct the mesh from bottom and add each layer layer-by-layer until the top surface, see the interfaces.dat file. 
Codes changed for reading external tomography: (starting from line 217 to 234 of read_material_table.f90)
   else if (icodemat(i) <= 0) then ! lucas: this is used for reading external tomography, where icodemat(i)=-1, i.e., indic =-1      ! tomographic material      number_of_materials_defined_by_tomo_file = number_of_materials_defined_by_tomo_file + 1
      if (number_of_materials_defined_by_tomo_file > 1) &        call stop_the_code( &'Just one material can be defined by a tomo file for now (we would need to write a nummaterial_velocity_file)')
      ! Assign dummy values for now (they will be read by the solver). Vs must be == 0 for acoustic media anyway      !-----lucas changed below-------------      rho_s_read(i) = val0read ! -1.0d0        cp(i) = val1read !-1.0d0       cs(i) = val2read      QKappa(i) = val5read! -1.0d0      Qmu(i) = val6read  !-1.0d0      aniso3(i) = val3read ! -1.0d0          aniso4(i) = val3read ! -1.0d0    !----lucas changed above---------------

Important (Specfem2D): plane wave starting from -t0, where t0= minval( t0_source(:) ), and t0_source=1.20d0 * hdur + tshift_src(i_source) , where hdur=1/f0. That is for point and plane wave source also, both start from -1.2/f0 (=-1.2*hdur).

Specfem2D: in case of change the codes but without any direct result changed, you need to execute  "make clean" and "make all" twice.   
when use external tomo set MODEL= default, since it can be reset to MODEL= 'tomo' in setup_mesh_external_models() before calling the read_external_model() in the setup_mesh_external_models().
--------------------------------------------------3D volume visialaization: https://uk.mathworks.com/help/matlab/visualize/exploring-volumes-with-slice-planes.html

Specfem2D/3D: Stream I/O (e.g., write or read with random pos=offset) needs fortran2003.

Specfem2D: Hessian Kernels computation (the wavefield storage method, not recommended)

1. forward: Got forward wavefield s(m), s(m+dm) when set simulation_type=1 and save_forward=true, which needs to change the codes for two different names of the output wavefields.
2.adjoint: got adjoint wavefiles s*(m) by using simulation_type=3, and save_forward=true, which needs to comment the reading wavefield routine out (in iteration_time:manage_no_backward_reconstruction_io). In this step, we only need the adjoint wavefields.
3. Hessian: set simulation_type=3 and save_forward=false, which needs to read s(m) and s(m+dm) for H^(1->2), and s(m) and s*(m) for H^(2->1) , where ds*=[s*(m+dm)-s*(m)] and ds=[s(m+dm)-s(m)].
Keep in mind: 1) when save s*(m), do not forget to change the par_file for material and set the SEM/adj correct.                              2) when compute Hessian, do not forget to set SAVE_FORWARD =false, and also change the SEM/adj to perturbated adjoint source.

File name: disp_m1_iframe_at  and disp_m2_iframe_at indicate s(m) and s(m+dm).   reverse_disp_m1_iframe_at  indicates s*(m).

homogeneous model: 1 1 2900.d0 8000.d0 4800.d0 0 0 9999 9999 0 0 0 0 0 0,   x:0-800km, z:0 to -360km, 400 elements in x and 360 elements in z direction.nbregions                       = 11 400  1  360 1rectangle 10% perturbation model (perturbation model: 8km in x and 6 km in z): x:0-800km, z:0 to -360km, 400 elements in x and 250 elements in z direction.nbregions                       = 2 1 400  1  360 1201 210 251 270 2  # 20 km in x ,and 20 km in z direction.interfaces.dat   -360km to 0km 
1 1 2900.d0 8000.d0 4800.d0 0 0 9999 9999 0 0 0 0 0 02 1 2900.d0 8800.d0 5280.d0 0 0 9999 9999 0 0 0 0 0 0  #%10 perturbation

adj_source for m1  and m2 is 67 to 74 seconds

Important: in computing the Hessian kernels, the meshing grid should be rigorous the same for model m and model update (m + dm); otherwise,  the du = u(m+dm) - u(m) is not equal to zero between the model perturbation to the source.

1. 20X20 km nbregions 21  400  1   360 1201  210  251  270 2
Tips: should be tested the size of the model perturbation and the percent of the perturbation.
when NSTEP_BETWEEN_COMPUTE_KERNELS   = 1, each saved wave field has 4.6 MB, so the total save wave fields are 4.6 x 40000=179.6875 GB.
computing the H2->1 problem: Solution: use new matrix no_backward_displ_buffer_lucas_tmp to save the value of du* (see read_forward_arrays.f90 and compute_kernels.f90)
computing the H1->2 problem: Solution: use the same mesher for u(m) and u(m+dm) , also can use new matrix for du (not tested yet)
Still one unnecessary question: why we can't assign new values to displ_elastic, i.e., not the results as expected.
Notes: for the adjoint source, at least one component has values. If both the x and z componenets are zero, the code will occur error: Program received signal SIGSEGV: Segmentation fault - invalid memory reference. This is only applied to the MacOS version.
Attenuation: when considered attenuation, only the parsimonious disk storage method is recomended. The full wavefield storage method takes a large amount of disk space. Also only the P-SV case (Specfem2D) is implemented so far for wavefield simulation with attenuation, the SH wave case can be considered later.




specfem2D_FK: read external tomography file (only elastic case or acoustic case)

solution: 1. set model domain by internal default model as the same model domain of the tomography.2. run the default model (no matter how many materials) to get kernels (needs the forward and adjoint simulation)3. use the matlab package klregularmesh.m to compute the gllmodel.data (needs to prepared the regular tomography file)4. use gllmodel.data as the input model (set gllmodel=true)5. change the specfem2D.f90 code in which gllmodel is used, e.g.,  change  !    reading model perturbation with respect to the initial 1D model  # about line 2290 of the code       read(218,'(24X,3F24.8)') ratio_rho,ratio_vp,ratio_vs        valrho = valrho*(1.0+ratio_rho)       valvp  = valvp *(1.0+ratio_vp)       valvs  = valvs *(1.0+ratio_vs)
  into   !    reading new tomo model outside       read(218,'(24X,3F24.8)') ratio_rho,ratio_vp,ratio_vs  # where ratio here already indicate the new model, not ratio or perturbation       valrho = ratio_rho       valvp  = ratio_vp       valvs  = ratio_vs
This one has one problem,:when coupling with elastic with acoustic, importing the gllmodel.data do not implement in the compute_forces_acoustic.f90. I highly recommend to use the following one!!!
-----------------------------------------------------------------------------------------------

specfem2D_FK: read external tomography file (any case, problems fixed and recommend!!!)

1. in constants.h to set OUTPUT_MODEL_VELOCITY_FILE = .true. 

2. in Par_file, set assign_external_model= false, and get ''DATA/model_velocity.dat_output'' (due to 1 step), which would be in format:

3. use the matlab package klregularmesh.m to prepare "DATA/model_velocity.dat_input" (here needs to change the klregularmesh.m, i.e., 

change

FID=fopen('rhop_alpha_beta.004.00');c=fscanf(FID,'%e %e %e %e %e\n',[5 inf]);fclose(FID);d=c';xnew=d(:,1);znew=d(:,2);

into

FID=fopen('model_velocity.dat_output'). c=fscanf(FID,'%f %f %f %f %f %f\n',[6 inf]);  fclose(FID);d=c';xnew=d(:,2);znew=d(:,3);-----------------------------------# or see klregularmesh_tomo.m for details.

4. set assign_external_model = true, READ_EXTERNAL_SEP_FILE = true, to run the forward simulation.

Failed: if loading the acoustic-elastic tomo model, error: external velocity model cannot be both fluid and solid inside the same spectral element. 

Solutions: 

(1) In specfem2D.F90 (from 4056 to 4064 lines). Checking fluid/solid edge topology, i.e., the edge index of the acoustic element (e.g., 1) should fit with the edge index of the elastic element (e.g., 3). Index number, 1: bottom, 2:right, 3:top, 4:left. See the specfem2D.F90 for variables 'fluid_solid_acoustic_iedge(inum)' and 'fluid_solid_elastic_iedge(inum)'. 

(2) In read_external_model.f90, when vsext(i,j,ispec) < TINYVAL, where TINYVAL=1.d-9, the elements assumed to be acoustic usually. However, one should be careful the vsext(i,j,ispec) used for the top row elastic element that connect with the acoustic elements. Because within these top row elastic elements, their top GLL points have the vsext(i,5,ispec)=0. So the code should be set these elements are elastic like:

--------------------------------------

else if(vsext(1,1,ispec) < TINYVAL) then !lucas, should be careful and should be specfic the element, not the GLL. This is the reason why the elastic wavefield can't go into the acoustic  (see from 193 to 196 lines )

          elastic(ispec) = .false.

          poroelastic(ispec) = .false.

          any_acoustic = .true.

--------------------------------------

This is to make sure these elements (even their top-surface GLL points with vsext(i,5,ispec)=0) are elastic. Keep in mind that we need to use the same mesh as set in Par_file for external tomography model.


Sp and Ps sensitivity kernels

1. the color zone (color stripe ) of the kernels depends on the length (shape or waveform) of the adjoint source.


CTD-SEM


1. forward by coupling two solver engines, done! 2. adjoint by coupling five solvers for elastic and anelastic, done!

Seisflows on Cray #PBS

see my gmail: seisflows on PBS: issue 1: fix the setup() issue in generate_data
In Archer,  for 3D FWI only, we have muticore.py can be used (under testing), it's better to use pbs_lg (to be tested), note the run_single used or not?In irids5 (slurm), for 2D FWI, I use slurm_lg.py (from the Seisflows), not use the multicore.py because it run on login nodes.In iridis4 (PBS) for 3D FWI, we have pbs_lg (under testing) or muticore.py (done!) In iridis5 (slurm), for 3D FWI (e.g., PI-LAB 3D), I will use the slurm_lg, not use the multicore.py because it works on login nodes.
slurm_lg:-------------submit()--> sbatchrun()-->sbatchrun_single()-->sbatch--fwd()-->sbatchadj()-->sbatchxcombine_sem()-->sbatchxclip_sem()-->sbatchxsmooth_sem()-->sbatch=====================multicore-------------, depends on the system (slurm or pbs)submit()--> run()-->run_single()-->--fwd()-->adj()-->xcombine_sem()-->xclip_sem()-->xsmooth_sem()-->================pbs_lg-------------submit()-->qsub #PBS run()-->qsub #PBS run_single()-->qsub #PBS (have this or not?)--fwd()-->qsub #PBS adj()-->qsub #PBS xcombine_sem()-->qsub #PBS xclip_sem()-->qsub #PBS xsmooth_sem()-->qsub #PBS =====================

sdr2cmt.py

#!/home/yujiang/anaconda3/envs/obspy/bin/pythonfrom pyrocko import moment_tensor as mtmmagnitude = 3.11  # Magnitude of the earthquakem0 = mtm.magnitude_to_moment(magnitude)  # convert the mag to momentstrike = 296dip = 78rake = 142mt = mtm.MomentTensor(strike=strike, dip=dip, rake=rake, scalar_moment=m0)m6 = [mt.mnn, mt.mee, mt.mdd, mt.mne, mt.mnd, mt.med]  # The six MT components#print(m6/mt.scalar_moment())  # normalized MT componentsprint(m6)print('Convert to Harvard CMT solution')Mrr= mt.mddMtt= mt.mnnMpp= mt.meeMrt= mt.mndMrp=-mt.medMtp=-mt.mnem6_CMTSOLUTION = [Mrr, Mtt, Mpp, Mrt, Mrp, Mtp] # get CMTSOLUTION for Specfem3Dprint('Print to CMTSOLUTION for Specfem3D') #CMTSOLUTION file values are in dyne.cm# where in pyrocko: scalar_moment – scalar moment in [Nm], so need to time 10^7 because 1 Newton.m = 1e7 dyne.cmMrr=Mrr*1e7 Mtt=Mtt*1e7 Mpp=Mpp*1e7 Mrt=Mrt*1e7Mrp=Mrp*1e7Mtp=Mtp*1e7print('Mrr:', Mrr)print('Mtt:', Mtt)print('Mpp:', Mpp)print('Mrt:', Mrt)print('Mrp:', Mrp)print('Mtp:', Mtp)#3.11   296     78    142 

Install specfem2D on Iridis5 ./configure FC=ifort CC=icc MPIFC=mpif90 --with-mpi --with-scotch-dir=/scratch/yx1r19/specfem2d-master/external_libs/scotch_5.1.12b  MPI_INC=/local/software/mpich/3.2.1/intel/include FCLIBS=" "
./configure FC=gfortran CC=gcc MPIFC=mpif90 --with-mpi --with-scotch-dir=/scratch/yx1r19/specfem2d_develop_lucas/external_libs/scotch_5.1.12b  MPI_INC=/local/software/mpich/3.2.1/gcc/include FCLIBS=" "

specfem3D: Qp and Qs (=Qmu) from external model, and Qkappa and Qmu from internal model.  In external case, need to change the codes to read Qkappa directly, instead of reading the Qp and recomputing the Qkappa (to do). 
Specfem2D:Qkappa and Qmu from external model, and Qkappa and Qmu from internal model.