This repository contains everything needed to run the simulation examples found in the MPSA-Newmark paper.
The model classes are standardized setups for solving the elastic wave equation with absorbing boundary conditions with PorePy. There are currently two model class setups:
- elastic_wave_equation_abc considers a general setup for solving the elastic wave equation with absorbing boundaries on all domain sides.
- elastic_wave_equation_abc_linear considers the case where all equations related to fractures in the domain are discarded. In practice this means that the model class is used for solving elastic wave propagation in rocks without fractures, or in rocks with open fractures (inner zero traction Neumann boundaries). In both cases the problem is linear, meaning that the Jacobian doesn't change between time steps. Hence, this model setup facilitates that the Jacobian is assembled only once. The linear model setup should be used together with the custom run-model function which is detailed in the section below.
solvers contains mixins for custom solvers. Specifically, the mixin that is there now will allow for PETSc usage whenever that is available. It also takes into consideration whether the Jacobian should be assembled or fetched, where the latter is the case if elastic_wave_equation_abc_linear is the model setup to be run.
run_models contains custom run-models functions. In the case of a the Jacobian being assembled only once, adaptations were needed to the run-model-function which originally lies within PorePy such that the residual was not assembled twice per time step.
The convergence analyses presented in the article are performed with homogeneous Dirichlet conditions on a 3D simplex grid:
- Convergence in space and time:
- Convergence in space:
- Convergence in time:
All the runscripts utilize manufactured_solution_dynamic_3D as the manufactured solution setup.
Convergence of the solution is performed in a quasi-1D setting
Convergence in space and time:
All the runscripts utilize model_convergence_ABC2 as the model class setup.
The energy decay analysis is performed both for successive refinement of the grid, as well as for varying wave incidence angles.
-
Successive grid refinement is done by running the script render_figure_energy_decay_space_refinement. The script which renders the figure uses the runscript runscript_ABC_energy_vary_dx for running the simulations.
-
Varying the wave incidence angle,
$\theta$ , is done by running the script runscript_ABC_energy_vary_theta. -
A quasi-1d version is found in runscript_ABC_energy_quasi_1d.
Simulation example runscripts are found within this folder.
- The simulation from Example 1.1, which considers a seismic source located inside an inner transversely isotropic domain, is run by runscript_example_1_1_source_in_inner_domain.
- The simulation from Example 1.2, which considers a seismic source located outside an inner transversely isotropic domain, is run by runscript_example_1_2_source_in_outer_domain.
- The simulation from Example 2, which considers a layered heterogeneous medium with an open fracture, is run by runscript_example_2_heterogeneous_fractured_domain.
A collection of utility material is found within the utils directory:
- anisotropy mixins contains mixins for anisotropic stiffness tensors.
- perturbed_geometry_mixins contains mixins for three types/configurations of perturbed geometry.
- stiffness tensors contains a fourth order stiffness tensor object for a transversely isotropic material.
- utility functions contains mostly functions related to analytical solution expressions and fetching subdomain-related quantities.
I refer to the files within the directory for more details about the specific contents.
Tests are covering:
- MPSA-Newmark convergence with homogeneous Dirichlet conditions in 2D and 3D.
- MPSA-Newmark convergence with absorbing boundary conditions.
- Construction of the transversely isotropic tensor.
- The utility function
inner_domain_cells
which is used in the construction of the transversely isotropic tensor.