This tutorial is not for the latest stable Salvus version (viewing 0.12.16).

Benchmark: Poisson’s equation

Reference solution: analytic
Physis: Poisson’s equation

In this notebook we benchmark the numerical solution of the Poisson equation in a simple test case for which an analytic solution exists.

The Poisson equation is given by

- \Delta u = f,

with some suitable boundary conditions. Note that there are different conventions for the sign of the Laplace operator on the left-hand-side. Salvus considers Poisson’s equation with a negative sign. The problem is only well-posed with a specification of boundary conditions, which are discussed below.

Specifically, we are following the approach of a “manufactured solution” in a homogeneous media and test different boundary conditions.

We consider $a > 1$ and the function

u(\mathbf{x}) = 2^{4a}\, \prod_{i=0}^d x_i^a\,(1-x_i)^a

on the square domain $\Omega = [0,1]^d$ . Clearly $u|_{\partial\Omega} = 0$ , i.e., $u$ satisfies homogenenous boundary conditions. Now, we set up Poisson’s equation with constant coefficients, such that $u$ is the unique solution. Hence, we compute $\Delta u$ to obtain the analytical right-hand-side $f$ as

f(\mathbf{x}) = 2^{4a}\, \sum_{j=0}^d\left[ x_i^{(a-2)}\,(1-x_i)^{(a-2)}\,\left[a((a-1)(1-2x_j)^2 - 2x_j(1-x_j) \right] \left[\prod_{i=0, i\neq j}^d x_i^a\,(1-x_i)^a\right]\right]

In the following we will set $a=10$ .

Imports and test config

%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")

Analytical solution and right-hand-side

def g(a, x):
    return (x**a) * (1 - x) ** a


def h(a, x):
    return a * (a - 1) * (1 - 2 * x) ** 2 - 2 * a * x * (1 - x)


def boundary_mask(p, tol=1e-15):
    bmask_0 = p[:, :, 0] < tol
    bmask_1 = p[:, :, 1] < tol
    bmask_2 = p[:, :, 0] > 1 - tol
    bmask_3 = p[:, :, 1] > 1 - tol
    bnd_mask = bmask_0 + bmask_1 + bmask_2 + bmask_3
    return bnd_mask


def right_hand_side(p, a=10):
    x0 = p[:, :, 0]
    x1 = p[:, :, 1]

    u_00 = (2 ** (4 * a)) * g(a, x1) * g(a - 2, x0) * h(a, x0)
    u_11 = (2 ** (4 * a)) * g(a, x0) * g(a - 2, x1) * h(a, x1)

    rhs = -(u_00 + u_11)

    assert np.max(np.abs(rhs[boundary_mask(p)])) < 1e-15
    return rhs


def exact_solution(p, a=10):
    x0 = p[:, :, 0]
    x1 = p[:, :, 1]
    u = (2 ** (4 * a)) * g(a, x0) * g(a, x1)

    assert np.max(np.abs(u[boundary_mask(p)])) < 1e-15

    return u

def create_mesh(dim, tensor_order, domain_size, frequency):
    if dim == 2:
        m = sn.simple_mesh.CartesianHomogeneousAcoustic2D(
            vp=1.0,
            rho=1.0,
            x_max=domain_size,
            y_max=domain_size,
            max_frequency=frequency,
        )
        m.advanced.tensor_order = tensor_order
        mesh = m.create_mesh()

    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=frequency,
        )
        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", f)
    mesh.attach_field("M1", f)

    mesh.attach_field("rhs", right_hand_side(mesh.points[mesh.connectivity]))

    mesh.attach_field("u_true", exact_solution(mesh.points[mesh.connectivity]))

    return mesh

Simulation

def simulate(mesh):
    mesh.write_h5("mesh.h5")
    sim = sn.simple_config.simulation.Poisson(mesh=mesh)

    sim.domain.polynomial_order = mesh.shape_order

    sim.physics.poisson_equation.right_hand_side.filename = "mesh.h5"
    sim.physics.poisson_equation.right_hand_side.format = "hdf5"
    sim.physics.poisson_equation.right_hand_side.field = "rhs"

    sim.physics.poisson_equation.solution.filename = "solution.h5"
    sim.physics.poisson_equation.solution.fields = ["solution", "residuals"]

    boundaries = sn.simple_config.boundary.HomogeneousDirichlet(
        side_sets=["x0", "x1", "y0", "y1"]
    )

    # Associate boundaries with our simulation.
    sim.add_boundary_conditions(boundaries)

    sim.solver.max_iterations = 8000
    sim.solver.absolute_tolerance = 0.0
    sim.solver.relative_tolerance = 1e-14
    sim.solver.preconditioner = True

    sim.validate()

    sn.api.run(
        input_file=sim,
        site_name=SALVUS_FLOW_SITE_NAME,
        output_folder="poisson_2d",
        overwrite=True,
        ranks=MPI_RANKS,
    )

Error analysis

def analysis(mesh):
    tensor_order = mesh.shape_order

    if tensor_order == 1:
        TOL = 1e-3
    elif tensor_order == 2:
        TOL = 1e-4
    else:
        TOL = 1e-5

    solution = sn.UnstructuredMesh.from_h5("poisson_2d/solution.h5")

    u_salvus = solution.elemental_fields["solution"]
    u_analytic = mesh.elemental_fields["u_true"]
    diff = 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)
    solution.attach_field("diff", diff)
    print(
        f"|| u ||_inf = {u_max}\n",
        f"|| u_exact ||_inf = {ref_max}\n",
        f"|| u - u_exact ||_inf = {np.abs(diff).max()}\n",
        f"|| u - u_exact ||_inf / || u_exact ||_inf = {np.abs(diff).max() / ref_max}\n",
        f"|| A u - f ||_inf = {np.abs(solution.elemental_fields['residuals']).max()}\n",
    )
    np.testing.assert_allclose(
        u_salvus, u_analytic, rtol=1e-6, atol=ref_max * TOL
    )
    np.testing.assert_allclose(
        solution.elemental_fields["residuals"], 0, rtol=1e-6, atol=1e-7
    )

    return solution

Scenarios

2D Cartesian, order 1

mesh = create_mesh(dim=2, tensor_order=1, domain_size=1.0, frequency=60.0)
simulate(mesh)
solution = analysis(mesh)
solution

SalvusJob `job_2403151021751917_740c52073d` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 1.8 MB of results to `poisson_2d`.
* Total run time: 0.96 seconds.
* Pure simulation time: 0.77 seconds.
|| u ||_inf = 1.0006818771362305
 || u_exact ||_inf = 1.0
 || u - u_exact ||_inf = 0.0006818771362304688
 || u - u_exact ||_inf / || u_exact ||_inf = 0.0006818771362304688
 || A u - f ||_inf = 5.177404061029733e-10

<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f7f09216e10>

2D Cartesian, order 2

mesh = create_mesh(dim=2, tensor_order=2, domain_size=1.0, frequency=60.0)
simulate(mesh)
solution = analysis(mesh)
solution

SalvusJob `job_2403151021024238_96bfad529c` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 4.8 MB of results to `poisson_2d`.
* Total run time: 1.94 seconds.
* Pure simulation time: 1.35 seconds.
|| u ||_inf = 0.9999993443489075
 || u_exact ||_inf = 1.0
 || u - u_exact ||_inf = 1.3003924218546459e-06
 || u - u_exact ||_inf / || u_exact ||_inf = 1.3003924218546459e-06
 || A u - f ||_inf = 2.1258116100764113e-10

<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f7f007dc550>

2D Cartesian, order 4

mesh = create_mesh(dim=2, tensor_order=4, domain_size=1.0, frequency=60.0)
simulate(mesh)
solution = analysis(mesh)
solution

SalvusJob `job_2403151021450906_d2a60cbd1c` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.12.16
  * Floating point size: 32

* Downloaded 15.4 MB of results to `poisson_2d`.
* Total run time: 3.55 seconds.
* Pure simulation time: 2.91 seconds.
|| u ||_inf = 0.999998927116394
 || u_exact ||_inf = 1.0
 || u - u_exact ||_inf = 2.1720745363440486e-06
 || u - u_exact ||_inf / || u_exact ||_inf = 2.1720745363440486e-06
 || A u - f ||_inf = 7.68556052133107e-11

<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f7f34dd14d0>