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.

Diffusion Equation

We consider a spatial domain $\Omega \subset \mathbf{R}^d$ (d = 2 or 3), a time interval $I = [0, T]$, and a diffusion equation of the following form:
$m_0(\mathbf{x}) \partial_t u(\mathbf{x},t) - \nabla \cdot \big(\mathcal{D}(\mathbf{x}) \nabla u(\mathbf{x},t)\big) = 0,$
with initial conditions
$u(\mathbf{x},0) = u_{\text{init}}(\mathbf{x})$.
Here, $u$ denotes the space- and time-dependent diffusive field and $u_{\text{init}}$ are describes external forces. $\partial_t$ denotes the first time derivative and $\nabla$ the spatial gradient operator. Furthermore, the scalar parameter $m_0$ and the symmetric second-order diffusion tensor $\mathcal{D}$ are space-dependent coefficients.
$\mathcal{D}$ can be related to a Wiener process using the relation
\mathcal{D} = \frac{1}{2} \sigma \sigma\,^T,
which direction-dependent smoothing lengths $\sigma_i$.
For the special case of $m_0 = 1$ and $T = 1$, $\sigma$ corresponds to the standard deviation of the Gaussian smoothing in meters.
In the isotropic case, $\mathcal{D}$ simplifies to a scalar value, in which case we may re-write the diffusion equation as
$m_0(\mathbf{x}) \partial_t u(\mathbf{x},t) - \nabla \cdot \big(m_1(\mathbf{x}) \nabla u(\mathbf{x},t)\big) = 0,$
with $m_1 = 0.5 * \sigma^2$ and the isotropic smoothing length $\sigma$.
Copy
# Standard Python packages
import matplotlib.pyplot as plt
import numpy as np
import os
import toml

# Salvus imports
import salvus.namespace as sn

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "token")
mesh = sn.simple_mesh.basic_mesh.CartesianHomogeneousAcoustic2D(
x_max=4.0, y_max=6.0, vp=1000.0, rho=1000.0, max_frequency=5000.0
).create_mesh()

# We don't need any material parameters.
mesh.elemental_fields = {}
input_mesh = mesh.copy()
input_mesh.attach_field("some_field", np.random.randn(mesh.npoint))
input_mesh.map_nodal_fields_to_element_nodal()
input_mesh.write_h5("initial_values.h5")
input_mesh