Mondaic
PricingAbout UsDocumentation
Mondaic
Version:
No results found.

This documentation is not for the latest stable Salvus version.

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
−Δu=f,- \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>1a > 1 and the function
u(x)=24a ∏i=0dxia (1−xi)au(\mathbf{x}) = 2^{4a}\, \prod_{i=0}^d x_i^a\,(1-x_i)^a
on the square domain Ω=[0,1]d\Omega = [0,1]^d. Clearly u∣∂Ω=0u|_{\partial\Omega} = 0, i.e., uu satisfies homogenenous boundary conditions. Now, we set up Poisson's equation with constant coefficients, such that uu is the unique solution. Hence, we compute Δu\Delta u to obtain the analytical right-hand-side ff as
f(x)=24a ∑j=0d[xi(a−2) (1−xi)(a−2) [a((a−1)(1−2xj)2−2xj(1−xj)][∏i=0,i≠jdxia (1−xi)a]]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=10a=10.

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

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>
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>
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>
Previous
Analytic test - Diffusion equation
Next
Gradient test simple physics
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