add acoustic wave problem

Author: GiovanniCanali

docs/source/_rst/_code.rst

@@ -215,6 +215,7 @@ Problems Zoo
 .. toctree::
     :titlesonly:
 
+    AcousticWaveProblem <problem/zoo/acoustic_wave.rst>
     AdvectionProblem <problem/zoo/advection.rst>
     AllenCahnProblem <problem/zoo/allen_cahn.rst>
     DiffusionReactionProblem <problem/zoo/diffusion_reaction.rst>

docs/source/_rst/problem/zoo/acoustic_wave.rst

@@ -0,0 +1,9 @@
+AcousticWaveProblem
+=====================
+.. currentmodule:: pina.problem.zoo.acoustic_wave
+
+.. automodule:: pina.problem.zoo.acoustic_wave
+
+.. autoclass:: AcousticWaveProblem
+    :members:
+    :show-inheritance:

pina/equation/__init__.py

@@ -13,6 +13,7 @@
     "DiffusionReaction",
     "Helmholtz",
     "Poisson",
+    "AcousticWave",
 ]
 
 from .equation import Equation
@@ -27,5 +28,6 @@
     DiffusionReaction,
     Helmholtz,
     Poisson,
+    AcousticWave,
 )
 from .system_equation import SystemEquation

pina/equation/equation_factory.py

@@ -6,9 +6,6 @@
 from ..operator import grad, div, laplacian
 from ..utils import check_consistency
 
-# Pylint warning disabled because the classes defined in this module
-# inherit from Equation and are meant to be simple containers for equations.
-
 
 class FixedValue(Equation):  # pylint: disable=R0903
     """
@@ -452,3 +449,60 @@ def equation(input_, output_):
             return lap - self.forcing_term(input_)
 
         super().__init__(equation)
+
+
+class AcousticWave(Equation):  # pylint: disable=R0903
+    r"""
+    Implementation of the N-dimensional isotropic acoustic wave equation.
+    The equation is defined as follows:
+
+    .. math::
+
+        \frac{\partial^2 u}{\partial t^2} - c^2 \Delta u = 0
+
+    or alternatively:
+
+    .. math::
+
+        \Box u = 0
+
+    Here, :math:`c` is the wave propagation speed, and :math:`\Box` is the
+    d'Alembert operator.
+    """
+
+    def __init__(self, c):
+        """
+        Initialization of the :class:`AcousticWaveEquation` class.
+
+        :param c: The wave propagation speed.
+        :type c: float | int
+        """
+        check_consistency(c, (float, int))
+        self.c = c
+
+        def equation(input_, output_):
+            """
+            Implementation of the acoustic wave equation.
+
+            :param LabelTensor input_: The input data of the problem.
+            :param LabelTensor output_: The output data of the problem.
+            :return: The residual of the acoustic wave equation.
+            :rtype: LabelTensor
+            :raises ValueError: If the ``input_`` labels do not contain the time
+                variable 't'.
+            """
+            # Ensure time is passed as input
+            if "t" not in input_.labels:
+                raise ValueError(
+                    "The ``input_`` labels must contain the time 't' variable."
+                )
+
+            # Compute the time second derivative and the spatial laplacian
+            u_tt = laplacian(output_, input_, d=["t"])
+            u_xx = laplacian(
+                output_, input_, d=[di for di in input_.labels if di != "t"]
+            )
+
+            return u_tt - self.c**2 * u_xx
+
+        super().__init__(equation)

pina/problem/zoo/__init__.py

@@ -8,6 +8,7 @@
     "Poisson2DSquareProblem",
     "DiffusionReactionProblem",
     "InversePoisson2DSquareProblem",
+    "AcousticWaveProblem",
 ]
 
 from .supervised_problem import SupervisedProblem
@@ -17,3 +18,4 @@
 from .poisson_2d_square import Poisson2DSquareProblem
 from .diffusion_reaction import DiffusionReactionProblem
 from .inverse_poisson_2d_square import InversePoisson2DSquareProblem
+from .acoustic_wave import AcousticWaveProblem

pina/problem/zoo/acoustic_wave.py

@@ -0,0 +1,97 @@
+"""Formulation of the acoustic wave problem."""
+
+import torch
+from ... import Condition
+from ...problem import SpatialProblem, TimeDependentProblem
+from ...utils import check_consistency
+from ...domain import CartesianDomain
+from ...equation import (
+    Equation,
+    SystemEquation,
+    FixedValue,
+    FixedGradient,
+    AcousticWave,
+)
+
+
+def initial_condition(input_, output_):
+    """
+    Definition of the initial condition of the acoustic wave problem.
+
+    :param LabelTensor input_: The input data of the problem.
+    :param LabelTensor output_: The output data of the problem.
+    :return: The residual of the initial condition.
+    :rtype: LabelTensor
+    """
+    arg = torch.pi * input_["x"]
+    return output_ - torch.sin(arg) - 0.5 * torch.sin(4 * arg)
+
+
+class AcousticWaveProblem(TimeDependentProblem, SpatialProblem):
+    r"""
+    Implementation of the acoustic wave problem in the spatial interval
+    :math:`[0, 1]` and temporal interval :math:`[0, 1]`.
+
+    .. seealso::
+
+        **Original reference**: Wang, Sifan, Xinling Yu, and
+        Paris Perdikaris. *When and why PINNs fail to train:
+        A neural tangent kernel perspective*. Journal of
+        Computational Physics 449 (2022): 110768.
+        DOI: `10.1016 <https://doi.org/10.1016/j.jcp.2021.110768>`_.
+
+    :Example:
+
+        >>> problem = AcousticWaveProblem(c=2.0)
+    """
+
+    output_variables = ["u"]
+    spatial_domain = CartesianDomain({"x": [0, 1]})
+    temporal_domain = CartesianDomain({"t": [0, 1]})
+
+    domains = {
+        "D": CartesianDomain({"x": [0, 1], "t": [0, 1]}),
+        "t0": CartesianDomain({"x": [0, 1], "t": 0.0}),
+        "g1": CartesianDomain({"x": 0.0, "t": [0, 1]}),
+        "g2": CartesianDomain({"x": 1.0, "t": [0, 1]}),
+    }
+
+    conditions = {
+        "g1": Condition(domain="g1", equation=FixedValue(value=0.0)),
+        "g2": Condition(domain="g2", equation=FixedValue(value=0.0)),
+        "t0": Condition(
+            domain="t0",
+            equation=SystemEquation(
+                [Equation(initial_condition), FixedGradient(value=0.0, d="t")]
+            ),
+        ),
+    }
+
+    def __init__(self, c=2.0):
+        """
+        Initialization of the :class:`AcousticWaveProblem` class.
+
+        :param c: The wave propagation speed. Default is 2.0.
+        :type c: float | int
+        """
+        super().__init__()
+        check_consistency(c, (float, int))
+        self.c = c
+
+        self.conditions["D"] = Condition(
+            domain="D", equation=AcousticWave(self.c)
+        )
+
+    def solution(self, pts):
+        """
+        Implementation of the analytical solution of the acoustic wave problem.
+
+        :param LabelTensor pts: Points where the solution is evaluated.
+        :return: The analytical solution of the acoustic wave problem.
+        :rtype: LabelTensor
+        """
+        arg_x = torch.pi * pts["x"]
+        arg_t = self.c * torch.pi * pts["t"]
+        term1 = torch.sin(arg_x) * torch.cos(arg_t)
+        term2 = 0.5 * torch.sin(4 * arg_x) * torch.cos(4 * arg_t)
+        return term1 + term2

tests/test_equation/test_equation_factory.py

@@ -8,6 +8,7 @@
     DiffusionReaction,
     Helmholtz,
     Poisson,
+    AcousticWave,
 )
 from pina import LabelTensor
 import torch
@@ -195,3 +196,22 @@ def test_poisson_equation(forcing_term):
     # Residual
     residual = equation.residual(pts, u)
     assert residual.shape == u.shape
+
+
+@pytest.mark.parametrize("c", [1.0, 10, -7.5])
+def test_acoustic_wave_equation(c):
+
+    # Constructor
+    equation = AcousticWave(c=c)
+
+    # Should fail if c is not a float or int
+    with pytest.raises(ValueError):
+        AcousticWave(c="invalid")
+
+    # Residual
+    residual = equation.residual(pts, u)
+    assert residual.shape == u.shape
+
+    # Should fail if the input has no 't' label
+    with pytest.raises(ValueError):
+        residual = equation.residual(pts["x", "y"], u)

tests/test_problem_zoo/test_acoustic_wave.py

@@ -0,0 +1,19 @@
+import pytest
+from pina.problem.zoo import AcousticWaveProblem
+from pina.problem import SpatialProblem, TimeDependentProblem
+
+
+@pytest.mark.parametrize("c", [0.1, 1])
+def test_constructor(c):
+
+    problem = AcousticWaveProblem(c=c)
+    problem.discretise_domain(n=10, mode="random", domains="all")
+    assert problem.are_all_domains_discretised
+    assert isinstance(problem, SpatialProblem)
+    assert isinstance(problem, TimeDependentProblem)
+    assert hasattr(problem, "conditions")
+    assert isinstance(problem.conditions, dict)
+
+    # Should fail if c is not a float or int
+    with pytest.raises(ValueError):
+        AcousticWaveProblem(c="invalid")