跳转至

Equation(方程) 模块

ppsci.equation

PDE

Base class for Partial Differential Equation

Source code in ppsci/equation/pde/base.py 
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
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 
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
@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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
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 
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
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()

FractionalPoisson

Bases: PDE

TODO: refine this docstring 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])
Source code in ppsci/equation/fpde/fractional_poisson.py 
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
class FractionalPoisson(PDE):
    r"""

    TODO: refine this docstring
    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

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 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
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()

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 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
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 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
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 
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
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 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
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 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
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 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
86
87
88
89
90
91
92
93
94
95
96
97
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 
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
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()