跳转至

Equation(方程) 模块

ppsci.equation

PDE

Base class for Partial Differential Equation.

Source code in ppsci/equation/pde/base.py
class PDE:
    """Base class for Partial Differential Equation."""

    def __init__(self):
        super().__init__()
        self.equations: Dict[str, Union[Callable, sp.Basic]] = {}
        # for PDE which has learnable parameter(s)
        self.learnable_parameters = nn.ParameterList()

        self.detach_keys: Optional[Tuple[str, ...]] = None

    @staticmethod
    def create_symbols(
        symbol_str: str,
    ) -> Union[sp.Symbol, Tuple[sp.Symbol, ...]]:
        """Create symbolic variables.

        Args:
            symbol_str (str): String contains symbols, such as "x", "x y z".

        Returns:
            Union[sympy.Symbol, Tuple[sympy.Symbol, ...]]: Created symbol(s).

        Examples:
            >>> import ppsci
            >>> pde = ppsci.equation.PDE()
            >>> symbol_x = pde.create_symbols('x')
            >>> symbols_xyz = pde.create_symbols('x y z')
            >>> print(symbol_x)
            x
            >>> print(symbols_xyz)
            (x, y, z)
        """
        return sp.symbols(symbol_str)

    def create_function(self, name: str, invars: Tuple[sp.Symbol, ...]) -> sp.Function:
        """Create named function depending on given invars.

        Args:
            name (str): Function name. such as "u", "v", and "f".
            invars (Tuple[sympy.Symbol, ...]): List of independent variable of function.

        Returns:
            sympy.Function: Named sympy function.

        Examples:
            >>> import ppsci
            >>> pde = ppsci.equation.PDE()
            >>> x, y, z = pde.create_symbols('x y z')
            >>> u = pde.create_function('u', (x, y))
            >>> f = pde.create_function('f', (x, y, z))
            >>> print(u)
            u(x, y)
            >>> print(f)
            f(x, y, z)
        """
        expr = sp.Function(name)(*invars)

        return expr

    def _apply_detach(self):
        """
        Wrap detached sub_expr into detach(sub_expr) to prevent gradient back-propagation, only for those items speicified in self.detach_keys.

        NOTE: This function is expected to be called after self.equations is ready in PDE.__init__.

        Examples:
            >>> import ppsci
            >>> ns = ppsci.equation.NavierStokes(1.0, 1.0, 2, False)
            >>> print(ns)
            NavierStokes
                continuity: Derivative(u(x, y), x) + Derivative(v(x, y), y)
                momentum_x: u(x, y)*Derivative(u(x, y), x) + v(x, y)*Derivative(u(x, y), y) + 1.0*Derivative(p(x, y), x) - 1.0*Derivative(u(x, y), (x, 2)) - 1.0*Derivative(u(x, y), (y, 2))
                momentum_y: u(x, y)*Derivative(v(x, y), x) + v(x, y)*Derivative(v(x, y), y) + 1.0*Derivative(p(x, y), y) - 1.0*Derivative(v(x, y), (x, 2)) - 1.0*Derivative(v(x, y), (y, 2))
            >>> detach_keys = ("u", "v__y")
            >>> ns = ppsci.equation.NavierStokes(1.0, 1.0, 2, False, detach_keys=detach_keys)
            >>> print(ns)
            NavierStokes
                continuity: detach(Derivative(v(x, y), y)) + Derivative(u(x, y), x)
                momentum_x: detach(u(x, y))*Derivative(u(x, y), x) + v(x, y)*Derivative(u(x, y), y) + 1.0*Derivative(p(x, y), x) - 1.0*Derivative(u(x, y), (x, 2)) - 1.0*Derivative(u(x, y), (y, 2))
                momentum_y: detach(u(x, y))*Derivative(v(x, y), x) + detach(Derivative(v(x, y), y))*v(x, y) + 1.0*Derivative(p(x, y), y) - 1.0*Derivative(v(x, y), (x, 2)) - 1.0*Derivative(v(x, y), (y, 2))
        """
        if self.detach_keys is None:
            return

        from copy import deepcopy

        from sympy.core.traversal import postorder_traversal

        from ppsci.utils.symbolic import _cvt_to_key

        for name, expr in self.equations.items():
            if not isinstance(expr, sp.Basic):
                continue
            # only process sympy expression
            expr_ = deepcopy(expr)
            for item in postorder_traversal(expr):
                if _cvt_to_key(item) in self.detach_keys:
                    # inplace all related sub_expr into detach(sub_expr)
                    expr_ = expr_.replace(item, sp.Function(DETACH_FUNC_NAME)(item))

                    # remove all detach wrapper for more-than-once wrapped items to prevent duplicated wrapping
                    expr_ = expr_.replace(
                        sp.Function(DETACH_FUNC_NAME)(
                            sp.Function(DETACH_FUNC_NAME)(item)
                        ),
                        sp.Function(DETACH_FUNC_NAME)(item),
                    )

                    # remove unccessary detach wrapping for the first arg of Derivative
                    for item_ in list(postorder_traversal(expr_)):
                        if isinstance(item_, sp.Derivative):
                            if item_.args[0].name == DETACH_FUNC_NAME:
                                expr_ = expr_.replace(
                                    item_,
                                    sp.Derivative(
                                        item_.args[0].args[0], *item_.args[1:]
                                    ),
                                )

            self.equations[name] = expr_

    def add_equation(self, name: str, equation: Callable):
        """Add an equation.

        Args:
            name (str): Name of equation
            equation (Callable): Computation function for equation.

        Examples:
            >>> import ppsci
            >>> import sympy
            >>> pde = ppsci.equation.PDE()
            >>> x, y = pde.create_symbols('x y')
            >>> u = x**2 + y**2
            >>> equation = sympy.diff(u, x) + sympy.diff(u, y)
            >>> pde.add_equation('linear_pde', equation)
            >>> print(pde)
            PDE
                linear_pde: 2*x + 2*y
        """
        self.equations.update({name: equation})

    def parameters(self) -> List[paddle.Tensor]:
        """Return learnable parameters contained in PDE.

        Returns:
            List[Tensor]: A list of learnable parameters.

        Examples:
            >>> import ppsci
            >>> pde = ppsci.equation.Vibration(2, -4, 0)
            >>> print(pde.parameters())
            [Parameter containing:
            Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
                   -4.), Parameter containing:
            Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
                   0.)]
        """
        return self.learnable_parameters.parameters()

    def state_dict(self) -> Dict[str, paddle.Tensor]:
        """Return named learnable parameters in dict.

        Returns:
            Dict[str, Tensor]: A dict of states(str) and learnable parameters(Tensor).

        Examples:
            >>> import ppsci
            >>> pde = ppsci.equation.Vibration(2, -4, 0)
            >>> print(pde.state_dict())
            OrderedDict([('0', Parameter containing:
            Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
                   -4.)), ('1', Parameter containing:
            Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
                   0.))])
        """
        return self.learnable_parameters.state_dict()

    def set_state_dict(
        self, state_dict: Dict[str, paddle.Tensor]
    ) -> Tuple[List[str], List[str]]:
        """Set state dict from dict.

        Args:
            state_dict (Dict[str, paddle.Tensor]): The state dict to be set.

        Returns:
            Tuple[List[str], List[str]]: List of missing_keys and unexpected_keys.
                Expected to be two empty tuples mostly.

        Examples:
            >>> import paddle
            >>> import ppsci
            >>> paddle.set_default_dtype("float64")
            >>> pde = ppsci.equation.Vibration(2, -4, 0)
            >>> state = pde.state_dict()
            >>> state['0'] = paddle.to_tensor(-3.1)
            >>> pde.set_state_dict(state)
            ([], [])
            >>> print(state)
            OrderedDict([('0', Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=True,
                   -3.10000000)), ('1', Parameter containing:
            Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
                   0.))])
        """
        return self.learnable_parameters.set_state_dict(state_dict)

    def __str__(self):
        return "\n".join(
            [self.__class__.__name__]
            + [f"    {name}: {eq}" for name, eq in self.equations.items()]
        )
add_equation(name, equation)

Add an equation.

Parameters:

Name Type Description Default
name str

Name of equation

required
equation Callable

Computation function for equation.

required

Examples:

>>> import ppsci
>>> import sympy
>>> pde = ppsci.equation.PDE()
>>> x, y = pde.create_symbols('x y')
>>> u = x**2 + y**2
>>> equation = sympy.diff(u, x) + sympy.diff(u, y)
>>> pde.add_equation('linear_pde', equation)
>>> print(pde)
PDE
    linear_pde: 2*x + 2*y
Source code in ppsci/equation/pde/base.py
def add_equation(self, name: str, equation: Callable):
    """Add an equation.

    Args:
        name (str): Name of equation
        equation (Callable): Computation function for equation.

    Examples:
        >>> import ppsci
        >>> import sympy
        >>> pde = ppsci.equation.PDE()
        >>> x, y = pde.create_symbols('x y')
        >>> u = x**2 + y**2
        >>> equation = sympy.diff(u, x) + sympy.diff(u, y)
        >>> pde.add_equation('linear_pde', equation)
        >>> print(pde)
        PDE
            linear_pde: 2*x + 2*y
    """
    self.equations.update({name: equation})
create_function(name, invars)

Create named function depending on given invars.

Parameters:

Name Type Description Default
name str

Function name. such as "u", "v", and "f".

required
invars Tuple[Symbol, ...]

List of independent variable of function.

required

Returns:

Type Description
Function

sympy.Function: Named sympy function.

Examples:

>>> import ppsci
>>> pde = ppsci.equation.PDE()
>>> x, y, z = pde.create_symbols('x y z')
>>> u = pde.create_function('u', (x, y))
>>> f = pde.create_function('f', (x, y, z))
>>> print(u)
u(x, y)
>>> print(f)
f(x, y, z)
Source code in ppsci/equation/pde/base.py
def create_function(self, name: str, invars: Tuple[sp.Symbol, ...]) -> sp.Function:
    """Create named function depending on given invars.

    Args:
        name (str): Function name. such as "u", "v", and "f".
        invars (Tuple[sympy.Symbol, ...]): List of independent variable of function.

    Returns:
        sympy.Function: Named sympy function.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.PDE()
        >>> x, y, z = pde.create_symbols('x y z')
        >>> u = pde.create_function('u', (x, y))
        >>> f = pde.create_function('f', (x, y, z))
        >>> print(u)
        u(x, y)
        >>> print(f)
        f(x, y, z)
    """
    expr = sp.Function(name)(*invars)

    return expr
create_symbols(symbol_str) staticmethod

Create symbolic variables.

Parameters:

Name Type Description Default
symbol_str str

String contains symbols, such as "x", "x y z".

required

Returns:

Type Description
Union[Symbol, Tuple[Symbol, ...]]

Union[sympy.Symbol, Tuple[sympy.Symbol, ...]]: Created symbol(s).

Examples:

>>> import ppsci
>>> pde = ppsci.equation.PDE()
>>> symbol_x = pde.create_symbols('x')
>>> symbols_xyz = pde.create_symbols('x y z')
>>> print(symbol_x)
x
>>> print(symbols_xyz)
(x, y, z)
Source code in ppsci/equation/pde/base.py
@staticmethod
def create_symbols(
    symbol_str: str,
) -> Union[sp.Symbol, Tuple[sp.Symbol, ...]]:
    """Create symbolic variables.

    Args:
        symbol_str (str): String contains symbols, such as "x", "x y z".

    Returns:
        Union[sympy.Symbol, Tuple[sympy.Symbol, ...]]: Created symbol(s).

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.PDE()
        >>> symbol_x = pde.create_symbols('x')
        >>> symbols_xyz = pde.create_symbols('x y z')
        >>> print(symbol_x)
        x
        >>> print(symbols_xyz)
        (x, y, z)
    """
    return sp.symbols(symbol_str)
parameters()

Return learnable parameters contained in PDE.

Returns:

Type Description
List[Tensor]

List[Tensor]: A list of learnable parameters.

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Vibration(2, -4, 0)
>>> print(pde.parameters())
[Parameter containing:
Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
       -4.), Parameter containing:
Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
       0.)]
Source code in ppsci/equation/pde/base.py
def parameters(self) -> List[paddle.Tensor]:
    """Return learnable parameters contained in PDE.

    Returns:
        List[Tensor]: A list of learnable parameters.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Vibration(2, -4, 0)
        >>> print(pde.parameters())
        [Parameter containing:
        Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
               -4.), Parameter containing:
        Tensor(shape=[], dtype=float32, place=Place(gpu:0), stop_gradient=False,
               0.)]
    """
    return self.learnable_parameters.parameters()
set_state_dict(state_dict)

Set state dict from dict.

Parameters:

Name Type Description Default
state_dict Dict[str, Tensor]

The state dict to be set.

required

Returns:

Type Description
Tuple[List[str], List[str]]

Tuple[List[str], List[str]]: List of missing_keys and unexpected_keys. Expected to be two empty tuples mostly.

Examples:

>>> import paddle
>>> import ppsci
>>> paddle.set_default_dtype("float64")
>>> pde = ppsci.equation.Vibration(2, -4, 0)
>>> state = pde.state_dict()
>>> state['0'] = paddle.to_tensor(-3.1)
>>> pde.set_state_dict(state)
([], [])
>>> print(state)
OrderedDict([('0', Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=True,
       -3.10000000)), ('1', Parameter containing:
Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
       0.))])
Source code in ppsci/equation/pde/base.py
def set_state_dict(
    self, state_dict: Dict[str, paddle.Tensor]
) -> Tuple[List[str], List[str]]:
    """Set state dict from dict.

    Args:
        state_dict (Dict[str, paddle.Tensor]): The state dict to be set.

    Returns:
        Tuple[List[str], List[str]]: List of missing_keys and unexpected_keys.
            Expected to be two empty tuples mostly.

    Examples:
        >>> import paddle
        >>> import ppsci
        >>> paddle.set_default_dtype("float64")
        >>> pde = ppsci.equation.Vibration(2, -4, 0)
        >>> state = pde.state_dict()
        >>> state['0'] = paddle.to_tensor(-3.1)
        >>> pde.set_state_dict(state)
        ([], [])
        >>> print(state)
        OrderedDict([('0', Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=True,
               -3.10000000)), ('1', Parameter containing:
        Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
               0.))])
    """
    return self.learnable_parameters.set_state_dict(state_dict)
state_dict()

Return named learnable parameters in dict.

Returns:

Type Description
Dict[str, Tensor]

Dict[str, Tensor]: A dict of states(str) and learnable parameters(Tensor).

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Vibration(2, -4, 0)
>>> print(pde.state_dict())
OrderedDict([('0', Parameter containing:
Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
       -4.)), ('1', Parameter containing:
Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
       0.))])
Source code in ppsci/equation/pde/base.py
def state_dict(self) -> Dict[str, paddle.Tensor]:
    """Return named learnable parameters in dict.

    Returns:
        Dict[str, Tensor]: A dict of states(str) and learnable parameters(Tensor).

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Vibration(2, -4, 0)
        >>> print(pde.state_dict())
        OrderedDict([('0', Parameter containing:
        Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
               -4.)), ('1', Parameter containing:
        Tensor(shape=[], dtype=float64, place=Place(gpu:0), stop_gradient=False,
               0.))])
    """
    return self.learnable_parameters.state_dict()

AllenCahn

Bases: PDE

Class for Allen-Cahn equation.

\[ \dfrac{\partial u}{\partial t} - \epsilon^2 \Delta u + 5u^3 - 5u = 0 \]

Parameters:

Name Type Description Default
eps float

Represents the characteristicscale of interfacial width, influencing the thickness and dynamics of phase boundaries.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.AllenCahn(eps=0.01)
Source code in ppsci/equation/pde/allen_cahn.py
class AllenCahn(base.PDE):
    r"""Class for Allen-Cahn equation.

    $$
    \dfrac{\partial u}{\partial t} - \epsilon^2 \Delta u + 5u^3 - 5u = 0
    $$

    Args:
        eps (float): Represents the characteristicscale of interfacial width,
            influencing the thickness and dynamics of phase boundaries.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.AllenCahn(eps=0.01)
    """

    def __init__(
        self,
        eps: float,
        detach_keys: Optional[Tuple[str, ...]] = None,
    ):
        super().__init__()
        self.detach_keys = detach_keys
        self.eps = eps
        # t, x = self.create_symbols("t x")
        # invars = (t, x, )
        # u = self.create_function("u", invars)
        # allen_cahn = u.diff(t) + 5 * u**3 - 5 * u - 0.0001 * u.diff(x, 2)

        # TODO: Pow(u,3) seems cause slightly larger L2 error than multiply(u*u*u)
        def allen_cahn(out):
            t, x = out["t"], out["x"]
            u = out["u"]
            u__t, u__x = jacobian(u, [t, x])
            u__x__x = jacobian(u__x, x)

            return u__t - (self.eps**2) * u__x__x + 5 * u * u * u - 5 * u

        self.add_equation("allen_cahn", allen_cahn)

Biharmonic

Bases: PDE

Class for biharmonic equation with supporting special load.

\[ \nabla^4 \varphi = \dfrac{q}{D} \]

Parameters:

Name Type Description Default
dim int

Dimension of equation.

required
q Union[float, str, Basic]

Load.

required
D Union[float, str]

Rigidity.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Biharmonic(2, -1.0, 1.0)
Source code in ppsci/equation/pde/biharmonic.py
class Biharmonic(base.PDE):
    r"""Class for biharmonic equation with supporting special load.

    $$
    \nabla^4 \varphi = \dfrac{q}{D}
    $$

    Args:
        dim (int): Dimension of equation.
        q (Union[float, str, sympy.Basic]): Load.
        D (Union[float, str]): Rigidity.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Biharmonic(2, -1.0, 1.0)
    """

    def __init__(
        self,
        dim: int,
        q: Union[float, str, sympy.Basic],
        D: Union[float, str],
        detach_keys: Optional[Tuple[str, ...]] = None,
    ):
        super().__init__()
        self.detach_keys = detach_keys

        invars = self.create_symbols("x y z")[:dim]
        u = self.create_function("u", invars)

        if isinstance(q, str):
            q = self.create_function("q", invars)
        if isinstance(D, str):
            D = self.create_function("D", invars)

        self.dim = dim
        self.q = q
        self.D = D

        biharmonic = -self.q / self.D
        for invar_i in invars:
            for invar_j in invars:
                biharmonic += u.diff(invar_i, 2).diff(invar_j, 2)

        self.add_equation("biharmonic", biharmonic)

        self._apply_detach()

FractionalPoisson

Bases: PDE

(TODO)Docstring of this class will be refined in the future.

Parameters:

Name Type Description Default
alpha float

Alpha.

required
geom Geometry

Computation geometry.

required
resolution Tuple[int, ...]

Resolution.

required

Examples:

>>> import ppsci
>>> geom_disk = ppsci.geometry.Disk([0, 0], 1)
>>> ALPHA = 0.5
>>> fpde = ppsci.equation.FractionalPoisson(ALPHA, geom_disk, [8, 100])
Source code in ppsci/equation/fpde/fractional_poisson.py
class FractionalPoisson(PDE):
    """(TODO)Docstring of this class will be refined in the future.

    Args:
        alpha (float): Alpha.
        geom (geometry.Geometry): Computation geometry.
        resolution (Tuple[int, ...]): Resolution.

    Examples:
        >>> import ppsci
        >>> geom_disk = ppsci.geometry.Disk([0, 0], 1)
        >>> ALPHA = 0.5
        >>> fpde = ppsci.equation.FractionalPoisson(ALPHA, geom_disk, [8, 100])
    """

    dtype = paddle.get_default_dtype()

    def __init__(
        self, alpha: float, geom: geometry.Geometry, resolution: Tuple[int, ...]
    ):
        super().__init__()
        self.alpha = alpha
        self.geom = geom
        self.resolution = resolution
        self._w_init = self._init_weights()

        def compute_fpde_func(out):
            x = paddle.concat((out["x"], out["y"]), axis=1)
            y = out["u"]
            indices, values, shape = self.int_mat
            int_mat = sparse.sparse_coo_tensor(
                [[p[0] for p in indices], [p[1] for p in indices]],
                values,
                shape,
                stop_gradient=False,
            )
            lhs = sparse.matmul(int_mat, y)
            lhs = lhs[:, 0]
            lhs *= (
                special.gamma((1 - self.alpha) / 2)
                * special.gamma((2 + self.alpha) / 2)
                / (2 * np.pi**1.5)
            )
            x = x[: paddle.numel(lhs)]
            rhs = (
                2**self.alpha
                * special.gamma(2 + self.alpha / 2)
                * special.gamma(1 + self.alpha / 2)
                * (1 - (1 + self.alpha / 2) * paddle.sum(x**2, axis=1))
            )
            res = lhs - rhs
            return res

        self.add_equation("fpde", compute_fpde_func)

    def _init_weights(self):
        n = self._dynamic_dist2npts(self.geom.diam) + 1
        w = [1.0]
        for j in range(1, n):
            w.append(w[-1] * (j - 1 - self.alpha) / j)
        return np.array(w, dtype=self.dtype)

    def get_x(self, x_f):
        if hasattr(self, "train_x"):
            return self.train_x

        self.x0 = x_f
        if np.any(self.geom.on_boundary(self.x0)):
            raise ValueError("x0 contains boundary points.")

        if self.geom.ndim == 1:
            dirns, dirn_w = [-1, 1], [1, 1]
        elif self.geom.ndim == 2:
            gauss_x, gauss_w = np.polynomial.legendre.leggauss(self.resolution[0])
            gauss_x, gauss_w = gauss_x.astype(self.dtype), gauss_w.astype(self.dtype)
            thetas = np.pi * gauss_x + np.pi
            dirns = np.vstack((np.cos(thetas), np.sin(thetas))).T
            dirn_w = np.pi * gauss_w
        elif self.geom.ndim == 3:
            gauss_x, gauss_w = np.polynomial.legendre.leggauss(max(self.resolution[:2]))
            gauss_x, gauss_w = gauss_x.astype(self.dtype), gauss_w.astype(self.dtype)
            thetas = (np.pi * gauss_x[: self.resolution[0]] + np.pi) / 2
            phis = np.pi * gauss_x[: self.resolution[1]] + np.pi
            dirns, dirn_w = [], []
            for i in range(self.resolution[0]):
                for j in range(self.resolution[1]):
                    dirns.append(
                        [
                            np.sin(thetas[i]) * np.cos(phis[j]),
                            np.sin(thetas[i]) * np.sin(phis[j]),
                            np.cos(thetas[i]),
                        ]
                    )
                    dirn_w.append(gauss_w[i] * gauss_w[j] * np.sin(thetas[i]))
            dirn_w = np.pi**2 / 2 * np.array(dirn_w)

        x, self.w = [], []
        for x0i in self.x0:
            xi = list(
                map(
                    lambda dirn: self.background_points(
                        x0i, dirn, self._dynamic_dist2npts, 0
                    ),
                    dirns,
                )
            )
            wi = list(
                map(
                    lambda i: dirn_w[i]
                    * np.linalg.norm(xi[i][1] - xi[i][0]) ** (-self.alpha)
                    * self.get_weight(len(xi[i]) - 1),
                    range(len(dirns)),
                )
            )
            # first order
            # xi, wi = zip(self.modify_first_order(xij, wij) for xij, wij in zip(xi, wi))
            xi, wi = zip(*map(self.modify_first_order, xi, wi))
            # second order
            # xi, wi = zip(*map(self.modify_second_order, xi, wi))
            # third order
            # xi, wi = zip(*map(self.modify_third_order, xi, wi))
            x.append(np.vstack(xi))
            self.w.append(np.hstack(wi))
        self.x = np.vstack([self.x0] + x)
        self.int_mat = self._get_int_matrix(self.x0)
        self.train_x = misc.convert_to_dict(self.x, ("x", "y"))
        return self.train_x

    def get_weight(self, n):
        return self._w_init[: n + 1]

    def background_points(self, x, dirn, dist2npt, shift):
        dirn = dirn / np.linalg.norm(dirn)
        dx = self.distance2boundary_unitdirn(x, -dirn)
        n = max(dist2npt(dx), 1)
        h = dx / n
        pts = x - np.arange(-shift, n - shift + 1, dtype=self.dtype)[:, None] * h * dirn
        return pts

    def distance2boundary_unitdirn(self, x, dirn):
        # https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
        xc = x - self.geom.center
        xc = xc
        ad = np.dot(xc, dirn)
        return (
            -ad + (ad**2 - np.sum(xc * xc, axis=-1) + self.geom.radius**2) ** 0.5
        ).astype(self.dtype)

    def modify_first_order(self, x, w):
        x = np.vstack(([2 * x[0] - x[1]], x[:-1]))
        if not self.geom.is_inside(x[0:1])[0]:
            return x[1:], w[1:]
        return x, w

    def _dynamic_dist2npts(self, dx):
        return int(math.ceil(self.resolution[-1] * dx))

    def _get_int_matrix(self, x: np.ndarray) -> np.ndarray:
        dense_shape = (x.shape[0], self.x.shape[0])
        indices, values = [], []
        beg = x.shape[0]
        for i in range(x.shape[0]):
            for _ in range(self.w[i].shape[0]):
                indices.append([i, beg])
                beg += 1
            values = np.hstack((values, self.w[i]))
        return indices, values.astype(self.dtype), dense_shape

HeatExchanger

Bases: PDE

Class for heat exchanger equation.

\[ \begin{aligned} & L\left(\frac{q_m c_p}{v}\right)_{\mathrm{c}} \frac{\partial T_{\mathrm{c}}}{\partial \tau}-L\left(q_m c_p\right)_{\mathrm{c}} \frac{\partial T_{\mathrm{c}}}{\partial x}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{c}}\left(T_{\mathrm{w}}-T_{\mathrm{c}}\right), \\ & L\left(\frac{q_m c_p}{v}\right)_{\mathrm{h}} \frac{\partial T_{\mathrm{h}}}{\partial \tau}+L\left(q_m c_p\right)_{\mathrm{h}} \frac{\partial T_{\mathrm{h}}}{\partial x}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{h}}\left(T_{\mathrm{w}}-T_{\mathrm{h}}\right), \\ & \left(M c_p\right)_{\mathrm{w}} \frac{\partial T_{\mathrm{w}}}{\partial \tau}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{h}}\left(T_{\mathrm{h}}-T_{\mathrm{w}}\right)+\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{c}}\left(T_{\mathrm{c}}-T_{\mathrm{w}}\right). \end{aligned} \]

where:

  • \(T\) is temperature,
  • \(q_m\) is mass flow rate,
  • \(c_p\) represents specific heat capacity,
  • \(v\) denotes flow velocity,
  • \(L\) stands for flow length,
  • \(\eta_{\mathrm{o}}\) signifies fin surface efficiency,
  • \(\alpha\) stands for heat transfer coefficient,
  • \(A\) indicates heat transfer area,
  • \(M\) represents the mass of the heat transfer structure,
  • \(\tau\) correspond to time,
  • \(x\) correspond flow direction,
  • Subscripts \(\mathrm{h}\), \(\mathrm{c}\), and \(\mathrm{w}\) denote the hot fluid side, cold fluid side, and heat transfer wall, respectively.

Parameters:

Name Type Description Default
alpha_h Union[float, str]

\(\frac{(\eta_o\alpha A)_h}{L(c_p)_h}\)

required
alpha_c Union[float, str]

\(\frac{(\eta_o\alpha A)_c}{L(c_p)_c}\)

required
v_h Union[float, str]

\(v_h\)

required
v_c Union[float, str]

\(v_c\)

required
w_h Union[float, str]

\(\frac{(\eta_o\alpha A)_h}{M(c_p)_w}\)

required
w_c Union[float, str]

\(\frac{(\eta_o\alpha A)_c}{M(c_p)_w}\)

required

Examples:

>>> import ppsci
>>> pde = ppsci.equation.HeatExchanger(1.0,1.0,1.0,1.0,1.0,1.0)
Source code in ppsci/equation/pde/heat_exchanger.py
class HeatExchanger(base.PDE):
    r"""Class for heat exchanger equation.

    $$
    \begin{aligned}
    & L\left(\frac{q_m c_p}{v}\right)_{\mathrm{c}} \frac{\partial T_{\mathrm{c}}}{\partial \tau}-L\left(q_m c_p\right)_{\mathrm{c}} \frac{\partial T_{\mathrm{c}}}{\partial x}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{c}}\left(T_{\mathrm{w}}-T_{\mathrm{c}}\right), \\
    & L\left(\frac{q_m c_p}{v}\right)_{\mathrm{h}} \frac{\partial T_{\mathrm{h}}}{\partial \tau}+L\left(q_m c_p\right)_{\mathrm{h}} \frac{\partial T_{\mathrm{h}}}{\partial x}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{h}}\left(T_{\mathrm{w}}-T_{\mathrm{h}}\right), \\
    & \left(M c_p\right)_{\mathrm{w}} \frac{\partial T_{\mathrm{w}}}{\partial \tau}=\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{h}}\left(T_{\mathrm{h}}-T_{\mathrm{w}}\right)+\left(\eta_{\mathrm{o}} \alpha A\right)_{\mathrm{c}}\left(T_{\mathrm{c}}-T_{\mathrm{w}}\right).
    \end{aligned}
    $$

    where:

    - $T$ is temperature,
    - $q_m$ is mass flow rate,
    - $c_p$ represents specific heat capacity,
    - $v$ denotes flow velocity,
    - $L$ stands for flow length,
    - $\eta_{\mathrm{o}}$ signifies fin surface efficiency,
    - $\alpha$ stands for heat transfer coefficient,
    - $A$ indicates heat transfer area,
    - $M$ represents the mass of the heat transfer structure,
    - $\tau$ correspond to time,
    - $x$ correspond flow direction,
    - Subscripts $\mathrm{h}$, $\mathrm{c}$, and $\mathrm{w}$ denote the hot fluid side, cold fluid side, and heat transfer wall, respectively.

    Args:
        alpha_h: $\frac{(\eta_o\alpha A)_h}{L(c_p)_h}$
        alpha_c: $\frac{(\eta_o\alpha A)_c}{L(c_p)_c}$
        v_h: $v_h$
        v_c: $v_c$
        w_h: $\frac{(\eta_o\alpha A)_h}{M(c_p)_w}$
        w_c: $\frac{(\eta_o\alpha A)_c}{M(c_p)_w}$

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.HeatExchanger(1.0,1.0,1.0,1.0,1.0,1.0)
    """

    def __init__(
        self,
        alpha_h: Union[float, str],
        alpha_c: Union[float, str],
        v_h: Union[float, str],
        v_c: Union[float, str],
        w_h: Union[float, str],
        w_c: Union[float, str],
    ):
        super().__init__()
        x, t, qm_h, qm_c = self.create_symbols("x t qm_h qm_c")

        T_h = self.create_function("T_h", (x, t, qm_h))
        T_c = self.create_function("T_c", (x, t, qm_c))
        T_w = self.create_function("T_w", (x, t))

        T_h_x = T_h.diff(x)
        T_h_t = T_h.diff(t)
        T_c_x = T_c.diff(x)
        T_c_t = T_c.diff(t)
        T_w_t = T_w.diff(t)

        beta_h = (alpha_h * v_h) / qm_h
        beta_c = (alpha_c * v_c) / qm_c

        heat_boundary = T_h_t + v_h * T_h_x - beta_h * (T_w - T_h)
        cold_boundary = T_c_t - v_c * T_c_x - beta_c * (T_w - T_c)
        wall = T_w_t - w_h * (T_h - T_w) - w_c * (T_c - T_w)

        self.add_equation("heat_boundary", heat_boundary)
        self.add_equation("cold_boundary", cold_boundary)
        self.add_equation("wall", wall)

        self._apply_detach()

Laplace

Bases: PDE

Class for laplace equation.

\[ \nabla^2 \varphi = 0 \]

Parameters:

Name Type Description Default
dim int

Dimension of equation.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Laplace(2)
Source code in ppsci/equation/pde/laplace.py
class Laplace(base.PDE):
    r"""Class for laplace equation.

    $$
    \nabla^2 \varphi = 0
    $$

    Args:
        dim (int): Dimension of equation.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Laplace(2)
    """

    def __init__(self, dim: int, detach_keys: Optional[Tuple[str, ...]] = None):
        super().__init__()
        self.detach_keys = detach_keys

        invars = self.create_symbols("x y z")[:dim]
        u = self.create_function("u", invars)

        self.dim = dim

        laplace = 0
        for invar in invars:
            laplace += u.diff(invar, 2)

        self.add_equation("laplace", laplace)

        self._apply_detach()

LinearElasticity

Bases: PDE

Linear elasticity equations.

Use either (E, nu) or (lambda_, mu) to define the material properties.

\[ \begin{cases} stress\_disp_{xx} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial u}{\partial x} - \sigma_{xx} \\ stress\_disp_{yy} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial v}{\partial y} - \sigma_{yy} \\ stress\_disp_{zz} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial w}{\partial z} - \sigma_{zz} \\ stress\_disp_{xy} = \mu(\dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x}) - \sigma_{xy} \\ stress\_disp_{xz} = \mu(\dfrac{\partial u}{\partial z} + \dfrac{\partial w}{\partial x}) - \sigma_{xz} \\ stress\_disp_{yz} = \mu(\dfrac{\partial v}{\partial z} + \dfrac{\partial w}{\partial y}) - \sigma_{yz} \\ equilibrium_{x} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xx}}{\partial x} + \dfrac{\partial \sigma_{xy}}{\partial y} + \dfrac{\partial \sigma_{xz}}{\partial z}) \\ equilibrium_{y} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xy}}{\partial x} + \dfrac{\partial \sigma_{yy}}{\partial y} + \dfrac{\partial \sigma_{yz}}{\partial z}) \\ equilibrium_{z} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xz}}{\partial x} + \dfrac{\partial \sigma_{yz}}{\partial y} + \dfrac{\partial \sigma_{zz}}{\partial z}) \\ \end{cases} \]

Parameters:

Name Type Description Default
E Optional[Union[float, str]]

The Young's modulus. Defaults to None.

None
nu Optional[Union[float, str]]

The Poisson's ratio. Defaults to None.

None
lambda_ Optional[Union[float, str]]

Lamé's first parameter. Defaults to None.

None
mu Optional[Union[float, str]]

Lamé's second parameter (shear modulus). Defaults to None.

None
rho Union[float, str]

Mass density. Defaults to 1.

1
dim int

Dimension of the linear elasticity (2 or 3). Defaults to 3.

3
time bool

Whether contains time data. Defaults to False.

False
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.LinearElasticity(
...     E=None, nu=None, lambda_=1e4, mu=100, dim=3
... )
Source code in ppsci/equation/pde/linear_elasticity.py
class LinearElasticity(base.PDE):
    r"""Linear elasticity equations.

    Use either (E, nu) or (lambda_, mu) to define the material properties.

    $$
    \begin{cases}
        stress\_disp_{xx} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial u}{\partial x} - \sigma_{xx} \\
        stress\_disp_{yy} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial v}{\partial y} - \sigma_{yy} \\
        stress\_disp_{zz} = \lambda(\dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z}) + 2\mu \dfrac{\partial w}{\partial z} - \sigma_{zz} \\
        stress\_disp_{xy} = \mu(\dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x}) - \sigma_{xy} \\
        stress\_disp_{xz} = \mu(\dfrac{\partial u}{\partial z} + \dfrac{\partial w}{\partial x}) - \sigma_{xz} \\
        stress\_disp_{yz} = \mu(\dfrac{\partial v}{\partial z} + \dfrac{\partial w}{\partial y}) - \sigma_{yz} \\
        equilibrium_{x} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xx}}{\partial x} + \dfrac{\partial \sigma_{xy}}{\partial y} + \dfrac{\partial \sigma_{xz}}{\partial z}) \\
        equilibrium_{y} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xy}}{\partial x} + \dfrac{\partial \sigma_{yy}}{\partial y} + \dfrac{\partial \sigma_{yz}}{\partial z}) \\
        equilibrium_{z} = \rho \dfrac{\partial^2 u}{\partial t^2} - (\dfrac{\partial \sigma_{xz}}{\partial x} + \dfrac{\partial \sigma_{yz}}{\partial y} + \dfrac{\partial \sigma_{zz}}{\partial z}) \\
    \end{cases}
    $$

    Args:
        E (Optional[Union[float, str]]): The Young's modulus. Defaults to None.
        nu (Optional[Union[float, str]]): The Poisson's ratio. Defaults to None.
        lambda_ (Optional[Union[float, str]]): Lamé's first parameter. Defaults to None.
        mu (Optional[Union[float, str]]): Lamé's second parameter (shear modulus). Defaults to None.
        rho (Union[float, str], optional): Mass density. Defaults to 1.
        dim (int, optional): Dimension of the linear elasticity (2 or 3). Defaults to 3.
        time (bool, optional): Whether contains time data. Defaults to False.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.LinearElasticity(
        ...     E=None, nu=None, lambda_=1e4, mu=100, dim=3
        ... )
    """

    def __init__(
        self,
        E: Optional[Union[float, str]] = None,
        nu: Optional[Union[float, str]] = None,
        lambda_: Optional[Union[float, str]] = None,
        mu: Optional[Union[float, str]] = None,
        rho: Union[float, str] = 1,
        dim: int = 3,
        time: bool = False,
        detach_keys: Optional[Tuple[str, ...]] = None,
    ):
        super().__init__()
        self.detach_keys = detach_keys
        self.dim = dim
        self.time = time

        t, x, y, z = self.create_symbols("t x y z")
        normal_x, normal_y, normal_z = self.create_symbols("normal_x normal_y normal_z")
        invars = (x, y)
        if time:
            invars = (t,) + invars
        if self.dim == 3:
            invars += (z,)

        u = self.create_function("u", invars)
        v = self.create_function("v", invars)
        w = self.create_function("w", invars) if dim == 3 else sp.Number(0)

        sigma_xx = self.create_function("sigma_xx", invars)
        sigma_yy = self.create_function("sigma_yy", invars)
        sigma_xy = self.create_function("sigma_xy", invars)
        sigma_zz = (
            self.create_function("sigma_zz", invars) if dim == 3 else sp.Number(0)
        )
        sigma_xz = (
            self.create_function("sigma_xz", invars) if dim == 3 else sp.Number(0)
        )
        sigma_yz = (
            self.create_function("sigma_yz", invars) if dim == 3 else sp.Number(0)
        )

        # compute lambda and mu
        if lambda_ is None:
            if isinstance(nu, str):
                nu = self.create_function(nu, invars)
            if isinstance(E, str):
                E = self.create_function(E, invars)
            lambda_ = nu * E / ((1 + nu) * (1 - 2 * nu))
            mu = E / (2 * (1 + nu))
        else:
            if isinstance(lambda_, str):
                lambda_ = self.create_function(lambda_, invars)
            if isinstance(mu, str):
                mu = self.create_function(mu, invars)

        if isinstance(rho, str):
            rho = self.create_function(rho, invars)

        self.E = E
        self.nu = nu
        self.lambda_ = lambda_
        self.mu = mu
        self.rho = rho

        # compute stress equations
        stress_disp_xx = (
            lambda_ * (u.diff(x) + v.diff(y) + w.diff(z))
            + 2 * mu * u.diff(x)
            - sigma_xx
        )
        stress_disp_yy = (
            lambda_ * (u.diff(x) + v.diff(y) + w.diff(z))
            + 2 * mu * v.diff(y)
            - sigma_yy
        )
        stress_disp_zz = (
            lambda_ * (u.diff(x) + v.diff(y) + w.diff(z))
            + 2 * mu * w.diff(z)
            - sigma_zz
        )
        stress_disp_xy = mu * (u.diff(y) + v.diff(x)) - sigma_xy
        stress_disp_xz = mu * (u.diff(z) + w.diff(x)) - sigma_xz
        stress_disp_yz = mu * (v.diff(z) + w.diff(y)) - sigma_yz

        # compute equilibrium equations
        equilibrium_x = rho * ((u.diff(t)).diff(t)) - (
            sigma_xx.diff(x) + sigma_xy.diff(y) + sigma_xz.diff(z)
        )
        equilibrium_y = rho * ((v.diff(t)).diff(t)) - (
            sigma_xy.diff(x) + sigma_yy.diff(y) + sigma_yz.diff(z)
        )
        equilibrium_z = rho * ((w.diff(t)).diff(t)) - (
            sigma_xz.diff(x) + sigma_yz.diff(y) + sigma_zz.diff(z)
        )

        # compute traction equations
        traction_x = normal_x * sigma_xx + normal_y * sigma_xy + normal_z * sigma_xz
        traction_y = normal_x * sigma_xy + normal_y * sigma_yy + normal_z * sigma_yz
        traction_z = normal_x * sigma_xz + normal_y * sigma_yz + normal_z * sigma_zz

        # add stress equations
        self.add_equation("stress_disp_xx", stress_disp_xx)
        self.add_equation("stress_disp_yy", stress_disp_yy)
        self.add_equation("stress_disp_xy", stress_disp_xy)
        if self.dim == 3:
            self.add_equation("stress_disp_zz", stress_disp_zz)
            self.add_equation("stress_disp_xz", stress_disp_xz)
            self.add_equation("stress_disp_yz", stress_disp_yz)

        # add equilibrium equations
        self.add_equation("equilibrium_x", equilibrium_x)
        self.add_equation("equilibrium_y", equilibrium_y)
        if self.dim == 3:
            self.add_equation("equilibrium_z", equilibrium_z)

        # add traction equations
        self.add_equation("traction_x", traction_x)
        self.add_equation("traction_y", traction_y)
        if self.dim == 3:
            self.add_equation("traction_z", traction_z)

        self._apply_detach()

NavierStokes

Bases: PDE

Class for navier-stokes equation.

\[ \begin{cases} \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z} = 0 \\ \dfrac{\partial u}{\partial t} + u\dfrac{\partial u}{\partial x} + v\dfrac{\partial u}{\partial y} + w\dfrac{\partial u}{\partial z} = - \dfrac{1}{\rho}\dfrac{\partial p}{\partial x} + \nu( \dfrac{\partial ^2 u}{\partial x ^2} + \dfrac{\partial ^2 u}{\partial y ^2} + \dfrac{\partial ^2 u}{\partial z ^2} ) \\ \dfrac{\partial v}{\partial t} + u\dfrac{\partial v}{\partial x} + v\dfrac{\partial v}{\partial y} + w\dfrac{\partial v}{\partial z} = - \dfrac{1}{\rho}\dfrac{\partial p}{\partial y} + \nu( \dfrac{\partial ^2 v}{\partial x ^2} + \dfrac{\partial ^2 v}{\partial y ^2} + \dfrac{\partial ^2 v}{\partial z ^2} ) \\ \dfrac{\partial w}{\partial t} + u\dfrac{\partial w}{\partial x} + v\dfrac{\partial w}{\partial y} + w\dfrac{\partial w}{\partial z} = - \dfrac{1}{\rho}\dfrac{\partial p}{\partial z} + \nu( \dfrac{\partial ^2 w}{\partial x ^2} + \dfrac{\partial ^2 w}{\partial y ^2} + \dfrac{\partial ^2 w}{\partial z ^2} ) \\ \end{cases} \]

Parameters:

Name Type Description Default
nu Union[float, str]

Dynamic viscosity.

required
rho Union[float, str]

Density.

required
dim int

Dimension of equation.

required
time bool

Whether the equation is time-dependent.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.NavierStokes(0.1, 1.0, 3, False)
Source code in ppsci/equation/pde/navier_stokes.py
class NavierStokes(base.PDE):
    r"""Class for navier-stokes equation.

    $$
    \begin{cases}
        \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} + \dfrac{\partial w}{\partial z} = 0 \\
        \dfrac{\partial u}{\partial t} + u\dfrac{\partial u}{\partial x} + v\dfrac{\partial u}{\partial y} + w\dfrac{\partial u}{\partial z} =
            - \dfrac{1}{\rho}\dfrac{\partial p}{\partial x}
            + \nu(
                \dfrac{\partial ^2 u}{\partial x ^2}
                + \dfrac{\partial ^2 u}{\partial y ^2}
                + \dfrac{\partial ^2 u}{\partial z ^2}
            ) \\
        \dfrac{\partial v}{\partial t} + u\dfrac{\partial v}{\partial x} + v\dfrac{\partial v}{\partial y} + w\dfrac{\partial v}{\partial z} =
            - \dfrac{1}{\rho}\dfrac{\partial p}{\partial y}
            + \nu(
                \dfrac{\partial ^2 v}{\partial x ^2}
                + \dfrac{\partial ^2 v}{\partial y ^2}
                + \dfrac{\partial ^2 v}{\partial z ^2}
            ) \\
        \dfrac{\partial w}{\partial t} + u\dfrac{\partial w}{\partial x} + v\dfrac{\partial w}{\partial y} + w\dfrac{\partial w}{\partial z} =
            - \dfrac{1}{\rho}\dfrac{\partial p}{\partial z}
            + \nu(
                \dfrac{\partial ^2 w}{\partial x ^2}
                + \dfrac{\partial ^2 w}{\partial y ^2}
                + \dfrac{\partial ^2 w}{\partial z ^2}
            ) \\
    \end{cases}
    $$

    Args:
        nu (Union[float, str]): Dynamic viscosity.
        rho (Union[float, str]): Density.
        dim (int): Dimension of equation.
        time (bool): Whether the equation is time-dependent.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.NavierStokes(0.1, 1.0, 3, False)
    """

    def __init__(
        self,
        nu: Union[float, str],
        rho: Union[float, str],
        dim: int,
        time: bool,
        detach_keys: Optional[Tuple[str, ...]] = None,
    ):
        super().__init__()
        self.detach_keys = detach_keys
        self.dim = dim
        self.time = time

        t, x, y, z = self.create_symbols("t x y z")
        invars = (x, y)
        if time:
            invars = (t,) + invars
        if dim == 3:
            invars += (z,)

        if isinstance(nu, str):
            nu = sp_parser.parse_expr(nu)
            if isinstance(nu, sp.Symbol):
                invars += (nu,)

        if isinstance(rho, str):
            rho = sp_parser.parse_expr(rho)
            if isinstance(rho, sp.Symbol):
                invars += (rho,)

        self.nu = nu
        self.rho = rho

        u = self.create_function("u", invars)
        v = self.create_function("v", invars)
        w = self.create_function("w", invars) if dim == 3 else sp.Number(0)
        p = self.create_function("p", invars)

        continuity = u.diff(x) + v.diff(y) + w.diff(z)
        momentum_x = (
            u.diff(t)
            + u * u.diff(x)
            + v * u.diff(y)
            + w * u.diff(z)
            - (
                (nu * u.diff(x)).diff(x)
                + (nu * u.diff(y)).diff(y)
                + (nu * u.diff(z)).diff(z)
            )
            + 1 / rho * p.diff(x)
        )
        momentum_y = (
            v.diff(t)
            + u * v.diff(x)
            + v * v.diff(y)
            + w * v.diff(z)
            - (
                (nu * v.diff(x)).diff(x)
                + (nu * v.diff(y)).diff(y)
                + (nu * v.diff(z)).diff(z)
            )
            + 1 / rho * p.diff(y)
        )
        momentum_z = (
            w.diff(t)
            + u * w.diff(x)
            + v * w.diff(y)
            + w * w.diff(z)
            - (
                (nu * w.diff(x)).diff(x)
                + (nu * w.diff(y)).diff(y)
                + (nu * w.diff(z)).diff(z)
            )
            + 1 / rho * p.diff(z)
        )
        self.add_equation("continuity", continuity)
        self.add_equation("momentum_x", momentum_x)
        self.add_equation("momentum_y", momentum_y)
        if self.dim == 3:
            self.add_equation("momentum_z", momentum_z)

        self._apply_detach()

NormalDotVec

Bases: PDE

Normal Dot Vector.

\[ \mathbf{n} \cdot \mathbf{v} = 0 \]

Parameters:

Name Type Description Default
vec_keys Tuple[str, ...]

Keys for vectors, such as ("u", "v", "w") for velocity vector.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.NormalDotVec(("u", "v", "w"))
Source code in ppsci/equation/pde/normal_dot_vec.py
class NormalDotVec(base.PDE):
    r"""Normal Dot Vector.

    $$
    \mathbf{n} \cdot \mathbf{v} = 0
    $$

    Args:
        vec_keys (Tuple[str, ...]): Keys for vectors, such as ("u", "v", "w") for
            velocity vector.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.NormalDotVec(("u", "v", "w"))
    """

    def __init__(
        self, vec_keys: Tuple[str, ...], detach_keys: Optional[Tuple[str, ...]] = None
    ):
        super().__init__()
        self.detach_keys = detach_keys
        if not vec_keys:
            raise ValueError(f"len(vec_keys)({len(vec_keys)}) should be larger than 0.")

        self.vec_keys = vec_keys
        vec_vars = self.create_symbols(" ".join(vec_keys))
        normals = self.create_symbols("normal_x normal_y normal_z")

        normal_dot_vec = 0
        for (normal, vec) in zip(normals, vec_vars):
            normal_dot_vec += normal * vec

        self.add_equation("normal_dot_vec", normal_dot_vec)

        self._apply_detach()

Poisson

Bases: PDE

Class for poisson equation.

\[ \nabla^2 \varphi = C \]

Parameters:

Name Type Description Default
dim int

Dimension of equation.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Poisson(2)
Source code in ppsci/equation/pde/poisson.py
class Poisson(base.PDE):
    r"""Class for poisson equation.

    $$
    \nabla^2 \varphi = C
    $$

    Args:
        dim (int): Dimension of equation.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Poisson(2)
    """

    def __init__(self, dim: int, detach_keys: Optional[Tuple[str, ...]] = None):
        super().__init__()
        self.detach_keys = detach_keys
        invars = self.create_symbols("x y z")[:dim]
        p = self.create_function("p", invars)
        self.dim = dim

        poisson = 0
        for invar in invars:
            poisson += p.diff(invar, 2)

        self.add_equation("poisson", poisson)

        self._apply_detach()

Vibration

Bases: PDE

Vortex induced vibration equation.

\[ \rho \dfrac{\partial^2 \eta}{\partial t^2} + e^{k1} \dfrac{\partial \eta}{\partial t} + e^{k2} \eta = f \]

Parameters:

Name Type Description Default
rho float

Generalized mass.

required
k1 float

Learnable parameter for modal damping.

required
k2 float

Learnable parameter for generalized stiffness.

required

Examples:

>>> import ppsci
>>> pde = ppsci.equation.Vibration(1.0, 4.0, -1.0)
Source code in ppsci/equation/pde/viv.py
class Vibration(base.PDE):
    r"""Vortex induced vibration equation.

    $$
    \rho \dfrac{\partial^2 \eta}{\partial t^2} + e^{k1} \dfrac{\partial \eta}{\partial t} + e^{k2} \eta = f
    $$

    Args:
        rho (float): Generalized mass.
        k1 (float): Learnable parameter for modal damping.
        k2 (float): Learnable parameter for generalized stiffness.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.Vibration(1.0, 4.0, -1.0)
    """

    def __init__(self, rho: float, k1: float, k2: float):
        super().__init__()
        self.rho = rho
        self.k1 = paddle.create_parameter(
            shape=[],
            dtype=paddle.get_default_dtype(),
            default_initializer=initializer.Constant(k1),
        )
        self.k2 = paddle.create_parameter(
            shape=[],
            dtype=paddle.get_default_dtype(),
            default_initializer=initializer.Constant(k2),
        )
        self.learnable_parameters.append(self.k1)
        self.learnable_parameters.append(self.k2)

        t_f = self.create_symbols("t_f")
        eta = self.create_function("eta", (t_f,))
        k1 = self.create_symbols(self.k1.name)
        k2 = self.create_symbols(self.k2.name)
        f = self.rho * eta.diff(t_f, 2) + sp.exp(k1) * eta.diff(t_f) + sp.exp(k2) * eta
        self.add_equation("f", f)

        self._apply_detach()

Volterra

Bases: PDE

A second kind of volterra integral equation with Gaussian quadrature algorithm.

\[ x(t) - f(t)=\int_a^t K(t, s) x(s) d s \]

Volterra integral equation

Gaussian quadrature

Parameters:

Name Type Description Default
bound float

Lower bound a for Volterra integral equation.

required
num_points int

Sampled points in integral interval.

required
quad_deg int

Number of quadrature.

required
kernel_func Callable

Kernel func K(t,s).

required
func Callable

x(t) - f(t) in Volterra integral equation.

required

Examples:

>>> import ppsci
>>> import numpy as np
>>> vol_eq = ppsci.equation.Volterra(
...     0, 12, 20, lambda t, s: np.exp(s - t), lambda out: out["u"],
... )
Source code in ppsci/equation/ide/volterra.py
class Volterra(PDE):
    r"""A second kind of volterra integral equation with Gaussian quadrature algorithm.

    $$
    x(t) - f(t)=\int_a^t K(t, s) x(s) d s
    $$

    [Volterra integral equation](https://en.wikipedia.org/wiki/Volterra_integral_equation)

    [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval)

    Args:
        bound (float): Lower bound `a` for Volterra integral equation.
        num_points (int): Sampled points in integral interval.
        quad_deg (int): Number of quadrature.
        kernel_func (Callable): Kernel func `K(t,s)`.
        func (Callable): `x(t) - f(t)` in Volterra integral equation.

    Examples:
        >>> import ppsci
        >>> import numpy as np
        >>> vol_eq = ppsci.equation.Volterra(
        ...     0, 12, 20, lambda t, s: np.exp(s - t), lambda out: out["u"],
        ... )
    """

    dtype = paddle.get_default_dtype()

    def __init__(
        self,
        bound: float,
        num_points: int,
        quad_deg: int,
        kernel_func: Callable,
        func: Callable,
    ):
        super().__init__()
        self.bound = bound
        self.num_points = num_points
        self.quad_deg = quad_deg
        self.kernel_func = kernel_func
        self.func = func

        self.quad_x, self.quad_w = np.polynomial.legendre.leggauss(quad_deg)
        self.quad_x = self.quad_x.astype(Volterra.dtype).reshape([-1, 1])  # [Q, 1]
        self.quad_x = paddle.to_tensor(self.quad_x)  # [Q, 1]

        self.quad_w = self.quad_w.astype(Volterra.dtype)  # [Q, ]

        def compute_volterra_func(out):
            x, u = out["x"], out["u"]
            lhs = self.func(out)

            int_mat = paddle.to_tensor(self._get_int_matrix(x), stop_gradient=False)
            rhs = paddle.mm(int_mat, u)  # (N, 1)

            volterra = lhs[: len(rhs)] - rhs
            return volterra

        self.add_equation("volterra", compute_volterra_func)

    def get_quad_points(self, t: paddle.Tensor) -> paddle.Tensor:
        """Scale and transform quad_x from [-1, 1] to range [a, b].

        reference: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval

        Args:
            t (paddle.Tensor): Tensor array of upper bounds 't' for integral.

        Returns:
            paddle.Tensor: Transformed points in desired range with shape of [N, Q].
        """
        a, b = self.bound, t
        return ((b - a) / 2) @ self.quad_x.T + (b + a) / 2

    def _get_quad_weights(self, t: float) -> np.ndarray:
        """Scale weights to range according to given t and lower bound of integral.

        reference: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval

        Args:
            t (float): Array of upper bound 't' for integral.

        Returns:
            np.ndarray: Transformed weights in desired range with shape of [Q, ].
        """
        a, b = self.bound, t
        return (b - a) / 2 * self.quad_w

    def _get_int_matrix(self, x: np.ndarray) -> np.ndarray:
        int_mat = np.zeros(
            (self.num_points, self.num_points + (self.num_points * self.quad_deg)),
            dtype=Volterra.dtype,
        )
        for i in range(self.num_points):
            xi = float(x[i])
            beg = self.num_points + self.quad_deg * i
            end = self.num_points + self.quad_deg * (i + 1)
            K = np.ravel(
                self.kernel_func(np.full((self.quad_deg, 1), xi), x[beg:end].numpy())
            )
            int_mat[i, beg:end] = self._get_quad_weights(xi) * K
        return int_mat
get_quad_points(t)

Scale and transform quad_x from [-1, 1] to range [a, b].

reference: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval

Parameters:

Name Type Description Default
t Tensor

Tensor array of upper bounds 't' for integral.

required

Returns:

Type Description
Tensor

paddle.Tensor: Transformed points in desired range with shape of [N, Q].

Source code in ppsci/equation/ide/volterra.py
def get_quad_points(self, t: paddle.Tensor) -> paddle.Tensor:
    """Scale and transform quad_x from [-1, 1] to range [a, b].

    reference: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval

    Args:
        t (paddle.Tensor): Tensor array of upper bounds 't' for integral.

    Returns:
        paddle.Tensor: Transformed points in desired range with shape of [N, Q].
    """
    a, b = self.bound, t
    return ((b - a) / 2) @ self.quad_x.T + (b + a) / 2

NLSMB

Bases: PDE

Class for nonlinear Schrodinger-Maxwell-Bloch equation.

\[ \begin{cases} \dfrac{\partial E}{\partial x} = i \alpha_1 \dfrac{\partial^2 E}{\partial t ^2} - i \alpha_2 |E|^2 E+2 p \\ \dfrac{\partial p}{\partial t} = 2 i \omega_0 p+2 E \eta \\ \dfrac{\partial \eta}{\partial t} = -(E p^* + E^* p) \end{cases} \]

Parameters:

Name Type Description Default
alpha_1 Union[float, str]

Group velocity dispersion.

required
alpha_2 Union[float, str]

Kerr nonlinearity.

required
omega_0 Union[float, str]

The offset of resonance frequency.

required
time bool

Whether the equation is time-dependent.

required
detach_keys Optional[Tuple[str, ...]]

Keys used for detach during computing. Defaults to None.

None

Examples:

>>> import ppsci
>>> pde = ppsci.equation.NLSMB(0.5, -1.0, 0.5, True)
Source code in ppsci/equation/pde/nls_m_b.py
class NLSMB(base.PDE):
    r"""Class for nonlinear Schrodinger-Maxwell-Bloch equation.

    $$
    \begin{cases}
        \dfrac{\partial E}{\partial x} = i \alpha_1 \dfrac{\partial^2 E}{\partial t ^2} - i \alpha_2 |E|^2 E+2 p \\
        \dfrac{\partial p}{\partial t} = 2 i \omega_0 p+2 E \eta \\
        \dfrac{\partial \eta}{\partial t} = -(E p^* + E^* p)
    \end{cases}
    $$

    Args:
        alpha_1 (Union[float, str]): Group velocity dispersion.
        alpha_2 (Union[float, str]): Kerr nonlinearity.
        omega_0 (Union[float, str]): The offset of resonance frequency.
        time (bool): Whether the equation is time-dependent.
        detach_keys (Optional[Tuple[str, ...]]): Keys used for detach during computing.
            Defaults to None.

    Examples:
        >>> import ppsci
        >>> pde = ppsci.equation.NLSMB(0.5, -1.0, 0.5, True)
    """

    def __init__(
        self,
        alpha_1: Union[float, str],
        alpha_2: Union[float, str],
        omega_0: Union[float, str],
        time: bool,
        detach_keys: Optional[Tuple[str, ...]] = None,
    ):
        super().__init__()
        self.detach_keys = detach_keys
        self.time = time

        t, x = self.create_symbols("t x")
        invars = (x,)
        if time:
            invars = (t,) + invars

        self.alpha_1 = alpha_1
        self.alpha_2 = alpha_2
        self.omega_0 = omega_0

        Eu = self.create_function("Eu", invars)
        Ev = self.create_function("Ev", invars)
        pu = self.create_function("pu", invars)
        pv = self.create_function("pv", invars)
        eta = self.create_function("eta", invars)

        pu_t = pu.diff(t)
        pv_t = pv.diff(t)
        eta_t = eta.diff(t)

        Eu_x = Eu.diff(x)
        Ev_x = Ev.diff(x)

        Eu_tt = Eu.diff(t).diff(t)
        Ev_tt = Ev.diff(t).diff(t)

        Schrodinger_1 = (
            alpha_1 * Eu_tt - alpha_2 * Eu * (Eu**2 + Ev**2) + 2 * pv - Ev_x
        )
        Schrodinger_2 = (
            alpha_1 * Ev_tt - alpha_2 * Ev * (Eu**2 + Ev**2) - 2 * pu + Eu_x
        )
        Maxwell_1 = 2 * Ev * eta - pv_t + 2 * pu * omega_0
        Maxwell_2 = -2 * Eu * eta + pu_t + 2 * pv * omega_0
        Bloch = 2 * pv * Ev + 2 * pu * Eu + eta_t

        self.add_equation("Schrodinger_1", Schrodinger_1)
        self.add_equation("Schrodinger_2", Schrodinger_2)
        self.add_equation("Maxwell_1", Maxwell_1)
        self.add_equation("Maxwell_2", Maxwell_2)
        self.add_equation("Bloch", Bloch)

        self._apply_detach()