elastic.py

Go to the documentation of this file.

# This file is part of the ExaHyPE2 project. For conditions of distribution and
# use, please see the copyright notice at www.peano-framework.org
import sympy
 
"""
We import the required modules for this project.
We always need the 'peano4' module as this is our core project.
Since we are creating a SymHyPE/ExaHyPE 2 application, we additionally
need to import the 'exahype2' and 'symhype' modules.
"""
import peano4
import exahype2
import symhype
 
"""
We specify the space dimensions here.
We support either 2- or 3-dimensional problems.
"""
dimensions = 2
 
"""
The number of finite volumes per axis in one patch.
"""
patch_size = 16
 
"""
The number of levels the mesh is refined.
"""
depth = 5
 
"""
The simulation end time.
"""
end_time = 1.0
 
"""
Choose domain size and offset.
"""
size = [20.0, 20.0, 20.0]
offset = [-5.0, 0.0, 0.0]
 
"""
Choose how often a snapshot is written.
"""
time_in_between_two_snapshots = end_time / 10
# time_in_between_two_snapshots = 0 # Set to 0 to disable I/O
 
"""
Switch between 'Release', 'Debug', 'Asserts', 'Trace', 'Stats'.
"""
compile_mode = "Release"
 
"""
We first create a new ExaHyPE 2 project.
For this, we specify the (nested) namespaces, the name of our main file and our executable name.
"""
my_project = exahype2.Project(
    namespace=["tutorials", "symhype", "elastic"],
    project_name="Elastic",
    directory=".",
    executable="Elastic",
)
 
"""
Add the Finite Volumes solver using named arguments.
This is the way you can add further PDE terms.
This requires the 'BlockStructured' toolbox and 'ExaHyPE' to be built.
"""
my_solver = exahype2.solvers.fv.godunov.GlobalAdaptiveTimeStep(
    name="FVSolver",
    patch_size=patch_size,
    unknowns=dimensions
    + (
        3 if dimensions == 2 else 6
    ),  # 3 (v_x, v_y, v_z) + 6 (σ_xx, σ_yy, σ_zz, σ_xy, σ_xz, σ_yz)
    auxiliary_variables=0,
    min_volume_h=(1.1 * min(size[0:dimensions]) / (3.0**depth)),
    max_volume_h=(1.1 * min(size[0:dimensions]) / (3.0**depth)),
    time_step_relaxation=0.5,
)
 
"""
We want to define our PDE symbolically.
"""
my_pde = symhype.FirstOrderConservativePDEFormulation(
    unknowns=dimensions + (3 if dimensions == 2 else 6),
    auxiliary_variables=0,
    dimensions=dimensions,
)
 
# Define the unknowns.
# v: velocity components [v_x, v_y, v_z]
v = my_pde.name_Q_entries(0, dimensions, "v")
# sigma: stress components [σ_xx, σ_yy, σ_zz, σ_xy, σ_xz, σ_yz]
sigma = my_pde.name_Q_entries(dimensions, 3 if dimensions == 2 else 6, "sigma")
 
# Define material parameters.
rho = sympy.symbols("rho")
cp = sympy.symbols("cp")
cs = sympy.symbols("cs")
 
# Lamé parameters.
mu = rho * cs * cs
lamb = rho * cp * cp - 2.0 * mu
 
if dimensions == 3:
    """
    Define the fluxes in each spatial direction.
    The overall conservative form is:
      ∂ₜQ + ∂ₓ(F_x(Q)) + ∂ᵧ(F_y(Q)) + ∂𝓏(F_z(Q)) = S
    where Q = [v_x, v_y, v_z, σ_xx, σ_yy, σ_zz, σ_xy, σ_xz, σ_yz]^T.
    """
    # Flux [unknowns, dimensions]
    # --- x-direction (direction index 0) ---
    # Momentum equations:
    my_pde.F[0, 0] = -1.0 / rho * sigma[0]  # v_x: -σ_xx/ρ
    my_pde.F[1, 0] = -1.0 / rho * sigma[3]  # v_y: -σ_xy/ρ
    my_pde.F[2, 0] = -1.0 / rho * sigma[4]  # v_z: -σ_xz/ρ
 
    # Stress equations:
    my_pde.F[3, 0] = -(lamb + 2.0 * mu) * v[0]  # σ_xx: -(λ+2μ)v_x
    my_pde.F[4, 0] = -lamb * v[0]  # σ_yy: -λ v_x
    my_pde.F[5, 0] = -lamb * v[0]  # σ_zz: -λ v_x
    my_pde.F[6, 0] = -mu * v[1]  # σ_xy: -μ v_y
    my_pde.F[7, 0] = -mu * v[2]  # σ_xz: -μ v_z
    my_pde.F[8, 0] = 0  # σ_yz: 0
 
    # --- y-direction (direction index 1) ---
    # Momentum equations:
    my_pde.F[0, 1] = -1.0 / rho * sigma[3]  # v_x: -σ_xy/ρ
    my_pde.F[1, 1] = -1.0 / rho * sigma[1]  # v_y: -σ_yy/ρ
    my_pde.F[2, 1] = -1.0 / rho * sigma[5]  # v_z: -σ_yz/ρ
 
    # Stress equations:
    my_pde.F[3, 1] = -lamb * v[1]  # σ_xx: -λ v_y
    my_pde.F[4, 1] = -(lamb + 2.0 * mu) * v[1]  # σ_yy: -(λ+2μ)v_y
    my_pde.F[5, 1] = -lamb * v[1]  # σ_zz: -λ v_y
    my_pde.F[6, 1] = -mu * v[0]  # σ_xy: -μ v_x
    my_pde.F[7, 1] = 0  # σ_xz: 0
    my_pde.F[8, 1] = -mu * v[2]  # σ_yz: -μ v_z
 
    # --- z-direction (direction index 2) ---
    # Momentum equations:
    my_pde.F[0, 2] = -1.0 / rho * sigma[4]  # v_x: -σ_xz/ρ
    my_pde.F[1, 2] = -1.0 / rho * sigma[5]  # v_y: -σ_yz/ρ
    my_pde.F[2, 2] = -1.0 / rho * sigma[2]  # v_z: -σ_zz/ρ
 
    # Stress equations:
    my_pde.F[3, 2] = -lamb * v[2]  # σ_xx: -λ v_z
    my_pde.F[4, 2] = -lamb * v[2]  # σ_yy: -λ v_z
    my_pde.F[5, 2] = -(lamb + 2.0 * mu) * v[2]  # σ_zz: -(λ+2μ)v_z
    my_pde.F[6, 2] = 0  # σ_xy: 0
    my_pde.F[7, 2] = -mu * v[0]  # σ_xz: -μ v_x
    my_pde.F[8, 2] = -mu * v[1]  # σ_yz: -μ v_y
 
    """
    Define the eigenvalues.
    For a 3D elastic system the characteristic speeds in a given
    direction are:
      { cp, cs, cs, 0, 0, 0, -cs, -cs, -cp }
    where:
      cp = sqrt((λ+2μ)/ρ)   is the P-wave (compressional) speed,
      cs = sqrt(μ/ρ)        is the S-wave (shear) speed.
    """
    # Eigenvalues [unknowns, dimensions]
    my_pde.eigenvalues[0] = cp
    my_pde.eigenvalues[1] = cs
    my_pde.eigenvalues[2] = cs
    my_pde.eigenvalues[3] = 0
    my_pde.eigenvalues[4] = 0
    my_pde.eigenvalues[5] = 0
    my_pde.eigenvalues[6] = -cs
    my_pde.eigenvalues[7] = -cs
    my_pde.eigenvalues[8] = -cp
else:
    # Flux [unknowns, dimensions]
    my_pde.F[0, 0] = -1.0 / rho * sigma[0]
    my_pde.F[1, 0] = -1.0 / rho * sigma[2]
    my_pde.F[2, 0] = -(2.0 * mu + lamb) * v[0]
    my_pde.F[3, 0] = -lamb * v[0]
    my_pde.F[4, 0] = -mu * v[1]
 
    my_pde.F[0, 1] = -1.0 / rho * sigma[2]
    my_pde.F[1, 1] = -1.0 / rho * sigma[1]
    my_pde.F[2, 1] = -lamb * v[1]
    my_pde.F[3, 1] = -(2.0 * mu + lamb) * v[1]
    my_pde.F[4, 1] = -mu * v[0]
 
    # Eigenvalues [unknowns, dimensions]
    my_pde.eigenvalues[0] = cp
    my_pde.eigenvalues[1] = cs
    my_pde.eigenvalues[2] = 0
    my_pde.eigenvalues[3] = -cs
    my_pde.eigenvalues[4] = -cp
 
"""
Since 'my_pde' only holds the PDE without initial- or boundary conditions,
we still need to properly define initial- and boundary conditions.
This gives us then a complete description of a 'scenario'.
"""
for i in range(dimensions + (3 if dimensions == 2 else 6)):
    my_pde.initial_values[i] = 0
    my_pde.boundary_values[i] = 0
 
"""
Define a source term.
Here we add a time-dependent point force applied to one of the
shear stress components (σ_xy in this case).
"""
t0 = sympy.symbols("t0")
M0 = sympy.symbols("M0")
t = sympy.symbols("t")
force = M0 * t / (t0 * t0) * sympy.exp(-t / t0)
 
max_h = sympy.symbols("MaxAdmissibleVolumeH")
 
if dimensions == 3:
    point_source = sympy.sqrt(
        (10 - my_pde.x[0]) ** 2 + (10 - my_pde.x[1]) ** 2 + (10 - my_pde.x[2]) ** 2
    )
 
    # Apply the force to σ_xy, which is the 7th component in the overall Q (index 6).
    my_pde.sources[0] = 0
    my_pde.sources[1] = 0
    my_pde.sources[2] = 0
    my_pde.sources[3] = 0
    my_pde.sources[4] = 0
    my_pde.sources[5] = 0
    my_pde.sources[6] = sympy.Piecewise((force, point_source <= max_h), (0, True))
    my_pde.sources[7] = 0
    my_pde.sources[8] = 0
else:
    point_source = sympy.sqrt((10 - my_pde.x[0]) ** 2 + (10 - my_pde.x[1]) ** 2)
 
    my_pde.sources[0] = 0
    my_pde.sources[1] = 0
    my_pde.sources[2] = 0
    my_pde.sources[3] = 0
    my_pde.sources[4] = sympy.Piecewise((force, point_source <= max_h), (0, True))
 
# Substitute numerical values for the material and source parameters.
my_pde.substitute_expression(rho, 2.7)
my_pde.substitute_expression(cp, 6.0)
my_pde.substitute_expression(cs, 3.464)
my_pde.substitute_expression(t0, 0.1)
my_pde.substitute_expression(M0, 1000.0)
 
"""
Specify which implementation our solvers uses.
Here we want to set the implementation we get from our symbolically defined PDE,
i.e., we get the C++ implementation which is generated by SymHyPE.
"""
my_solver.set_implementation(
    initial_conditions=my_pde.implementation_of_initial_conditions(),
    boundary_conditions=my_pde.implementation_of_boundary_conditions(),
    flux=my_pde.implementation_of_flux(),
    max_eigenvalue=my_pde.implementation_of_max_eigenvalue(),
    source_term=my_pde.implementation_of_sources(),
)
 
"""
To see which variables (unknowns + auxiliary variables) we can visualise,
let's add a plot description for the variables to our solver.
"""
my_solver.plot_description = my_pde.unknown_identifier_for_plotter()
 
"""
Add the solver to our project
"""
my_project.add_solver(my_solver)
 
"""
Configure some global parameters
"""
my_project.set_global_simulation_parameters(
    dimensions=dimensions,
    size=size[0:dimensions],
    offset=offset[0:dimensions],
    min_end_time=end_time,
    max_end_time=end_time,
    first_plot_time_stamp=0.0,
    time_in_between_plots=time_in_between_two_snapshots,
    periodic_BC=[False, False, False],
)
 
"""
This defines where the output files should go.
If you omit this, output files are automatically put into the application's folder.
"""
my_project.set_output_path("solutions")
 
"""
Configure load balancer for parallel execution.
"""
my_project.set_load_balancer("new ::exahype2::LoadBalancingConfiguration")
 
"""
We need to set the location of our core libraries ('Peano4').
This helps us to resolve any dependencies.
Additionally, we specify the build mode which you can also change to a different mode.
"""
my_project.set_Peano4_installation(
    "../../../", mode=peano4.output.string_to_mode(compile_mode)
)
 
"""
We generate and grab the underlying core project of 'Peano4'.
This gives us access to some functions we want to use to finalise and build this project.
"""
my_project = my_project.generate_Peano4_project(verbose=False)
 
"""
Finally, we want to build our project.
First, all of the necessary glue code is generated in the application folder,
then 'make' is invoked automatically which compiles the generated code and links against our core libraries
and toolboxes which have been built before.
You can also always invoke 'make' yourself to compile, or cleanup with 'make clean'.
"""
my_project.build(make=True, make_clean_first=True, throw_away_data_after_build=True)
 
print("\nPDE is: ")
print(my_pde.__str__())
print(my_solver)
print(my_project)