Version:

This documentation is not for the latest stable Salvus version.

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.

# Benchmark: Diffusion equation

• 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 $\mathcal{D}(\mathbf{x}) = \left(\begin{array}{cc}k & 0\\0 & k \end{array}\right)$ with $k > 0$:

\begin{aligned} \partial_t u(\mathbf{x},t) - \nabla \cdot \mathcal{D}(\mathbf{x})\, \nabla u(\mathbf{x},t) &= 0,\\ u(\mathbf{x},0) &= \delta(\mathbf{x}-\mathbf{\bar{x}}). \end{aligned}

On $\mathbf{R}^d$, ($d=2,3$), this equation has a unique solution given by:

$u_{\text{exact}} = \frac{1}{(k\pi\,t)^{d/2}} \exp\left( - \frac{(\mathbf{x}-\mathbf{\bar{x}})^T (\mathbf{x}-\mathbf{\bar{x}})}{4\,k\,t}\right),$

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 s to 0.02 s.
• While the analytic solution is defined on the unbounded domain $\mathbf{R}^d$, we restrict the computational to a box or sphere. To avoid artifacts from the artificial boundaries, we use a fairly short simulation time.

## Imports and test config

Copy
%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.01

## Initial values and exact solution

def 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_exact
def 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,
)
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
)
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,
)
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 mesh

## Simulation

def 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,
)

## Error analysis

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 solution

## Scenarios

### 2D Cartesian, order 1

mesh = 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)
solution
Job job_2009012124572585_a9ed916d0d running on local_f64 with 4 rank(s).
Site information:
* Salvus version: 0.11.15
* Floating point size: 64
-> Current Task: Time loop complete* Downloaded 1.8 MB of results to diffusion.
* Total run time: 1.55 seconds.
* Pure simulation time: 1.22 seconds.
simulation time: 0.01
|| u ||_inf = 31.959188667283925
|| u_exact ||_inf = 31.830988618379067
|| u - u_exact ||_inf = 0.12820004890485848
|| u - u_exact ||_inf / || u_exact ||_inf = 0.004027523318293556


/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
yield from self._mapping

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f40a1c36090>

### 2D Cartesian, order 2

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)
solution
Job job_2009012124619078_6c559a235e running on local_f64 with 4 rank(s).
Site information:
* Salvus version: 0.11.15
* Floating point size: 64


* Downloaded 4.8 MB of results to diffusion.
* Total run time: 3.51 seconds.
* Pure simulation time: 3.07 seconds.
simulation time: 0.01
|| u ||_inf = 31.84917675381394
|| u_exact ||_inf = 31.830988618379067
|| u - u_exact ||_inf = 0.018220665804165037
|| u - u_exact ||_inf / || u_exact ||_inf = 0.0005724190983387964


/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
yield from self._mapping

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f409b41a4d0>

### 2D Cartesian, order 4

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)
solution
Job job_2009012124100013_8a0005d5fe running on local_f64 with 4 rank(s).
Site information:
* Salvus version: 0.11.15
* Floating point size: 64


* Downloaded 15.3 MB of results to diffusion.
* Total run time: 18.19 seconds.
* Pure simulation time: 17.60 seconds.
simulation time: 0.01
|| u ||_inf = 31.83316360903441
|| u_exact ||_inf = 31.830988618379067
|| u - u_exact ||_inf = 0.002174990655344544
|| u - u_exact ||_inf / || u_exact ||_inf = 6.83293466445687e-05


/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
yield from self._mapping

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f409b0f7690>

### 2D spherical, order 4

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)
solution
Job job_2009012125663910_afbb34223b running on local_f64 with 4 rank(s).
Site information:
* Salvus version: 0.11.15
* Floating point size: 64


* Downloaded 2.0 MB of results to diffusion.
* Total run time: 6.41 seconds.
* Pure simulation time: 5.35 seconds.
simulation time: 0.01
|| u ||_inf = 31.832201944166
|| u_exact ||_inf = 31.830988618379067
|| u - u_exact ||_inf = 0.0012274499227693525
|| u - u_exact ||_inf / || u_exact ||_inf = 3.856147660021557e-05


/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
yield from self._mapping

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f409b108410>
PAGE CONTENTS