Dear QDYN developer,

I'm a newer in dynamic simulation. I feel very appricate for you to share the QDYN soft.
I find the results calculated by serial code and parallel code are not identical. Such as "JP_2d_s.m", the Vmax~time relation are extremely different when NPROCS=1,3,4. The same issue also existed in other examples, such as JP_3d.m, JP_3d_l.m.

The fortran compiler is gfortran, with only MPI no openMP.

Looking forward for your reply. Thank you very much!

MPI and serial results are different

import matplotlib 
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
from scipy.ndimage import gaussian_filter

from mpl_toolkits.axes_grid1 import make_axes_locatable

import time 
# Append paths
sys.path.append("/home/gauss/qdyn-release-3.0.0/qdyn/")
from pyqdyn import qdyn

# Instantiate the QDYN class object
p = qdyn()

# Predefine parameters
t_yr = 3600 * 24 * 365.0    # Seconds per year
L = 8e3                     # Length of fault along-strike
L_corner = 1.5e3
W = 4e3                     # Length of fault along-dip

h1 = 2000 
h2 = 4000 

resolution = 7              # Mesh resolution / process zone width
V_SS = 1e-6
V_PL = 1e-10
V_0 = V_PL
dip = 45
f_0 = 0.6
mu = 30e9 
vs = 3000 
nu = mu/vs * 0.5

a0 = 0.002
a1 = 0.025
b = 0.015
dc = 1e-3

sigma = 2E6 


# Get the settings dict
set_dict = p.set_dict

""" Step 1: Define simulation/mesh parameters """
# Global simulation parameters
set_dict["MESHDIM"] = 2        # Simulation dimensionality (1D fault in 2D medium)
set_dict["FAULT_TYPE"] = 2     # Thrust fault
set_dict["TMAX"] = 30*t_yr      # Maximum siut every N steps
set_dict["NTOUT_OT"] = 10			# Temporal interval (number of time steps) for time series output
set_dict["NTOUT_OX"] = 100			# Temporal interval (number of time steps) for snapshot output
set_dict["NXOUT_OX"] = 2			# Spatial interval (number of elements in x-direction) for snapshot output
set_dict["NWOUT_OX"] = 2            # Spatial interval (number of elements in y-direction) for snapshot output
set_dict["V_PL"] = V_PL        # Plate velocity
set_dict["MU"] = mu          # Shear modulus
set_dict["VS"] = vs          # Shear wave speed
set_dict["SIGMA"] = sigma       # Effective normal stress [Pa]
set_dict["FEAT_STRESS_COUPL"] = 1	# Normal stress coupling
set_dict["ACC"] = 1e-7         # Solver accuracy
set_dict["SOLVER"] = 2         # RK
set_dict["Z_CORNER"] = -W*np.sin(dip*np.pi/180)    # Base of the fault (depth taken <0); NOTE: Z_CORNER must be < -W !
set_dict["DIP_W"] = dip         # Dip of the fault
set_dict["NPROC"] = 1         # Parallel
set_dict["MPI_PATH"] = "/usr/bin/mpirun"

# Setting some (default) RSF parameter values
set_dict["SET_DICT_RSF"]["A"] = a0    # Direct effect (will be overwritten later)
set_dict["SET_DICT_RSF"]["B"] = b      # Evolution effect
set_dict["SET_DICT_RSF"]["DC"] = dc     # Characteristic slip distance
set_dict["SET_DICT_RSF"]["V_SS"] = V_SS    # Reference velocity [m/s]
set_dict["SET_DICT_RSF"]["V_0"] = set_dict["V_PL"]     # Initial velocity [m/s]
set_dict["SET_DICT_RSF"]["TH_0"] = 0.99 * set_dict["SET_DICT_RSF"]["DC"] / set_dict["V_PL"]    # Initial (steady-)state [s]
set_dict["SET_DICT_RSF"]["RNS_LAW"] = 2    #  regularized form

# Process zone width [m]
Lb = set_dict["MU"] * set_dict["SET_DICT_RSF"]["DC"] / (set_dict["SET_DICT_RSF"]["B"] * set_dict["SIGMA"])
# Nucleation length [m]
Lc = set_dict["MU"] * set_dict["SET_DICT_RSF"]["DC"] / ((set_dict["SET_DICT_RSF"]["B"] - set_dict["SET_DICT_RSF"]["A"]) * set_dict["SIGMA"])

print(f"Process zone size: {Lb} m \t Nucleation length: {Lc} m")

# Find next power of two for number of mesh elements
Nx = int(np.power(2, np.ceil(np.log2(resolution * L / Lb))))
Nw = int(np.power(2, np.ceil(np.log2(resolution * W / Lb))))

# Spatial coordinate for mesh
xx = np.linspace(-L/2, L/2, Nx, dtype=float)
ww = np.linspace(0, W, Nw, dtype=float)

zz = set_dict["Z_CORNER"] + (ww ) * np.sin(set_dict["DIP_W"] * np.pi / 180.)
XX, ZZ = np.meshgrid(xx, zz)




# Set mesh size and fault length
set_dict["NX"] = Nx
set_dict["NW"] = Nw
set_dict["L"] = L
set_dict["W"] = W 
set_dict["DW"] = W / Nw
# Set time series output node to the middle of the fault
set_dict["IC"] = Nx * (3*Nw // 4) + Nx // 2

IOT1 = [Nx * (3*Nw // 4) + Nx*ii // 6 for ii in range(1,6,3) ] 
IOT2 = [Nx * (1*Nw // 4) + Nx*ii // 6 for ii in range(1,6,3) ] 


set_dict["IOT"] = np.zeros(Nx*Nw)
set_dict["IOT"][IOT1] = 1
set_dict["IOT"][IOT2] = 1


# """ Step 2: Set (default) parameter values and generate mesh """
p.settings(set_dict)
p.render_mesh()

# """ Step 3: override default mesh values """
# # Distribute direct effect a over mesh according to some arbitrary function
scale = 1e3



aa = a0 + (a1 - a0) / (h2-h1) * (ww - h1)
aa[ww<h1] = a0 
aa[ww>h2] = a1
aa = np.repeat(aa, Nx).reshape(Nw,Nx)
aa[p.mesh_dict["X"].reshape(XX.shape)<L_corner] = a1
aa[p.mesh_dict["X"].reshape(XX.shape)>L-L_corner] = a1
aa=aa[::-1]
aa = gaussian_filter(aa, sigma=(5,5))


sigma_n = np.ones(zz.size) * sigma


p.mesh_dict["A"] = aa.ravel()



# Create a 3D plot
fig = plt.figure(figsize = (9,5), clear = True)
ax1 = fig.add_subplot(121, projection='3d')


color_dimension = aa-b # change to desired fourth dimension

minn, maxx = color_dimension.min(), color_dimension.max()
norm = matplotlib.colors.Normalize(minn, maxx)
m = plt.cm.ScalarMappable(norm = norm, cmap='jet')
m.set_array([])
fcolors = m.to_rgba(color_dimension)

im1 = ax1.plot_surface( p.mesh_dict["X"].reshape(XX.shape)/scale, 
                  p.mesh_dict["Y"].reshape(XX.shape)/scale, 
                  p.mesh_dict["Z"].reshape(XX.shape)/scale, 
                  rstride=1, cstride=1, facecolors = fcolors, 
                  linewidth=0, antialiased=False, alpha = 0.5)

ax1.scatter3D(p.mesh_dict["X"][IOT1]/scale, p.mesh_dict["Y"][IOT1]/scale, p.mesh_dict["Z"][IOT1]/scale, color = "k", marker="+",)
ax1.scatter3D(p.mesh_dict["X"][IOT2]/scale, p.mesh_dict["Y"][IOT2]/scale, p.mesh_dict["Z"][IOT2]/scale, color = "k", marker="+")

plt.colorbar(mappable=im1,
              fraction=0.03, pad=0.2, label="a-b", orientation = "horizontal", )
ax1.set_xlabel("X [km]")
ax1.set_ylabel("Y [km]")
ax1.set_zlabel("Z [km]")     

ax3 = fig.add_subplot(122) 


ax3.plot(sigma_n*1e-6, zz*1e-3, )

ax3.set_xlabel("$\\sigma_n$ [MPa]")
ax3.set_ylabel("z [km]")
plt.tight_layout()



target_path = os.path.join(os.getcwd(), "sol_mpi{:.0F}_test_qdyn_r300_a{:.1E}_b{:.1E}_d{:.1E}_dip{:.0f}_zcorner{:.1f}_FINITE{:.0f}"\
                            .format(set_dict["NPROC"],a0,b,dc,dip,set_dict["Z_CORNER"]*1e-3, set_dict["FINITE"])
                            )

if not os.path.exists(target_path):
    os.makedirs(target_path)
os.chdir(target_path)

fig.savefig("setup.png", dpi = 200, bbox_inches="tight")

plt.clf()

# # # # # # # # Write input to qdyn.in
p.write_input()

# MPI 
ti=time.time()
p.run()
tf = time.time()
print(f"Elaspsed time: {tf-ti}")

The serial run results:

This is the mpi run (NPROC=4) that fails to integrate:

I used the below option for the Fortran compilers:
#Only MPI no openMP
F90 = mpif90
OPT = -O3 -Wall

Also, dtau_dt gives 0 values in the output.

Dear developers of QDYN software,

I am a Ph.D. candidate of Earthquake-Seismology at the International Institute of Earthquake Engineering and Seismology (IIEES), Tehran, Iran. I am new in the field of earthquake simulation, so my questions are probably quite naive.

I began learning QDYN by running the "JP_2D_s.m" and “ JP_3d.m” files (without deep asperities) located in the ~/example/Tohoku folder.
To have the same parameters in simulations, I have decreased the solver accuracy of the 3D model to “1e-14”, same as that of 2D model.
My questions are motivated by checking some of the outputs of 2D and 3D simulations.
It seems to me that the results for 2D simulation are different from those of 3D ones. For example, I have compared the recurrence time-periods of the M9.1 event and maximum slip rates (ot.v) for the models.
On the other hand, the results of the 3D-model for different sets of along-strike rupture lengths (p.L) are not stable.

I would be very grateful if you could kindly explain how to fix these problems.

Looking forward to hearing from you.

Best regards,

Mohammad Talebi

I found a problem while I used the pre-compiled executables for Windows.
TOHOKU example JP_3d_l_asp.m script
This is the error I get:
forrtl: severe (170): Program Exception - stack overflow
Image PC Routine Line Source
qdyn.exe 00235E87 Unknown Unknown Unknown
qdyn.exe 001D39E3 SOLVER_mp_SOLVE 29 solver.f90
qdyn.exe 001DA0E9 MAIN 15 main.f90
qdyn.exe 0028D683 Unknown Unknown Unknown
qdyn.exe 00236089 Unknown Unknown Unknown
qdyn.exe 00235F4F Unknown Unknown Unknown
KERNEL32.DLL 74473744 Unknown Unknown Unknown
ntdll.dll 772E9E54 Unknown Unknown Unknown
ntdll.dll 772E9E1F Unknown Unknown Unknown

And the JP_2d_s.m script can run successfully.

I searched this problem on google, it seemed like the program had to increase the stach size.

Dear developers,
how could I set up several asperities ?

In output.f90::ot_init, the following warning was found:

!JPA WARNING VMAX and IOT outputs not implemented yet in parallel

This clearly needs to be addressed...

Hello,

I am trying to reproduce the simulations documented in the paper by Kaneko and 
Ampuero (GRL, 2011) titled: "A mechanism for preseismic steady rupture fronts 
observed in laboratory experiments".  

I believe that I have become relatively acquainted with the newest version of 
QDYN.  I have setup the problem using the parameters provided in the paper as 
follows:

p = qdyn('set');

%fromulation
p.MESHDIM=2;%dimesion (3D)
p.THETA_LAW=2;%slip law = 2,aging law = 1, aging (no heal) = 0
p.RNS_LAW=1;%1=rate-and-state with cutoff
p.V1=100; %cutoff velocities
p.V2=1e-7;%cutoff velocities

%geometry/nodes
p.L=150e-2;
p.W=50e-2;
p.NX=512*4;
p.NW=1;
p.Z_CORNER=-150e-2;
p.N=p.NX*p.NW;
p.DW(1:p.NW)=p.W/p.NW;
p.DIP_W(1:p.NW)=0;


%Initial condition
p.SIGMA=Normal_Stress; %input via function in Matlab
p.V_SS=10e-6;
p.TH_SS=p.DC/p.V_SS;
p.IC=ceil(p.N/2);
p.V_0 =1.01*p.V_SS ;


%BC's
P.FINITE=1;%0=periodic, 1=surrounded by steady slip



%Constitutive parameters
    %I have just discretized the a and b constitutive parameters here so
    %that the left and right hand side of the fault is VS (37.5 cm, 
    %a-b=0.01) and the middle section is VW (75 cm, a-b=-0.0008)

edge_effect=512; 
b_1=0.00*ones(1,edge_effect);
b_2=0.0108.*ones(1,p.N-2*edge_effect);
p.B=[b_1,b_2,b_1];
p.A=ones(1,p.N).*0.01;
p.DC=0.1e-6;
p.MU=32e9;




%time stepping
twm=1.5*0.297E-09;
year = 3600*24*365;
p.TMAX=twm*year;
p.ACC = 1e-14;



%output
p.OX_SEQ=0; %type of fortran file output (0 = only fort.19)
p.NTOUT=1;%temporal interval for snapshot outputs 
p.NXOUT=1;%spatial interval
p.NSTOP=1; %0=stop @TMAX, 1=end of localization , 2=first slip rate peak

My question is:

1. In the paper the boundary conditions are specified as follows:
 "slip rate delta_dot_Load = 5 cm/year is prescribed on the two 75 cm long outer portions of the fault".  How is this specified in the code?  Is this periodic boundary conditions?

I am performing similar laboratory experiments to Nieslen et al. (2010) and 
observe slow propagating fronts prior to nucleation.  I am working with 
Professor Glaser at the University of California in Berkeley. Please feel free 
to email me at [email protected] also.

Regards,

Paul Selvadurai

Original issue reported on code.google.com by [email protected] on 9 Apr 2013 at 6:06

Dear QDYN developers,
In our recent research, we found dynamic weakening could be greatly enhanced by mechanochemical reaction (e.g. flash heating promted by dehydration reaction and amorphization) for clay-bearing fault gouge. For further research, we want to see how this effect would contribute to the kinematics in a fault system. So I am writing to query if the mechanochemical terms (see the Eq. (6) and (7) and associated terms in the attached file) could be added in thermal pressurisation model parameters. Thanks.
Best,
Bowen
Yu et al., 2023.pdf

Currently DW (spacing along-dip in 3D simulations) is a vector, and may spatially vary. When FFT_TYPE == 2 (FFT taken along-dip rather than along-strike), DW must be constant (DW = W / NW). If contributions from the fault curvature are ignored, DW should be W / NW at all times, regardless of the value(s) of DIP_W.

Sorry, Nothing is wrong with the derivatives. My bad!

The restarted simulation differs from the seamless result. I put an example from the QDYN documentation 'Single asperity simulations (3D)' by slightly changing the input.

import matplotlib.pyplot as plt
import numpy as np

import os
import sys

ver = 'qdyn3_es1' # version of the qdyn
# Add QDYN source directory to PATH
qdyn_dir = os.path.join(os.path.join(os.path.expanduser("~"),'qdyn_v'), ver) 
src_dir = os.path.join(qdyn_dir, 'qdyn')
utils_dir= os.path.join(src_dir, 'utils')

# Append src directory to Python path
sys.path.append(src_dir)
# Import QDYN
from pyqdyn import qdyn


# Working Directory
working_dir = os.path.join(os.path.join(os.path.expanduser("~"),'Documents'), f'test3_{ver}') 
if not os.path.exists( working_dir ) :
    os.system(f'mkdir {working_dir}')


# Instantiate the QDYN class object
p = qdyn()

# Predefine parameters
t_yr = 3600 * 24 * 365.0    # Seconds per year
L = 5e3                     # Length of fault along-strike
W = 5e3                     # Length of fault along-dip
resolution = 5              # Mesh resolution / process zone width

# Get the settings dict
set_dict = p.set_dict

""" Step 1: Define simulation/mesh parameters """
# Global simulation parameters
set_dict["MESHDIM"] = 2        # Simulation dimensionality (1D fault in 2D medium)
set_dict["FAULT_TYPE"] = 2     # Thrust fault
set_dict["TMAX"] = 350*t_yr      # Maximum simulation time [s]
set_dict["NTOUT_LOG"] = 10          # Temporal interval (number of time steps) for time series output
set_dict["NTOUT_OT"] = 10			# Temporal interval (number of time steps) for time series output
set_dict["NTOUT_OX"] = 10			# Temporal interval (number of time steps) for snapshot output
set_dict["NXOUT_OX"] = 1			# Spatial interval (number of elements in x-direction) for snapshot output
set_dict["NWOUT_OX"] = 1            # Spatial interval (number of elements in y-direction) for snapshot output
set_dict["V_PL"] = 1e-11        # Plate velocity
set_dict["MU"] = 3e10          # Shear modulus
set_dict["SIGMA"] = 1e7        # Effective normal stress [Pa]
set_dict["ACC"] = 1e-7         # Solver accuracy
set_dict["SOLVER"] = 2         # Solver type (Runge-Kutta)
set_dict["Z_CORNER"] = -0.25e4    # Base of the fault (depth taken <0); NOTE: Z_CORNER must be < -W !
set_dict["DIP_W"] = 30         # Dip of the fault
set_dict["FEAT_STRESS_COUPL"] = 1	# Normal stress coupling


# Setting some (default) RSF parameter values
set_dict["SET_DICT_RSF"]["A"] = 0.2e-2    # Direct effect (will be overwritten later)
set_dict["SET_DICT_RSF"]["B"] = 1e-2      # Evolution effect
set_dict["SET_DICT_RSF"]["DC"] = 0.8e-3     # Characteristic slip distance
set_dict["SET_DICT_RSF"]["V_SS"] = set_dict["V_PL"]    # Reference velocity [m/s]
set_dict["SET_DICT_RSF"]["V_0"] = set_dict["V_PL"]     # Initial velocity [m/s]
set_dict["SET_DICT_RSF"]["TH_0"] = 0.99 * set_dict["SET_DICT_RSF"]["DC"] / set_dict["V_PL"]    # Initial (steady-)state [s]

# Process zone width [m]
Lb = set_dict["MU"] * set_dict["SET_DICT_RSF"]["DC"] / (set_dict["SET_DICT_RSF"]["B"] * set_dict["SIGMA"])
# Nucleation length [m]
Lc = set_dict["MU"] * set_dict["SET_DICT_RSF"]["DC"] / ((set_dict["SET_DICT_RSF"]["B"] - set_dict["SET_DICT_RSF"]["A"]) * set_dict["SIGMA"])

print(f"Process zone size: {Lb} m \t Nucleation length: {Lc} m")

# Find next power of two for number of mesh elements
Nx = int(np.power(2, np.ceil(np.log2(resolution * L / Lb))))
Nw = int(np.power(2, np.ceil(np.log2(resolution * W / Lb))))

# Spatial coordinate for mesh
x = np.linspace(-L/2, L/2, Nx, dtype=float)
z = np.linspace(-W/2, W/2, Nw, dtype=float)
X, Z = np.meshgrid(x, z)
z = -(set_dict["Z_CORNER"] + (z + W/2) * np.cos(set_dict["DIP_W"] * np.pi / 180.))

# Set mesh size and fault length
set_dict["NX"] = Nx
set_dict["NW"] = Nw
set_dict["L"] = L
set_dict["W"] = W 
set_dict["DW"] = W / Nw
# Set time series output node to the middle of the fault
set_dict["IC"] = Nx * (Nw // 2) + Nx // 2

# Working dir
os.chdir(working_dir)
working_dir = os.path.join(working_dir,
                               '{}__a{:.1E}_b{:.1E}_d{:.1E}_TMAX{:.1E}'.format(
                                   set_dict["FEAT_STRESS_COUPL"],
                                   set_dict["SET_DICT_RSF"]["A"],
                                   set_dict["SET_DICT_RSF"]["B"],
                                   set_dict["SET_DICT_RSF"]["DC"],
                                   set_dict["TMAX"]/t_yr,
                                  
                                   )
    )

if not os.path.exists(working_dir):
    os.system(f'mkdir {working_dir}')
os.chdir(working_dir)


for mod in [False, True]:
    
    if mod == True: # restart mod
        
        os.system('cp -rvf {} {}_oo'.format(working_dir, working_dir))
        
        set_dict["TMAX"] += (350*t_yr)
        
        
    """ Step 2: Set (default) parameter values and generate mesh """
    p.settings(set_dict)
    p.render_mesh()
    
    """ Step 3: override default mesh values """
    # Distribute direct effect over mesh according to some arbitrary function
    scale = 1e3
    p.mesh_dict["A"] = 2 * set_dict["SET_DICT_RSF"]["B"] * (1 - 0.9*np.exp(- (X**2 + Z**2) / (2 * scale**2))).ravel()
    
    
        
    # Write input to qdyn.in
    p.write_input()
    
    # QDYN RUN
    # os.system('{}'.format(os.path.join(src_dir,'qdyn')))
    p.run(restart=mod)

The above code runs for 350 years and then restarts for an additional 350 years. I also run a continuous 700 years solution. Then I compared the results. Here is a 'output_vmax' file result.

There is a constant shift between the seamless and

restarted solutions.

The amount of shift between the solutions differs w.r.t the initial values and the inputs.

Currently QDYN internally generates a mesh, except for 3D MPI simulations. Since the Python wrapper already pre-computes a mesh with default settings, we might as well shift the burden from FORTRAN to Python/Matlab, and always read the mesh from the input file.

To-do

• Verify that the default mesh generator in Python outputs the same results as the one in FORTRAN
• Create Python routines for standard fault configurations
• Create an internal check that the along-strike direction is co-linear for each row of the mesh, and that the number of elements in that direction is an integer power of 2 (both are required by the FFT)
• Design an integration test for an unambiguous fault geometry (i.e. mesh in == mesh out)

What steps will reproduce the problem?
1.  By setting p.FINITE = 1 (constant slip rate) the program fails.
2.
3.

What is the expected output? What do you see instead?
ERROR in matlab when trying to call the qdyn.exe script

What version of the product are you using? On what operating system?
Windows 8 and Windows 7 / latest version on QDYN

Please provide any additional information below.


Hello,
This is the error I get:

 Start reading input: ... 
 Input complete 
 Initializing mesh ... 
 1D fault, uniform grid 
 Initializing parameters: ... 
 impedance =    1226666.66666667      
 Intializing kernel: ... 
 FFT applied 
forrtl: severe (29): file not found, unit 57, file 
C:\Users\pa.selvadurai\Desktop\SLIP\~\2D_RUPTURE\STATIC\Matlab\kernel_I_32768.ta
b  
Image              PC        Routine            Line        Source              
qdyn.exe           00A5623A  Unknown               Unknown  Unknown 
qdyn.exe           009AD13A  Unknown               Unknown  Unknown 
qdyn.exe           009AC33B  Unknown               Unknown  Unknown 
qdyn.exe           009BC176  Unknown               Unknown  Unknown 
qdyn.exe           0098A5AC  _FAULT_STRESS_mp_         109  fault_stress.f90 
qdyn.exe           00989D85  _FAULT_STRESS_mp_          59  fault_stress.f90 
qdyn.exe           009AA022  _INITIALIZE_mp_IN          81  initialize.f90 
qdyn.exe           009AA0D9  _MAIN__                    14  main.f90 
qdyn.exe           00A5D683  Unknown               Unknown  Unknown 
qdyn.exe           00A06089  Unknown               Unknown  Unknown 
qdyn.exe           00A05F4F  Unknown               Unknown  Unknown 
KERNEL32.DLL       754B8543  Unknown               Unknown  Unknown 
ntdll.dll          77DAAC69  Unknown               Unknown  Unknown 
ntdll.dll          77DAAC3C  Unknown               Unknown  Unknown 
Error using textread (line 166)
File not found.

Error in qdyn>read_qdyn_out (line 445)
[ot.t, ot.locl, ot.cl, ot.p, ot.pdot, ...

Error in qdyn (line 328)
    [ot,ox]= read_qdyn_out(NAME);

Error in SLIP_LAW (line 188)
[p,ot0,ox0]  = qdyn('run',p);

Original issue reported on code.google.com by [email protected] on 16 Apr 2013 at 6:09

See the documentation. Under "Discretization and accuracy parameters", the table entry NPROCS should be NPROC

Hi Team,

why is number of mesh elements N:

Find next power of two for number of mesh elements

N = int(np.power(2, np.ceil(np.log2(resolution * L / Lb)))) ?

Where is equation for it?

Thank you!

Since the QDYN-SPECFEM bridge (QSB) hasn't been used in a very long time, it is unlikely that it still works. We need to update and test the code, and integrate it in the standard test suite to ensure its operability. A first step would be to create a simple example that we can run to test the code.

The spatial resolution of the snapshot output (OX) is controlled by NXOUT, but this only applies to the x-coordinate. The output resolution in the y-direction is the mesh resolution. It would be better to also have a NYOUT. Then, instead of:

do ixout=1, nn, pb%ox%nxout
    ...

one should have:

do ixout=1, nx, pb%ox%nxout
    do iyout=1, ny, pb%ox%nyout
        nout = (ixout - 1) * ny + iyout
        ...

in the OX module of output.f90 (currently being revised in #43)

Dear co-developer I noticed that the fault type (strike-slip, thrust, etc) is hard-coded as constant (in constants.f90), i.e. each time users of QDYN want to use a different fault type they will have to recompile the code (and it is not properly documented in the manual). Is there a good reason for doing so (e.g. computation efficiency, etc.)? If not we need to fix that as a run-time input.

Yingdi

As @gpercy and @jpampuero mentioned to me, the current FFT routine requires a vector of size 2^N, though there are no internal checks for this requirement, possibly resulting in corrupt/inaccurate data. Possible workarounds:

1. Add a sanity check while the input file is being read, to make sure that the mesh size is a power of 2
2. Allow for arbitrary mesh sizes, but pad the data vector with zeros

The latter is more user-friendly, while the former is computationally more efficient.

The variable iot_buf is intended to be used only by the master processor, but this seems not be the case in the module output.f90. This might be the reason behind the corruption of arrays in the module fault_stress.f90 when trying to run a simulation with multiple processors. One possible fix is to enclose L322 to L354 in the master branch between an if (is_MPI_master()) then... endif

The timeseries output is no longer up to date, and contains numerous code duplicates (see also #42). The nested code structure needs to be unravelled and generalised for serial and parallel execution.

Since QDYN supports Python, as an open-source alternative to MATLAB, support for Octave is not really needed. Also, QDYN hasn't been tested with Octave for years, so that its functioning cannot be guaranteed. Dropping Octave support will remove some clutter in output.f90 and reduces the complexity of refurbishing the time series output (OT) module (see #45).

Dear QDYN developers,

I found another issue during calculations. For the example 'JP_3d_l_asp.m', 'ox1' is empty for serial code and parallel code.

Looking forward for your reply. Thank you very much!

Hi, I am trying to install QDYN following the instructions in the documentation, but I'm running with some problems in the installation. I'm using the intel_oneAPI compiler (ifort) and openMP. My first question is related to setting the parameters in the section “User Settings” of constants.f90. I left all the variables with the default values. Should I make any change there?
Secondly, when I run make I get the following error:

ifort -O3 -ip -ipo -parallel -qopenmp -par-threshold:50   -o qdyn unittests_rsf.o friction.o mesh.o constants.o fftsg.o ode_bs.o diffusion_solver.o friction_cns.o problem_class.o parallel.o main.o solver.o unittests.o derivs_all.o input.o utils.o okada.o unittests_aux.o fault_stress.o initialize.o output.o ode_rk45.o ode_rk45_2.o
ipo: warning #11003: no IR in object file unittests_rsf.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file friction.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file mesh.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file constants.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file fftsg.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file ode_bs.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file diffusion_solver.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file friction_cns.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file problem_class.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file parallel.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file main.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file solver.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file unittests.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file derivs_all.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file input.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file utils.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file okada.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file unittests_aux.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file fault_stress.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file initialize.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file output.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file ode_rk45.o; was the source file compiled with -ipo
ipo: warning #11003: no IR in object file ode_rk45_2.o; was the source file compiled with -ipo
parallel.o: In function `my_mpi_mp_finalize_mpi_':
parallel.f90:(.text+0x1b): undefined reference to `mpi_barrier_'
parallel.f90:(.text+0x24): undefined reference to `mpi_finalize_'
parallel.o: In function `my_mpi_mp_gather_allv_':
parallel.f90:(.text+0x14c): undefined reference to `mpi_allgatherv_'
parallel.o: In function `my_mpi_mp_gather_allvdouble_':
parallel.f90:(.text+0x18c): undefined reference to `mpi_allgatherv_'
parallel.o: In function `my_mpi_mp_gather_allvdouble_root_':
parallel.f90:(.text+0x1c2): undefined reference to `mpi_gatherv_'
parallel.o: In function `my_mpi_mp_gather_allvi_root_':
parallel.f90:(.text+0x1f2): undefined reference to `mpi_gatherv_'
parallel.o: In function `my_mpi_mp_gather_alli_':
parallel.f90:(.text+0x221): undefined reference to `mpi_allgather_'
parallel.o: In function `my_mpi_mp_synchronize_all_':
parallel.f90:(.text+0x23b): undefined reference to `mpi_barrier_'
parallel.o: In function `my_mpi_mp_sum_allreduce_':
parallel.f90:(.text+0x301): undefined reference to `mpi_allreduce_'
parallel.o: In function `my_mpi_mp_max_allproc_':
parallel.f90:(.text+0x3b6): undefined reference to `mpi_allreduce_'
parallel.o: In function `my_mpi_mp_min_allproc_':
parallel.f90:(.text+0x3e6): undefined reference to `mpi_allreduce_'
parallel.o: In function `my_mpi_mp_init_mpi_':
parallel.f90:(.text+0x3fa): undefined reference to `mpi_init_'
parallel.f90:(.text+0x40e): undefined reference to `mpi_comm_size_'
parallel.f90:(.text+0x422): undefined reference to `mpi_comm_rank_'
parallel.f90:(.text+0x43c): undefined reference to `mpi_finalize_'
make: *** [qdyn] Error 1

Do you know what could be the source of error? I attached the Makefile for reference.
I am not sure if it's useful, but the ifort version is 2021.2.0 20210228 and the Intel MPI library Version is 2021.2 Build 20210302.
Thanks in advance!
Constanza
Makefile.txt

Dear QDYN developers

Is there any scaling relations between the parameter w (half-width of shear distribution [m]) and the length of element, or any requirement should be meet? All of the thermal pressurization parameters should be set for the whole model or only part of the model?

Looking forward for your reply! Thank you!

I downloaded the QDYN software on GitHub, and installed OpenMPI, and compiled the QDYN successfully with the Makefile (shown below).

GNU, serial

#Only MPI no openMP
F90 = mpif90
#OPT = -O3 -Wall

For parallel runs with MPI+OpenMP:

#F90 = mpif90
OPT = -O3 -Wall -fopenmp
#OPT = -O2 -w
#-- Intel Fortran --
#F90 = ifort
#OPT = -O3 -ip -ipo

For parallel runs with OpenMP:

#OPT = -O3 -ip -ipo -parallel -openmp -par-threshold:50
#-no-inline-factor
#OPT = -O3 -openmp
#OPT = -O0 -warn all -C
#OPT = -O0 -warn all -parallel -par-report3 -C
#OPT = -O0 -g
#OPT = -O3 -ip -ipo -arch SSE2 -tpp7 -xN # pablo.ethz.ch
#OPT = -O3 -ip -ipo -xB # lapfabio
#OPT = -O3 -ip -ipo -unroll

#-- Lahey Fortran95 --
#F90 = lf95
#OPT = -O
#OPT = --chk aesux

#-- Digital --
#F90 = f95
#OPT = -fast

3. Set the compiler option that enables preprocessing in your compiler

ifort

#PREPROC = -fpp

gfortran

PREPROC = -cpp

I only changed the OPT = -O3 -Wall -fopenmp and PREPROC = -cpp parameters, all other parameters were not changed.

It can run successfully with 1D fault examples, but almost all the 2D fault examples cannot run successfully, I do not know why.

example test_3dfft.m
This is the error I get:
filename =test_2d_fft_ab0.3L80nx64W80nw64z-100.mat
Lb_over_dx = 1.9200
Lc=3428.5714 L/Lc=23.3333 W/Lc=23.3333
Linf=3118.1377 L/Linf=25.6563 W/Linf=25.6563
MESHDIM = 2
Number of processors = 1
Start reading input: ...
Mesh input complete
Flags input complete
Input complete
Initializing mesh ...
2D fault, uniform grid along-strike
Impedance = 5000000.0000000000
Intializing kernel: ...
Generating 3D kernel...
OouraFFT applied along-strike
Kernel intialized
Initialization completed
Values at selected point of the fault:
it, dt (secs), time (yrs), v_max (m/s), sigma_max (MPa)
0 0.000E+00 0.000E+00 0.101E-08 0.500E+00
*** glibc detected *** /home/x/qdyn/src/qdyn: free(): invalid next size (fast): 0x0000000000e8e9a0 ***
======= Backtrace: =========
/lib64/libc.so.6[0x3a15876166]
/lib64/libc.so.6[0x3a15878ca3]
/home/x/qdyn/src/qdyn[0x408700]
/home/x/qdyn/src/qdyn[0x408c20]
/home/x/qdyn/src/qdyn[0x401b63]
/home/x/qdyn/src/qdyn[0x4301af]
/home/x/qdyn/src/qdyn[0x43104f]
/home/x/qdyn/src/qdyn[0x41b4f7]
/home/x/qdyn/src/qdyn[0x4319ea]
/lib64/libc.so.6(__libc_start_main+0xfd)[0x3a1581ed1d]
/home/x/qdyn/src/qdyn[0x401789]
======= Memory map: ========
00400000-00435000 r-xp 00000000 08:21 715194528 /home/x/qdyn/src/qdyn
00635000-00636000 rw-p 00035000 08:21 715194528 /home/x/qdyn/src/qdyn
00636000-00637000 rw-p 00000000 00:00 0
00e52000-0106d000 rw-p 00000000 00:00 0 [heap]
3a15400000-3a15420000 r-xp 00000000 fd:00 786857 /lib64/ld-2.12.so
3a1561f000-3a15620000 r--p 0001f000 fd:00 786857 /lib64/ld-2.12.so
3a15620000-3a15621000 rw-p 00020000 fd:00 786857 /lib64/ld-2.12.so
3a15621000-3a15622000 rw-p 00000000 00:00 0
3a15800000-3a158e0000 r-xp 00000000 fd:00 786879 /lib64/libc-2.12.so
3a158e0000-3a158e1000 r-xp 000e0000 fd:00 786879 /lib64/libc-2.12.so
3a158e1000-3a158e5000 r-xp 000e1000 fd:00 786879 /lib64/libc-2.12.so
3a158e5000-3a158e6000 r-xp 000e5000 fd:00 786879 /lib64/libc-2.12.so
3a158e6000-3a158e9000 r-xp 000e6000 fd:00 786879 /lib64/libc-2.12.so
3a158e9000-3a158ea000 r-xp 000e9000 fd:00 786879 /lib64/libc-2.12.so
3a158ea000-3a158eb000 r-xp 000ea000 fd:00 786879 /lib64/libc-2.12.so
3a158eb000-3a158ec000 r-xp 000eb000 fd:00 786879 /lib64/libc-2.12.so
3a158ec000-3a1598b000 r-xp 000ec000 fd:00 786879 /lib64/libc-2.12.so
3a1598b000-3a15b8a000 ---p 0018b000 fd:00 786879 /lib64/libc-2.12.so
3a15b8a000-3a15b8e000 r--p 0018a000 fd:00 786879 /lib64/libc-2.12.so
3a15b8e000-3a15b8f000 rw-p 0018e000 fd:00 786879 /lib64/libc-2.12.so
3a15b8f000-3a15b94000 rw-p 00000000 00:00 0
3a15c00000-3a15c83000 r-xp 00000000 fd:00 786897 /lib64/libm-2.12.so
3a15c83000-3a15e82000 ---p 00083000 fd:00 786897 /lib64/libm-2.12.so
3a15e82000-3a15e83000 r--p 00082000 fd:00 786897 /lib64/libm-2.12.so
3a15e83000-3a15e84000 rw-p 00083000 fd:00 786897 /lib64/libm-2.12.so
3a16000000-3a16017000 r-xp 00000000 fd:00 786890 /lib64/libpthread-2.12.so
3a16017000-3a16217000 ---p 00017000 fd:00 786890 /lib64/libpthread-2.12.so
3a16217000-3a16218000 r--p 00017000 fd:00 786890 /lib64/libpthread-2.12.so
3a16218000-3a16219000 rw-p 00018000 fd:00 786890 /lib64/libpthread-2.12.so
3a16219000-3a1621d000 rw-p 00000000 00:00 0
3a16400000-3a16402000 r-xp 00000000 fd:00 786946 /lib64/libdl-2.12.so
3a16402000-3a16602000 ---p 00002000 fd:00 786946 /lib64/libdl-2.12.so
3a16602000-3a16603000 r--p 00002000 fd:00 786946 /lib64/libdl-2.12.so
3a16603000-3a16604000 rw-p 00003000 fd:00 786946 /lib64/libdl-2.12.so
3a16c00000-3a16c07000 r-xp 00000000 fd:00 786891 /lib64/librt-2.12.so
3a16c07000-3a16e06000 ---p 00007000 fd:00 786891 /lib64/librt-2.12.so
3a16e06000-3a16e07000 r--p 00006000 fd:00 786891 /lib64/librt-2.12.so
3a16e07000-3a16e08000 rw-p 00007000 fd:00 786891 /lib64/librt-2.12.so
3a26000000-3a26002000 r-xp 00000000 fd:00 786971 /lib64/libutil-2.12.so
3a26002000-3a26201000 ---p 00002000 fd:00 786971 /lib64/libutil-2.12.so
3a26201000-3a26202000 r--p 00001000 fd:00 786971 /lib64/libutil-2.12.so
3a26202000-3a26203000 rw-p 00002000 fd:00 786971 /lib64/libutil-2.12.so
3e4de00000-3e4de0d000 r-xp 00000000 fd:00 1744494 /usr/lib64/libgomp.so.1.0.0
3e4de0d000-3e4e00c000 ---p 0000d000 fd:00 1744494 /usr/lib64/libgomp.so.1.0.0
3e4e00c000-3e4e00d000 rw-p 0000c000 fd:00 1744494 /usr/lib64/libgomp.so.1.0.0
7fb25c000000-7fb25c021000 rw-p 00000000 00:00 0
7fb25c021000-7fb260000000 ---p 00000000 00:00 0
7fb264000000-7fb264021000 rw-p 00000000 00:00 0
7fb264021000-7fb268000000 ---p 00000000 00:00 0
7fb268000000-7fb268021000 rw-p 00000000 00:00 0
7fb268021000-7fb26c000000 ---p 00000000 00:00 0
7fb26c000000-7fb26c021000 rw-p 00000000 00:00 0
7fb26c021000-7fb270000000 ---p 00000000 00:00 0
7fb270000000-7fb270021000 rw-p 00000000 00:00 0
7fb270021000-7fb274000000 ---p 00000000 00:00 0
7fb274000000-7fb274021000 rw-p 00000000 00:00 0
7fb274021000-7fb278000000 ---p 00000000 00:00 0
7fb278000000-7fb278021000 rw-p 00000000 00:00 0
7fb278021000-7fb27c000000 ---p 00000000 00:00 0
7fb27c000000-7fb27c021000 rw-p 00000000 00:00 0
7fb27c021000-7fb280000000 ---p 00000000 00:00 0
7fb280000000-7fb280021000 rw-p 00000000 00:00 0
7fb280021000-7fb284000000 ---p 00000000 00:00 0
7fb284000000-7fb284021000 rw-p 00000000 00:00 0
7fb284021000-7fb288000000 ---p 00000000 00:00 0
7fb288000000-7fb288021000 rw-p 00000000 00:00 0
7fb288021000-7fb28c000000 ---p 00000000 00:00 0
7fb28c000000-7fb28c021000 rw-p 00000000 00:00 0
7fb28c021000-7fb290000000 ---p 00000000 00:00 0
7fb290000000-7fb290021000 rw-p 00000000 00:00 0
7fb290021000-7fb294000000 ---p 00000000 00:00 0
7fb296ff8000-7fb296ff9000 ---p 00000000 00:00 0
7fb296ff9000-7fb2979f9000 rw-p 00000000 00:00 0
7fb2979f9000-7fb2979fa000 ---p 00000000 00:00 0
7fb2979fa000-7fb2983fa000 rw-p 00000000 00:00 0
7fb2983fa000-7fb2983fb000 ---p 00000000 00:00 0
7fb2983fb000-7fb298dfb000 rw-p 00000000 00:00 0
7fb298dfb000-7fb298dfc000 ---p 00000000 00:00 0
7fb298dfc000-7fb2997fc000 rw-p 00000000 00:00 0
7fb2997fc000-7fb2997fd000 ---p 00000000 00:00 0
7fb2997fd000-7fb29a1fd000 rw-p 00000000 00:00 0
7fb29a1fd000-7fb29a1fe000 ---p 00000000 00:00 0
7fb29a1fe000-7fb29abfe000 rw-p 00000000 00:00 0
7fb29abfe000-7fb29abff000 ---p 00000000 00:00 0
7fb29abff000-7fb29b5ff000 rw-p 00000000 00:00 0
7fb29b5ff000-7fb29b600000 ---p 00000000 00:00 0
7fb29b600000-7fb29c000000 rw-p 00000000 00:00 0
7fb29c000000-7fb29c021000 rw-p 00000000 00:00 0
7fb29c021000-7fb2a0000000 ---p 00000000 00:00 0
7fb2a0695000-7fb2a0696000 ---p 00000000 00:00 0
7fb2a0696000-7fb2a1096000 rw-p 00000000 00:00 0
7fb2a1096000-7fb2a109b000 r-xp 00000000 08:21 55607473 /home/x/soft/openmpi/lib/openmpi/mca_osc_sm.so
7fb2a109b000-7fb2a129a000 ---p 00005000 08:21 55607473 /home/x/soft/openmpi/lib/openmpi/mca_osc_sm.so
7fb2a129a000-7fb2a129c000 rw-p 00004000 08:21 55607473 /home/x/soft/openmpi/lib/openmpi/mca_osc_sm.so
7fb2a129c000-7fb2a12b4000 r-xp 00000000 08:21 55607477 /home/x/soft/openmpi/lib/openmpi/mca_osc_rdma.so
7fb2a12b4000-7fb2a14b4000 ---p 00018000 08:21 55607477 /home/x/soft/openmpi/lib/openmpi/mca_osc_rdma.so
7fb2a14b4000-7fb2a14b5000 rw-p 00018000 08:21 55607477 /home/x/soft/openmpi/lib/openmpi/mca_osc_rdma.so
7fb2a14b5000-7fb2a14c9000 r-xp 00000000 08:21 55607475 /home/x/soft/openmpi/lib/openmpi/mca_osc_pt2pt.so
7fb2a14c9000-7fb2a16c8000 ---p 00014000 08:21 55607475 /home/x/soft/openmpi/lib/openmpi/mca_osc_pt2pt.so
7fb2a16c8000-7fb2a16ca000 rw-p 00013000 08:21 55607475 /home/x/soft/openmpi/lib/openmpi/mca_osc_pt2pt.so
7fb2a1706000-7fb2a1707000 ---p 00000000 00:00 0
7fb2a1707000-7fb2a2107000 rw-p 00000000 00:00 0
7fb2a2107000-7fb2a210d000 r-xp 00000000 08:21 55607491 /home/x/soft/openmpi/lib/openmpi/mca_vprotocol_pessimist.so
7fb2a210d000-7fb2a230d000 ---p 00006000 08:21 55607491 /home/x/soft/openmpi/lib/openmpi/mca_vprotocol_pessimist.so
7fb2a230d000-7fb2a230e000 rw-p 00006000 08:21 55607491 /home/x/soft/openmpi/lib/openmpi/mca_vprotocol_pessimist.so
7fb2a2bfe000-7fb2a2bff000 ---p 00000000 00:00 0
7fb2a2bff000-7fb2a35ff000 rw-p 00000000 00:00 0
7fb2a35ff000-7fb2a3600000 ---p 00000000 00:00 0
7fb2a3600000-7fb2a4000000 rw-p 00000000 00:00 0
7fb2a4000000-7fb2a4021000 rw-p 00000000 00:00 0
7fb2a4021000-7fb2a8000000 ---p 00000000 00:00 0
7fb2a8140000-7fb2a8141000 ---p 00000000 00:00 0
7fb2a8141000-7fb2a8b41000 rw-p 00000000 00:00 0
7fb2a8b41000-7fb2a8b42000 ---p 00000000 00:00 0
7fb2a8b42000-7fb2a9542000 rw-p 00000000 00:00 0
7fb2a9dac000-7fb2a9dad000 ---p 00000000 00:00 0
7fb2a9dad000-7fb2aa7ad000 rw-p 00000000 00:00 0
7fb2aa7ad000-7fb2aa7ae000 r-xp 00000000 08:21 55607319 /home/x/soft/openmpi/lib/openmpi/mca_patcher_overwrite.so
7fb2aa7ae000-7fb2aa9ae000 ---p 00001000 08:21 55607319 /home/x/soft/openmpi/lib/openmpi/mca_patcher_overwrite.so
7fb2aa9ae000-7fb2aa9af000 rw-p 00001000 08:21 55607319 /home/x/soft/openmpi/lib/openmpi/mca_patcher_overwrite.so
7fb2aa9af000-7fb2aa9b3000 rw-p 00000000 00:00 0
7fb2aa9b3000-7fb2aa9e8000 r-xp 00000000 fd:02 2364996 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libquadmath.so.0
7fb2aa9e8000-7fb2aabe7000 ---p 00035000 fd:02 2364996 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libquadmath.so.0
7fb2aabe7000-7fb2aabe8000 rw-p 00034000 fd:02 2364996 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libquadmath.so.0
7fb2aabe8000-7fb2aabe9000 rw-p 00000000 00:00 0
7fb2aabe9000-7fb2aacbf000 r-xp 00000000 08:21 55050249 /home/x/soft/openmpi/lib/libopen-pal.so.20.1.0
7fb2aacbf000-7fb2aaebf000 ---p 000d6000 08:21 55050249 /home/x/soft/openmpi/lib/libopen-pal.so.20.1.0
7fb2aaebf000-7fb2aaec8000 rw-p 000d6000 08:21 55050249 /home/x/soft/openmpi/lib/libopen-pal.so.20.1.0
7fb2aaec8000-7fb2aaed0000 rw-p 00000000 00:00 0
7fb2aaed0000-7fb2aaf4a000 r-xp 00000000 08:21 55050257 /home/x/soft/openmpi/lib/libopen-rte.so.20.0.0
7fb2aaf4a000-7fb2ab149000 ---p 0007a000 08:21 55050257 /home/x/soft/openmpi/lib/libopen-rte.so.20.0.0
7fb2ab149000-7fb2ab14e000 rw-p 00079000 08:21 55050257 /home/x/soft/openmpi/lib/libopen-rte.so.20.0.0
7fb2ab14e000-7fb2ab150000 rw-p 00000000 00:00 0
7fb2ab150000-7fb2ab165000 r-xp 00000000 fd:02 2365006 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgcc_s.so.1
7fb2ab165000-7fb2ab364000 ---p 00015000 fd:02 2365006 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgcc_s.so.1
7fb2ab364000-7fb2ab365000 rw-p 00014000 fd:02 2365006 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgcc_s.so.1
7fb2ab365000-7fb2ab366000 rw-p 00000000 00:00 0
7fb2ab380000-7fb2ab381000 rw-p 00000000 00:00 0
7fb2ab381000-7fb2ab493000 r-xp 00000000 fd:02 2365007 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgfortran.so.3
7fb2ab493000-7fb2ab692000 ---p 00112000 fd:02 2365007 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgfortran.so.3
7fb2ab692000-7fb2ab694000 rw-p 00111000 fd:02 2365007 /usr/local/MATLAB/R2014a/sys/os/glnxa64/libgfortran.so.3
7fb2ab694000-7fb2ab751000 r-xp 00000000 08:21 55050263 /home/x/soft/openmpi/lib/libmpi.so.20.0.1
7fb2ab751000-7fb2ab951000 ---p 000bd000 08:21 55050263 /home/x/soft/openmpi/lib/libmpi.so.20.0.1
7fb2ab951000-7fb2ab961000 rw-p 000bd000 08:21 55050263 /home/x/soft/openmpi/lib/libmpi.so.20.0.1
7fb2ab961000-7fb2ab973000 rw-p 00000000 00:00 0
7fb2ab973000-7fb2ab9c3000 r-xp 00000000 08:21 55050683 /home/x/soft/openmpi/lib/libmpi_mpifh.so.20.0.0
7fb2ab9c3000-7fb2abbc3000 ---p 00050000 08:21 55050683 /home/x/soft/openmpi/lib/libmpi_mpifh.so.20.0.0
7fb2abbc3000-7fb2abbc4000 rw-p 00050000 08:21 55050683 /home/x/soft/openmpi/lib/libmpi_mpifh.so.20.0.0
7fb2abbc4000-7fb2abbc6000 rw-p 00000000 00:00 0
7fb2abbc6000-7fb2abbc8000 r-xp 00000000 08:21 55050687 /home/x/soft/openmpi/lib/libmpi_usempi.so.20.0.0
7fb2abbc8000-7fb2abdc7000 ---p 00002000 08:21 55050687 /home/x/soft/openmpi/lib/libmpi_usempi.so.20.0.0
7fb2abdc7000-7fb2abdc8000 rw-p 00001000 08:21 55050687 /home/x/soft/openmpi/lib/libmpi_usempi.so.20.0.0
7fb2abdc8000-7fb2abdc9000 rw-p 00000000 00:00 0
7fffba14f000-7fffba165000 rw-p 00000000 00:00 0 [stack]
7fffba1b3000-7fffba1b4000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]
[CJXU-GEO:32335] *** Process received signal ***
[CJXU-GEO:32335] Signal: Aborted (6)
[CJXU-GEO:32335] Signal code: (-6)
[CJXU-GEO:32335] [ 0] /lib64/libpthread.so.0[0x3a1600f710]
[CJXU-GEO:32335] [ 1] /lib64/libc.so.6(gsignal+0x35)[0x3a15832925]
[CJXU-GEO:32335] [ 2] /lib64/libc.so.6(abort+0x175)[0x3a15834105]
[CJXU-GEO:32335] [ 3] /lib64/libc.so.6[0x3a15870837]
[CJXU-GEO:32335] [ 4] /lib64/libc.so.6[0x3a15876166]
[CJXU-GEO:32335] [ 5] /lib64/libc.so.6[0x3a15878ca3]
[CJXU-GEO:32335] [ 6] /home/x/qdyn/src/qdyn[0x408700]
[CJXU-GEO:32335] [ 7] /home/x/qdyn/src/qdyn[0x408c20]
[CJXU-GEO:32335] [ 8] /home/x/qdyn/src/qdyn[0x401b63]
[CJXU-GEO:32335] [ 9] /home/x/qdyn/src/qdyn[0x4301af]
[CJXU-GEO:32335] [10] /home/x/qdyn/src/qdyn[0x43104f]
[CJXU-GEO:32335] [11] /home/x/qdyn/src/qdyn[0x41b4f7]
[CJXU-GEO:32335] [12] /home/x/qdyn/src/qdyn[0x4319ea]
[CJXU-GEO:32335] [13] /lib64/libc.so.6(__libc_start_main+0xfd)[0x3a1581ed1d]
[CJXU-GEO:32335] [14] /home/x/qdyn/src/qdyn[0x401789]
[CJXU-GEO:32335] *** End of error message ***
/home/x/qdyn/src/qdyn: Aborted

In derivs_all.f90:151 the time derivative of the normal stress needs to be added, i.e.:

    dmain_var = ( dtau_per + dtau_dt - sigma*dmu_dtheta*dth_dt - dtau_dP*dP_dt ) &
                     / ( sigma*dmu_dv + pb%zimpedance )

should become

    dmain_var = ( dtau_per + dtau_dt - sigma*dmu_dtheta*dth_dt - dtau_dP*(dP_dt - dsigma_dt) &
                     / ( sigma*dmu_dv + pb%zimpedance )

with dtau_dP = -mu (rename variable to make more sense...).

This only affects 3D simulations with FEAT_STRESS_COUPL = 1.

Dear Dr. Martijn van de Ende

These days I am attempting to incorporating a "new" friction law into QDYN that takes into account the regularized RSF law, flash heating, and the effects of mechanical work on weakening temperature Tw (see the equation set in the attached figure). The most difficult aspect of this approach was adding the differential equation "dmw/dt = tau * v" in order to calculate the real-time mechanical work mw and then Tw. In my modified version, the treatment of mw is identical to your treatment of the state variable theta. Could you tell me whether this implementation is correct?

I have successfully compiled the modified code, but an error titled occurs when I attempt to run the Single asperity (2D) tetorial code in python. The progression may halt at p.run(). I have tried to debug this in numerous ways, but have always failed. I would appreciate your assistance with that.

I have uploaded the updated filed assiciated with this issue. For your convenience, I have also categorized the positions I have made changes, as shown below:

1. In friction. f90: line18, 41-43, 78-79, 94-104, 113-114, 138-148
2. In problem_class. f90: line132 (neqs 3=>4), 117-118, 129-130, 151
3. In derivs_all. f90: line29, 38, 55, 99, 156
4. In solver. f90: line86, 161, 296
5. In utils. f90: line25, 29, 44, 59, 63, 76
6. In unittests_rsf. f90: line37, 95, 106
7. In unittests. f90: line171, 175. 188
8. In output. f90: line57, 93, 116, 372, 456, 458, 989
9. In input. f90: line58, 158
10. In pyqdyn. py: line114-119, 254, 432

Thanks for your time!

Bowen

modified src(bowen).zip

Originally posted by @bowenyu1 in #79 (comment)

Currently, the Runge-Kutta (RK) solver does not synchronise the error estimates or time steps across processors, which is required for parallel (MPI) simulations. ode_rk.f90 needs to be modified to include an MPI gather at the end of each error estimation step.

When a simulation is run with FINITE=3 (forced symmetry to simulate free surface), the resulting snapshot output does not maintain sorted order. The last mesh element in each snapshot has a position that should correspond with the first element. To illustrate this a bit better, I have plotted the mesh node position in a 1D fault against its index:

However, the stress/slip/velocity/etc. output do not exactly match the first element:

When plotting slip distributions, I get the following result (note discontinuity at the bottom):

Or when the elements are sorted by their position (so that last element becomes the first; note discontinuity at the top):

This problem does not emerge when FINITE=1. The (Python) code to reproduce the issue:

from pyqdyn import qdyn
import numpy as np
import matplotlib.pyplot as plt

p = qdyn()

N = np.power(2, 6)
t_yr = 3600*24*365.0
V = 1e-9

# Python dictionary with general settings
set_dict = {
    "N": N,
    "NXOUT": np.power(2, 1),
    "NTOUT": 100,
    "ACC": 1e-7,
    "MU": 30e9,
    "DTTRY": 1e-8,
    "TMAX": 200*t_yr,
    "MESHDIM": 1,
    "VS": 3000,
    "L": 1.0e3,
    "FINITE": 3,
    "IC": 0,
    "SOLVER": 1,
    "V_PL": V,
}

# Python dictionary with rate-and-state friction parameters
set_dict_RSF = {
    "RNS_LAW": 0,
    "THETA_LAW": 1,
    "DC": 10e-3,
    "V_0": V,
    "V_SS": V,
    "MU_SS": 0.6,
    "A": 1e-2,
    "B": 2e-2,
    "SIGMA": 50e6,
}

set_dict["TH_0"] = 0.99*set_dict_RSF["DC"]/set_dict_RSF["V_0"]
set_dict["SET_DICT_RSF"] = set_dict_RSF

p.settings(set_dict)
p.render_mesh()
p.write_input()
p.run()

I am a novice in the dynamic earthquake simulation. I have read some papers about this respect, but I find that reading and writing code maybe more important than reading papers. It is very lucky that I find the qdyn program, and I am curious about are there any other freely available code about dynamic earthquake simulation. I do not mean that the qdyn program is not good, to the opposite, it gives me a lot of help. I just want to learn more about the dynamic earthquake simulation. Would you like to give me some help or advices?

I suggest adding topics such as earthquakes, seismology, geophysics in the About section at https://github.com/ydluo/qdyn.

From output.f90::ot_write:

open(pb%ot%unit,access='APPEND',status='old',iostat=ios)
if (ios>0) stop 'Fatal error: ot_write: Error opening a fort.18 file'
!JPA add test for the first time we try to open this file but it does not exist yet
!JPA add test to prevent appending data to a file from a previous simulation

These checks should be implemented in output.f90::ot_init; it is not necessary to check every output time step.

At the moment we contribute to the documentation via a Google Doc, which is not an ideal solution in my opinion. We could switch to using Github Pages as a way of documenting our code, in which the documentation is generated from one or more Markdown files contained in the /docs directory (in the same way as our Readme works). I have copied a piece of the current documentation from the PDF to my fork, see this page for the result. It will be a bit of work to fully copy the entire documentation into a Markdown file (my guess: 2-3 hours), but I see the following (long-term) advantages:

1. Documentation is kept alongside the code within the repository, rather than on a Google Drive
2. Individual features/scripts/examples/tests could be documented individually, and linked to the main documentation page. In this way not everything has to go on one single page/document, which allows for a more clearly structured documentation
3. Changes in the documentation are kept as commits, so version history is maintained
4. Related: updates to the documentation (or lack thereof) can be tied to individual pull requests
5. Anyone with a Github account (not just QDYN core developers) can contribute to the documentation and make suggestions via pull requests
6. Markdown (HTML) is a bit more flexible in terms of formatting, I feel

Any thoughts or feelings?

Dear QDYN developers,

During the first SCEC-SEAS benchmark practice we discovered the QDYN slip output is updated with only first order accuracy by outputting v*dt instead of the high-order RK/BS method as used in computation.
It is only a output issue and does not affect the actual computation, yet it better be fixed.

Yingdi

The p.ot giving nan value for every parameter.
Python as well as Matlab Wrapper not giving output with p.ot.
Please help.

Dear co-developers I recently experienced some minor and confusion issues while running the current version of QDYN (R518) with older version of Matlab. The code runs properly but there is no screen output of any kind in Matlab (also no warning / error message). Everything is fine with newer version of Matlab though (Matlab 9.0 +).
Does anyone experience similar issue here? It does not impact the functionality of the code (can check the progress in fort.*) but it's just annoying.

YD

Since 3D simulations are a prominent feature of QDYN, they should be tested on a regular basis. For new users, it would be helpful to have a tutorial as a starting point for more advanced simulations.

To-do:

• Create Jupyter notebook with a simple 3D simulation (including dynamic OX output?)
• Add notebook to documentation (advanced tutorial)
• Generate frozen results for benchmark testing
• Update test suite

I was trying to plot the timeseries in Q-Dyn writing: "qdyn_plot.timeseries(p.ot[0], p.ot_vmax)", but each time I tried I received the following error:

KeyError                                  Traceback (most recent call last)
~/anaconda3/lib/python3.8/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   3079             try:
-> 3080                 return self._engine.get_loc(casted_key)
   3081             except KeyError as err:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'v_max'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
<ipython-input-47-6c9f3387d7df> in <module>
      1 # Spatio-temporal evolution of slip rates
      2 qdyn_plot.slip_profile(p.ox, warm_up=25*t_yr)
----> 3 qdyn_plot.timeseries(p.ot[0], p.ot_vmax)

~/Q-Dyn/qdyn-master/utils/post_processing/plot_functions.py in timeseries(ot, ot_vmax)
     20     axes[1].set_ylabel("state [s]")
     21 
---> 22     axes[2].plot(ot_vmax["t"] / t_year, ot_vmax["v_max"])
     23     axes[2].set_ylabel("max v [m/s]")
     24     axes[2].set_xlabel("time [yr]")

~/anaconda3/lib/python3.8/site-packages/pandas/core/frame.py in __getitem__(self, key)
   3022             if self.columns.nlevels > 1:
   3023                 return self._getitem_multilevel(key)
-> 3024             indexer = self.columns.get_loc(key)
   3025             if is_integer(indexer):
   3026                 indexer = [indexer]

~/anaconda3/lib/python3.8/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   3080                 return self._engine.get_loc(casted_key)
   3081             except KeyError as err:
-> 3082                 raise KeyError(key) from err
   3083 
   3084         if tolerance is not None:

KeyError: 'v_max'

Thank you

The output from ox_dyn (in output.f90::ox_write) is outdated. The time of the snapshot is incorrect, the values of dtau_dt are not initialised.

• Identify the required output quantities (should be the same as for the regular ox output
• Create one universal file writer for ox and ox_dyn output
• Create global ox_dyn output unit (currently a magic number 20001+x)
• Fix non-initialised output

I can't make on MacOS Mojave 10.14. I've tried all compilers(gfortran, ifort, mpif90). I'm new fortran developer. What are appropriate compiler flags?
parallel.f90:219:21:

219 | call MPI_ALLREDUCE(send, buffer, countval, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ier)
| 1
......
242 | call MPI_ALLREDUCE(sendbuf,recvbuf,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,ier)

In release 3.0.0, output_ox_last is only written at the timestep N where pb%it == pb%itstop (N=-1). It'd be useful to set N = -1 by default, which signals that a snapshot should be written at pb%itstop, but to have the option to write the output every N steps for N > 0 (including one at pb%itstop). The latter situation is useful when the simulation is running (i) through the wrapper and then it's interrupted; (ii) through a jobscript in the cluster and the job reaches the time limit before reaching TMAX.

What steps will reproduce the problem?
1.  Extract QDYN_R179.zip
2.  Addpath of new folder to Matlab
3.  Try to run any code in this folder 

What is the expected output? What do you see instead?
The code should compile.  I get the following error:

Lb =

   1.2094e+04


filename =

JP_2D_s2_dc0.5twm3000ts1000L2066.7827nx1W206.6783nw207dip14.mat


Lb_over_dx =

    0.0059

MESHDIM = 2
'~' is not recognized as an internal or external command, 
operable program or batch file. 
Error using textread (line 166)
File not found.

Error in qdyn>read_qdyn_out (line 469)
[ot.t, ot.locl, ot.cl, ot.p, ot.pdot, ...

Error in qdyn (line 352)
    [ot,ox]= read_qdyn_out(NAME);

Error in JP_2d_s (line 95)
[p,ot0,ox0]  = qdyn('run',p);

What version of the product are you using? On what operating system?
QDYN_r179 / Win7 /Matlab 2013a

Please provide any additional information below.

I can get code to compile when I have created the code in the 
TohokuEarthquakeExample.zip download from the website.  The differences that I 
notice between the two folders some differences:

a) the tohoku earthquake folder has a qdyne.exe file not located in the 
QDYN_R179.zip download.

b)the QDYN_R179.zip download has more .f90 (fortran extension files; e.g 
friction.f90 , main.f90, okada.f90).  Are these supposed to compiled into a 
qdyn.exe executable based on the inputs specified from matlab? 

Thank you kindly,

Paul

Original issue reported on code.google.com by [email protected] on 16 Apr 2013 at 9:25

In the single/double asperity example notebooks, L is defined as Length of fault / nucleation length. This is not correct, and should read Length of fault / asperity length. This should also be updated in the tutorial documentation.

The default values in pyqdyn.py for W and DIP_W are inconsistent with those in qdyn.m and the documentation. Set defaults to W = 50e3 and DIP_W = 90.

Introduction

Regardless of the time frame over which the HDF binary output is scheduled to be implemented, I think it would be useful at some point to explore this (as discussed in #42). Particularly for large-scale (3D) simulations, the data volume can impose some limitations on the post-processing (memory, disk space, IO time). Luckily for us, HDF has a design principle that is similar to an OS file system, so that we can exploit particular structures in the data.

Proposed data structure

With reference to the figure below, it is possible to store the various data types (time series, snapshots) in groups and subgroups, starting with the "root" group, and metadata can be assigned to each each (sub)group header. The data contained in a given group can be shared with other groups, which enables us to share/recycle invariant data (similar to symbolic links in file systems, or relational databases in SQL environments). Lastly, instead of garden-variety data (integers, floats), files (e.g. the qdyn.in file) can also be directly attached to a group.

Concretely, for snapshot output, the OX mesh locations can be stored once and shared with each snapshot (1, 2, ..., N). The time of each snapshot is simply a meta tag, which eliminates 4 out of 10 OX quantities to be stored for each snapshot. An identical scheme applies to dynamic OX output (sampled at a different rate). Similarly, for the time series output the time vector (but also v_max, tau_max, etc.; not indicated in the figure) is shared between OT and IOT, each sampled at a different location (1, 2, ... K). OT and IOT output are consequently generalised to a common data structure.

For reading the simulation output, one would no longer need to read the entire output file into memory. Instead, one selects only the relevant data columns of a given group. So if I require the slip of snapshots from t = 10 to t = 25, I only need to check the meta data of the snapshot groups, and select the slip column for each of those that satisfy the selection criteria. The mesh locations need to be extracted only once. In addition, I believe it is also possible to slice columns: if I need time series v_max data from t = 100 to t = 500, I could select only a portion of the v_max column.

Implementation challenges

The data structure proposed above (column-based) is very different from the current output structure (row-based). While this is not necessarily a problem, it requires a bit more thought to implement. The procedure for creating and modifying an HDF5 file in Fortran is a little convoluted (see e.g. this SO answer for a "simple" implementation example). While playing around, I already lost 2 hours just trying to get a tutorial code compiled, so I expect that it will take some time so set-up and debug everything.

Even though HDF5 overlays OpenMPI IO, I am not sure whether the data structure proposed above is suitable for parallel IO. Since IO likely does not incur a large overhead for most simulations, it might be better to stick with a conventional MPI gather and do the IO in serial mode.

My biggest fear is that the data will likely be corrupted if the HDF file (and its subspaces) is not closed properly. So when a simulation crashes or is manually terminated, all of the simulation data may be lost (which does not happen with ASCII file formats). Since we all work very hard to eliminate all bugs and instabilities, I'm not so much worried about crashing simulations, but it happens often that I terminate a simulation at an early stage (e.g. after the post-seismic phase) to inspect the data. In C++ and Python it is possible to intercept a keyboard interrupt and terminate things safely, but in Fortran this seems to be a bit tricky (from quickly browsing through this PDF). We'll have to see in practice how we could best implement a deconstructor triggered by various exceptions.