Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

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 D(x)=(k00k)\mathcal{D}(\mathbf{x}) = \left(\begin{array}{cc}k & 0\\0 & k \end{array}\right) with k>0k > 0:
∂tu(x,t)−∇⋅D(x) ∇u(x,t)=0,u(x,0)=δ(x−xˉ).\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 Rd\mathbf{R}^d, (d=2,3d=2,3), this equation has a unique solution given by:
uexact=1(kπ t)d/2exp⁡(−(x−xˉ)T(x−xˉ)4 k t),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 Rd\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
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,
        )
        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 mesh
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,
    )
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

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
SalvusJob `job_2411151453249023_6d6794bbcc` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.2
  * Floating point size: 32

* Downloaded 1.6 MB of results to `diffusion`.
* Total run time: 2.25 seconds.
* Pure simulation time: 1.87 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.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x762436358e10>
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
SalvusJob `job_2411151453748088_23ec17146e` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.2
  * Floating point size: 32
-> Current Task: Time loop complete* Downloaded 4.3 MB of results to `diffusion`.
* Total run time: 2.74 seconds.
* Pure simulation time: 2.57 seconds.
simulation time: 0.01
 || u ||_inf = 31.849178314208984
 || u_exact ||_inf = 31.830988618379067
 || u - u_exact ||_inf = 0.01822223726619754
 || u - u_exact ||_inf / || u_exact ||_inf = 0.0005724684672745635


<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x762476e28e90>
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
SalvusJob `job_2411151454905618_6d415d0a6f` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.2
  * Floating point size: 32
                                          

* Downloaded 14.0 MB of results to `diffusion`.
* Total run time: 4.08 seconds.
* Pure simulation time: 3.48 seconds.
simulation time: 0.01
 || u ||_inf = 31.83320426940918
 || u_exact ||_inf = 31.830988618379067
 || u - u_exact ||_inf = 0.0022156510301130083
 || u - u_exact ||_inf / || u_exact ||_inf = 6.960672999121685e-05


<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x762476ec7190>
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
SalvusJob `job_2411151454768521_4b1c6becaa` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 2024.1.2
  * Floating point size: 32

* Downloaded 2.5 MB of results to `diffusion`.
* Total run time: 1.46 seconds.
* Pure simulation time: 1.19 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.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x762477276b10>
Next
GETTING STARTED
We use cookies to analyze our traffic, but won't track you until explicit consent is given. Using our contact form will create relevant cookies and transmit data to third party services to keep track of your request. Read more about how we use personal information in our Privacy Policy.
PAGE CONTENTS