Mondaic
This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.
- Reference solution:
analytic - Physis: Diffusion equation
In this notebook we benchmark the numerical solution of the diffusion equation in a simple test case for which an analytic solution exists. Specifically, we are considering the fundamental solution of the diffusion equation in homogeneous media, which is often called diffusion kernel.
To this end, we consider the diffusion equation with a Dirac measure as initial values and a constant diffusion tensor with :
On , (), this equation has a unique solution given by:
which can also be interpreted as a Wiener process.
Remarks:
- Representing the delta source on the finite element mesh requires a proper projection onto the finite element basis. Instead, we simulate the time interval from
0.01 sto0.02 s. - While the analytic solution is defined on the unbounded domain , we restrict the computational to a box or sphere. To avoid artifacts from the artificial boundaries, we use a fairly short simulation time.
%matplotlib inline
%config Completer.use_jedi = False
import os
import numpy as np
import salvus.namespace as sn# Number of processes SalvusCompute will run with.
# Get it from the environment or default to 4.
MPI_RANKS = int(os.environ.get("NUM_MPI_RANKS", 4))
# Choose on which site to run this.
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")domain_size = 2.0
k = 0.5
FINAL_TIME = 0.01def wiener_process(time, dim, mean, points, k):
x_bar = points - mean
u_exact = (
1
/ (np.sqrt((k * np.pi * time) ** dim))
* np.exp(-np.einsum("ijk,ijk->ij", x_bar, x_bar) / (4 * k * time))
)
return u_exactdef create_mesh(dim, mode, tensor_order, domain_size, k):
if dim == 2 and mode == "Cartesian":
m = sn.simple_mesh.CartesianHomogeneousAcoustic2D(
vp=1.0,
rho=1.0,
x_max=domain_size,
y_max=domain_size,
max_frequency=30.0,
)
m.advanced.tensor_order = tensor_order
mesh = m.create_mesh()
elif dim == 2 and mode == "spherical":
m = sn.simple_mesh.basic_mesh.SphericalHomogeneousAcoustic2D(
radius=domain_size / 2, vp=2.0, rho=1.0, max_frequency=20.0
)
m.advanced.tensor_order = tensor_order
mesh = m.create_mesh()
mesh.points[:, 0] += 1.0
mesh.points[:, 1] += 1.0
elif dim == 3 and mode == "Cartesian":
m = sn.simple_mesh.CartesianHomogeneousAcoustic3D(
vp=1.0,
rho=1.0,
x_max=domain_size,
y_max=domain_size,
z_max=domain_size,
max_frequency=10.0,
)
m.advanced.tensor_order = tensor_order
mesh = m.create_mesh()
f = np.ones_like(mesh.elemental_fields["VP"])
del mesh.elemental_fields["VP"]
del mesh.elemental_fields["RHO"]
mesh.attach_field("M0", 1.0 * f)
mesh.attach_field("M1", k * f)
mesh.attach_field(
"uinit",
wiener_process(
0.01,
mesh.ndim,
[1.0] * mesh.ndim,
mesh.points[mesh.connectivity],
k,
),
)
return meshdef simulate(mesh, final_time):
mesh.write_h5("init.h5")
w = sn.simple_config.simulation.Diffusion(mesh=mesh)
w.domain.polynomial_order = mesh.shape_order
w.physics.diffusion_equation.initial_values.filename = "init.h5"
w.physics.diffusion_equation.initial_values.format = "hdf5"
w.physics.diffusion_equation.initial_values.field = "uinit"
w.physics.diffusion_equation.final_values.filename = "final.h5"
w.physics.diffusion_equation.end_time_in_seconds = final_time
w.validate()
sn.api.run(
input_file=w,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="diffusion",
overwrite=True,
ranks=MPI_RANKS,
)def analysis(mesh, final_time):
tensor_order = mesh.shape_order
if tensor_order == 1:
TOL = 1e-2
elif tensor_order == 2:
TOL = 1e-3
else:
TOL = 1e-4
solution = sn.UnstructuredMesh.from_h5("diffusion/final.h5")
solution.attach_field("init", mesh.elemental_fields["uinit"])
u_analytic = wiener_process(
0.01 + final_time,
mesh.ndim,
[1.0] * mesh.ndim,
solution.points[solution.connectivity],
k,
)
u_salvus = solution.elemental_fields["uinit"]
residuals = u_salvus - u_analytic
u_max = u_salvus.max()
ref_max = u_analytic.max()
solution.attach_field("analytic_solution", u_analytic)
solution.attach_field("solution", u_salvus)
del solution.elemental_fields["uinit"]
solution.attach_field("residuals", residuals)
print(
f"simulation time: {final_time}\n",
f"|| u ||_inf = {u_max}\n",
f"|| u_exact ||_inf = {ref_max}\n",
f"|| u - u_exact ||_inf = {np.abs(residuals).max()}\n",
f"|| u - u_exact ||_inf / || u_exact ||_inf = {np.abs(residuals).max() / ref_max}\n",
)
np.testing.assert_allclose(
u_salvus, u_analytic, rtol=1e-6, atol=ref_max * TOL
)
return solutionmesh = create_mesh(
dim=2, mode="Cartesian", tensor_order=1, domain_size=2.0, k=0.5
)
simulate(mesh, 0.01)
solution = analysis(mesh, 0.01)
solutionSalvusJob `job_2403151021078004_f718bfac8d` running on `local` with 4 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32 -> Current Task: Time loop complete* Downloaded 1.6 MB of results to `diffusion`. * Total run time: 0.71 seconds. * Pure simulation time: 0.53 seconds. simulation time: 0.01 || u ||_inf = 31.959190368652344 || u_exact ||_inf = 31.830988618379067 || u - u_exact ||_inf = 0.12820175027327707 || u - u_exact ||_inf / || u_exact ||_inf = 0.004027576768358805
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb420ec2fd0>
mesh = create_mesh(
dim=2, mode="Cartesian", tensor_order=2, domain_size=2.0, k=0.5
)
simulate(mesh, 0.01)
solution = analysis(mesh, 0.01)
solutionSalvusJob `job_2403151021084766_448fcef43a` running on `local` with 4 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32
* Downloaded 4.3 MB of results to `diffusion`. * Total run time: 0.91 seconds. * Pure simulation time: 0.75 seconds. simulation time: 0.01 || u ||_inf = 31.849172592163086 || u_exact ||_inf = 31.830988618379067 || u - u_exact ||_inf = 0.018218422568931913 || u - u_exact ||_inf / || u_exact ||_inf = 0.0005723486250255098
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb420ea3ed0>
mesh = create_mesh(
dim=2, mode="Cartesian", tensor_order=4, domain_size=2.0, k=0.5
)
simulate(mesh, 0.01)
solution = analysis(mesh, 0.01)
solutionSalvusJob `job_2403151021426846_eff28124e8` running on `local` with 4 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32
* Downloaded 14.0 MB of results to `diffusion`. * Total run time: 5.57 seconds. * Pure simulation time: 4.86 seconds. simulation time: 0.01 || u ||_inf = 31.833200454711914 || u_exact ||_inf = 31.830988618379067 || u - u_exact ||_inf = 0.0022120523721120833 || u - u_exact ||_inf / || u_exact ||_inf = 6.949367481583196e-05
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb455612910>
mesh = create_mesh(
dim=2, mode="spherical", tensor_order=4, domain_size=2.0, k=0.5
)
simulate(mesh, 0.01)
solution = analysis(mesh, 0.01)
solutionSalvusJob `job_2403151021901275_eae7cf9eb7` running on `local` with 4 rank(s). Site information: * Salvus version: 0.12.16 * Floating point size: 32
* Downloaded 2.5 MB of results to `diffusion`. * Total run time: 1.95 seconds. * Pure simulation time: 1.36 seconds. simulation time: 0.01 || u ||_inf = 31.83287811279297 || u_exact ||_inf = 31.830988618379067 || u - u_exact ||_inf = 0.0018894944139020708 || u - u_exact ||_inf / || u_exact ||_inf = 5.936021769713698e-05
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fb45598d710>
PAGE CONTENTS