QDYN release 3.0.0

[martijnende]

What's new in version 3.0.0?

This is a major upgrade of the QDYN code base, with (some) compatibility-breaking changes made in both the Fortran code and the Python wrapper. This version sees the release of some interesting new features that greatly extend the range of use cases of QDYN.

Improved 3D mesh builder (#71)

Up to now, the 3D fault mesh was generated by mesh.f90 based on the fault dip and element size provided in the input script. To allow for more flexibility in the design of the simulation, this functionality has been moved inside of pyqdyn.py. After initiating the mesh (qdyn.render_mesh()), a more advanced 3D mesh can be rendered with qdyn.compute_mesh_coords(). See the docstring of this function for its usage. New example notebooks demonstrating correct usage have been added (3D_mesh_builder.ipynb and intersecting_faults.ipynb)

Implementation of viscosity (#72)

Previously, simulations with rate-and-state friction (RSF) solved for the shear stress as a function of slip rate, whereas simulations with the Chen-Niemeijer-Spiers simulation solved for the slip rate as a function of the shear stress. From this version onward, all simulations solve for the slip rate regardless of the underlying friction law. This permits viscous creep to operate in parallel to frictional sliding governed by RSF, and time-dependent shear stress dissipation by creep. The role of viscous creep has been investigated by Beall et al. (2022; GRL), which is a suggested reference for citing the QDYN project in combination with creep. Set INV_VISC > 0 to specify the viscosity for each fault mesh element. Note that this parameter has dimensions of m/(Pa.s), such that V = INV_VISC * TAU for purely viscous creep.

As a consequence, the user must provide the initial state of stress (TAU) instead of the initial velocity (V_0). For backward-compatibility, the Python wrapper will automatically detect if TAU has been properly set; if not, it will attempt to convert V_0 into TAU using the friction law and the provided viscosity values (INV_VISC).

Simulation restarts (#83)

The idea of this implementation is to simplify the workflow to restart a simulation without having to starting at 0s. When restarting a model, the final state of the previous simulation (written to output_ox_last) is read and the simulation is continued from that point onward. All the outputs are appended to the files of the previous simulation instead of being rewritten, which allows the user to have a single folder with all the model runs.

A simulation restart is enabled by setting qdyn.run(restart=True) using the wrapper, or by calling qdyn restart when an input file has already been made. An example notebook (restart_simulation.ipynb) demonstrates the correct use of this feature.

Fault labels (#83)

Fault labels are introduced to facilitate postprocessing when working with multiple faults. When setting-up the fault geometry, each fault element can be assigned an integer representing the fault it belongs to. These fault labels appear in the simulation output, so that the output data can be easily separated.

Potency calculation (#83)

Up to now, QDYN computed a macroscopic value of potency and potency rate using all the nodes of the mesh. While this value holds physical significance when working with one fault, this is no longer the case when working with multiple faults. Using the fault labels, the potency (rate) is now computed for individual faults.

Logging and debugging

As part of the revisions to the I/O modules, a separate logger module has been added to relay information either to a screen or a log file. By default, all messages broadcasted during the simulation (e.g. initialisation messages or the current state of the simulation) are now directed to a log file. To receive log messages on the screen, set VERBOSE = 1.

By setting DEBUG = 1, many function calls are directly written to the log file each time step (and for each processing node). This can be useful information to identify exactly when a simulation crashes. Future developments will include a crash dump of the last N simulation steps when NaNs are encountered.

Test suite updated and moved to GitHub Actions (#74)

Since TravisCI no longer offers free CI/CD services to open source projects, the test suite was ported to GitHub Actions.

The tests involving stick-slip motion are now evaluated using a correlation coefficient rather than with a point-by-point comparison, which is highly sensitive to small time shifts.

Added global path to pyqdyn.py (#75)

To facilitate the import of pyqdyn.py functions, a qdyn package was created. This can be added to the global Python path by installing through pip:

cd [qdyn root directory]
pip install -e .

The QDYN wrapper is then imported with:

from qdyn import qdyn

This is an alternative to adding the path to pyqdyn.py each time (which still works).

⚠ Note: the qdyn/src directory was renamed to qdyn/qdyn.

Dropped MATLAB support (#76)

Due to a lack of development and testing capabilities, support for MATLAB was dropped. The MATLAB wrapper (qdyn.m) was outdated and no longer compatible with the QDYN codebase. Simultaneously, all of the .m examples were removed. MATLAB users can revert to a previous release to continue to use the MATLAB wrapper and examples.

Bug fixes since version 2.2.0

• Fixed sign bugs in the stress kernel that prevented the solver from converging for 3D normal fault simulations, and that caused a normal fault to slide in the wrong direction
• Correctly included the partial derivative of shear to normal stress in derivs_all.f90 (which now accounts for normal stress changes in 3D simulations) (Normal stress coupling needs to be added to partial diffs #73)
• Fixed an issue within the 3D FFT stress kernel initialisation (fault_stress.f90:init_kernel_3D_fft), whereby it was implicitly assumed that the fault was contiguous everywhere (which is not a requirement in the along-dip direction)
• When running parallel simulations, some arrays where subject to corruption because the array iot_buf was being manipulated by all the processors, instead of by the master processor only.
• Fixed NPROCS typo in the documentation (NPROCS in documentation should be NPROC #66)
• Fixed documentation error in single_asperity.ipynb example (Number of mesh elements #81, Typo in asperity notebooks #82)

Minor updates

• Added a version number to pyqdyn.py
• Added comments on sign conventions in fault_stress.f90
• The frequency of output generation is now set independently for the logger, time-series, and snapshots with NTOUT_LOG, NTOUT_OT, and NTOUT_OX, respectively.

[martijnende]

TODO

• Increase final snapshot time precision (The restarted solution differs from the continuous solution. #87)
• Check this comment and verify MPI execution (mpi and serial are not same #85)
• Check and merge Constanza's PR (Output velocity time-series for each individual fault #89)
• Include requirements.txt

[martijnende]

• fix repo URL (ydluo -> ydluo-c) for documentation