The restarted solution differs from the continuous solution.

[eyupsopaci]

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.

[eyupsopaci]

The inconsistency arises due to the output module. The digit precision for the output is inadequate in output.f90 Line 493.

  allocate(pb%ox%fmt(pb%ox%nox))
  ! Default output format
![check_restart_ooo](https://github.com/ydluo/qdyn/assets/65046328/a73353f5-3b19-4f94-ad29-733a60a8f8f8)

  pb%ox%fmt = "(e15.7)"

When it is set to 'e24.14', the solutions perfectly overlap (restart1 is with the longer digits).

[martijnende]

Thanks Eyup, I'll include this in the final release (#84)

[eyupsopaci]

Thanks Eyup, I'll include this in the final release (#84)

I changed the output the precision of the ox file in lines 493-506 since each variable requires different precision:

  allocate(pb%ox%fmt(pb%ox%nox))
  ! Default output format
  pb%ox%fmt = "(e15.7)"
  ! Simulation step is an integer
  pb%ox%fmt(1) = "(i15)"
  ! Time needs higher precision
  pb%ox%fmt(2) = "(e24.14)"
  ! Theta may require even higher precision
  pb%ox%fmt(7) = "(e26.16)"
  ! tau and sigma
  pb%ox%fmt(8) = "(e20.12)"
  pb%ox%fmt(11) = "(e20.12)"
  ! Fault label is an integer
  pb%ox%fmt(pb%ox%nox)= "(i15)" 

The above solution sustains a precise restart for a stiff setup, where theta values can exceed 1e9. The following figure shows the comparison between the seamless and restarted simulations.