From 6267eda37fca2715342bbd2dc75a2237db14e9aa Mon Sep 17 00:00:00 2001 From: weipengOO98 Date: Thu, 29 Jun 2023 10:46:16 +0800 Subject: [PATCH] no message --- README.md | 4 +- docs/conf.py | 4 +- docs/user/get_started/10_deepritz.md | 123 ++++++++++ .../get_started/9_navier_stokes_equation.md | 227 ++++++++++++++++++ docs/user/get_started/tutorial.rst | 4 + examples/deepritz/deepritz.py | 75 ++++++ examples/euler_beam/euler_beam.py | 3 +- examples/euler_beam/euler_beam_lbfgs.py | 79 ++++++ examples/ns_steady/NS_steady.py | 194 +++++++++++++++ examples/ns_unsteady/NS_unsteady.py | 188 +++++++++++++++ idrlnet/solver.py | 37 ++- idrlnet/variable.py | 17 ++ requirements.txt | 4 +- 13 files changed, 942 insertions(+), 17 deletions(-) create mode 100644 docs/user/get_started/10_deepritz.md create mode 100644 docs/user/get_started/9_navier_stokes_equation.md create mode 100644 examples/deepritz/deepritz.py create mode 100644 examples/euler_beam/euler_beam_lbfgs.py create mode 100644 examples/ns_steady/NS_steady.py create mode 100644 examples/ns_unsteady/NS_unsteady.py diff --git a/README.md b/README.md index c4f3ea0..208dbc9 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ pip install idrlnet ### From Source ``` -git clone https://osredm.com/idrl/idrlnet.git +git clone https://github.com/idrl-lab/idrlnet cd idrlnet pip install -e . ``` @@ -100,7 +100,7 @@ It is also easy to customize IDRLnet to meet new demands. First off, thanks for taking the time to contribute! -- **Reporting bugs.** To report a bug, simply open an issue(疑修) in the osredm "Issues" section. +- **Reporting bugs.** To report a bug, simply open an issue in the GitHub "Issues" section. - **Suggesting enhancements.** To submit an enhancement suggestion for IDRLnet, including completely new features and minor improvements to existing functionality, let us know by opening an issue. diff --git a/docs/conf.py b/docs/conf.py index 06309b2..08c21f9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -18,11 +18,11 @@ sys.path.insert(0, os.path.abspath("..")) # -- Project information ----------------------------------------------------- project = "idrlnet" -copyright = "2021, IDRL" +copyright = "2023, IDRL" author = "IDRL" # The full version, including alpha/beta/rc tags -release = "0.1.0" +release = "0.0.2-rc3" # -- General configuration --------------------------------------------------- diff --git a/docs/user/get_started/10_deepritz.md b/docs/user/get_started/10_deepritz.md new file mode 100644 index 0000000..dae6d80 --- /dev/null +++ b/docs/user/get_started/10_deepritz.md @@ -0,0 +1,123 @@ +# Deepritz + +This section repeats the Deepritz method presented by [Weinan E and Bing Yu](https://link.springer.com/article/10.1007/s40304-018-0127-z). + +Consider the 2d Poisson's equation such as the following: + +$$ +\begin{equation} +\begin{aligned} +-\Delta u=f, & \text { in } \Omega \\ +u=0, & \text { on } \partial \Omega +\end{aligned} +\end{equation} +$$ + +Based on the scattering theorem, its weak form is that both sides are multiplied by$ v \in H_0^1$(which can be interpreted as some function bounded by 0),to get + +$$ +\int f v=-\int v \Delta u=(\nabla u, \nabla v) +$$ + +The above equation holds for any $v \in H_0^1$. The bilinear part of the right-hand side of the equation with respect to $u,v$ is symmetric and yields the bilinear term: + +$$ +a(u, v)=\int \nabla u \cdot \nabla v +$$ + +By the Poincaré inequality, the $a(\cdot, \cdot)$ is a positive definite operator. By positive definite, we mean that there exists $\alpha >0$, such that + +$$ +a(u, u) \geq \alpha\|u\|^2, \quad \forall u \in H_0^1 +$$ + +The remaining term is a linear generalization of $v$, which is $l(v)$, which yields the equation: + +$$ +a(u, v) = l(v) +$$ + +For this equation, by discretizing $u,v$ in the same finite dimensional subspace, we can obtain a symmetric positive definite system of equations, which is the family of Galerkin methods, or we can transform it into a polarization problem to solve it. + +To find $u$ satisfies + +$$ +a(u, v) = l(v), \quad \forall v \in H_0^1 +$$ + +For a symmetric positive definite $a$ , which is equivalent to solving the variational minimization problem, that is, finding $u$, such that holds, where + +$$ +J(u) = \frac{1}{2} a(u, u) - l(u) +$$ + +Specifically + +$$ +\min _{u \in H_0^1} J(u)=\frac{1}{2} \int\|\nabla u\|_2^2-\int f v +$$ + +The DeepRitz method is similar to the PINN approach, replacing the neural network with u, and after sampling the region, just solve it with a solver like Adam. Written as + +$$ +\begin{equation} +\min _{\left.\hat{u}\right|_{\partial \Omega}=0} \hat{J}(\hat{u})=\frac{1}{2} \frac{S_{\Omega}}{N_{\Omega}} \sum\left\|\nabla \hat{u}\left(x_i, y_i\right)\right\|_2^2-\frac{S_{\Omega}}{N_{\partial \Omega}} \sum f\left(x_i, y_i\right) \hat{u}\left(x_i, y_i\right) +\end{equation} +$$ + +Note that the original $u \in H_0^1$, which is zero on the boundary, is transformed into an unconstrained problem by adding the penalty function term: + +$$ +\begin{equation} +\begin{gathered} +\min \hat{J}(\hat{u})=\frac{1}{2} \frac{S_{\Omega}}{N_{\Omega}} \sum\left\|\nabla \hat{u}\left(x_i, y_i\right)\right\|_2^2-\frac{S_{\Omega}}{N_{\Omega}} \sum f\left(x_i, y_i\right) \hat{u}\left(x_i, y_i\right)+\beta \frac{S_{\partial \Omega}}{N_{\partial \Omega}} \\ +\sum \hat{u}^2\left(x_i, y_i\right) +\end{gathered} +\end{equation} +$$ + +Consider the 2d Poisson's equation defined on $\Omega=[-1,1]\times[-1,1]$, which satisfies $f=2 \pi^2 \sin (\pi x) \sin (\pi y)$. + +### Define Sampling Methods and Constraints + +For the problem, boundary condition and PDE constraint are presented and use the Identity loss. + +```python +@sc.datanode(sigma=1000.0) +class Boundary(sc.SampleDomain): + def __init__(self): + self.points = geo.sample_boundary(100,) + self.constraints = {"u": 0.} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + + +@sc.datanode(loss_fn="Identity") +class Interior(sc.SampleDomain): + def __init__(self): + self.points = geo.sample_interior(1000) + self.constraints = {"integral_dxdy": 0,} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints +``` + +### Define Neural Networks and PDEs + +In the PDE definition section, based on the DeepRitz method we add two types of PDE nodes: + +```python +def f(x, y): + return 2 * sp.pi ** 2 * sp.sin(sp.pi * x) * sp.sin(sp.pi * y) + +dx_exp = sc.ExpressionNode( + expression=0.5*(u.diff(x) ** 2 + u.diff(y) ** 2) - u * f(x, y), name="dxdy" +) +net = sc.get_net_node(inputs=("x", "y"), outputs=("u",), name="net", arch=sc.Arch.mlp) + +integral = sc.ICNode("dxdy", dim=2, time=False) +``` + +The result is shown as follows: +![deepritz](https://github.com/xiangzixuebit/picture/raw/3d73005f3642f10400975659479e856fb99f6518/deepritz.png) diff --git a/docs/user/get_started/9_navier_stokes_equation.md b/docs/user/get_started/9_navier_stokes_equation.md new file mode 100644 index 0000000..284af9f --- /dev/null +++ b/docs/user/get_started/9_navier_stokes_equation.md @@ -0,0 +1,227 @@ +# Navier-Stokes equations + +This section repeats the Robust PINN method presented by [Peng et.al](https://deepai.org/publication/robust-regression-with-highly-corrupted-data-via-physics-informed-neural-networks). + +## Steady 2D NS equations + +The prototype problem of incompressible flow past a circular cylinder is considered. +![image](https://github.com/xiangzixuebit/picture/raw/3d73005f3642f10400975659479e856fb99f6518/NS1.png) + +The velocity vector is set to zero at all walls and the pressure is set to p = 0 at the outlet. The fluid density is taken as $\rho = 1kg/m^3$ and the dynamic viscosity is taken as $\mu = 2 · 10^{−2}kg/m^3$ . The velocity profile on the inlet is set as $u(0, y)=4 \frac{U_M}{H^2}(H-y) y$ with $U_M = 1m/s$ and $H = 0.41m$. + +The two-dimensional steady-state Navier-Stokes equation is equivalently transformed into the following equations: + +$$ +\begin{equation} +\begin{aligned} +\sigma^{11} &=-p+2 \mu u_x \\ +\sigma^{22} &=-p+2 \mu v_y \\ +\sigma^{12} &=\mu\left(u_y+v_x\right) \\ +p &=-\frac{1}{2}\left(\sigma^{11}+\sigma^{22}\right) \\ +\left(u u_x+v u_y\right) &=\mu\left(\sigma_x^{11}+\sigma_y^{12}\right) \\ +\left(u v_x+v v_y\right) &=\mu\left(\sigma_x^{12}+\sigma_y^{22}\right) +\end{aligned} +\end{equation} +$$ + +We construct a neural network with six outputs to satisfy the PDE constraints above: + +$$ +\begin{equation} +u, v, p, \sigma^{11}, \sigma^{12}, \sigma^{22}=net(x, y) +\end{equation} +$$ + +### Define Symbols and Geometric Objects + +For the 2d problem, we define two coordinate symbols`x`and`y`, six variables$ u, v, p, \sigma^{11}, \sigma^{12}, \sigma^{22}$ are defined. + +The geometry object is a simple rectangle and circle with the operator `-`. + +```python +x = Symbol('x') +y = Symbol('y') +rec = sc.Rectangle((0., 0.), (1.1, 0.41)) +cir = sc.Circle((0.2, 0.2), 0.05) +geo = rec - cir +u = sp.Function('u')(x, y) +v = sp.Function('v')(x, y) +p = sp.Function('p')(x, y) +s11 = sp.Function('s11')(x, y) +s22 = sp.Function('s22')(x, y) +s12 = sp.Function('s12')(x, y) +``` + +### Define Sampling Methods and Constraints + +For the problem, three boundary conditions , PDE constraint and external data are presented. We use the robust-PINN model inspired by the traditional LAD (Least Absolute Derivation) approach, where the L1 loss replaces the squared L2 data loss. + +```python +@sc.datanode +class Inlet(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = rec.sample_boundary(1000, sieve=(sp.Eq(x, 0.))) + constraints = sc.Variables({'u': 4 * (0.41 - y) * y / (0.41 * 0.41)}) + return points, constraints + +@sc.datanode +class Outlet(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_boundary(1000, sieve=(sp.Eq(x, 1.1))) + constraints = sc.Variables({'p': 0.}) + return points, constraints + +@sc.datanode +class Wall(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_boundary(1000, sieve=((x > 0.) & (x < 1.1))) + #print("points3", points) + constraints = sc.Variables({'u': 0., 'v': 0.}) + return points, constraints + +@sc.datanode(name='NS_external') +class Interior_domain(sc.SampleDomain): + def __init__(self): + self.density = 2000 + + def sampling(self, *args, **kwargs): + points = geo.sample_interior(2000) + constraints = {'f_s11': 0., 'f_s22': 0., 'f_s12': 0., 'f_u': 0., 'f_v': 0., 'f_p': 0.} + return points, constraints + +@sc.datanode(name='NS_domain', loss_fn='L1') +class NSExternal(sc.SampleDomain): + def __init__(self): + points = pd.read_csv('NSexternel_sample.csv') + self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns} + self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints +``` + +### Define Neural Networks and PDEs + +In the PDE definition part, we add these PDE nodes: + +```python +net = sc.MLP([2, 40, 40, 40, 40, 40, 40, 40, 40, 6], activation=sc.Activation.tanh) +net = sc.get_net_node(inputs=('x', 'y'), outputs=('u', 'v', 'p', 's11', 's22', 's12'), name='net', arch=sc.Arch.mlp) +pde1 = sc.ExpressionNode(name='f_s11', expression=-p + 2 * nu * u.diff(x) - s11) +pde2 = sc.ExpressionNode(name='f_s22', expression=-p + 2 * nu * v.diff(y) - s22) +pde3 = sc.ExpressionNode(name='f_s12', expression=nu * (u.diff(y) + v.diff(x)) - s12) +pde4 = sc.ExpressionNode(name='f_u', expression=u * u.diff(x) + v * u.diff(y) - nu * (s11.diff(x) + s12.diff(y))) +pde5 = sc.ExpressionNode(name='f_v', expression=u * v.diff(x) + v * v.diff(y) - nu * (s12.diff(x) + s22.diff(y))) +pde6 = sc.ExpressionNode(name='f_p', expression=p + (s11 + s22) / 2) +``` + +### Define A Solver + +Direct use of Adam optimization is less effective, so the LBFGS optimization method or a combination of both (Adam+LBFGS) is used for training: + +```python +s = sc.Solver(sample_domains=(Inlet(), Outlet(), Wall(), Interior_domain(), NSExternal()), + netnodes=[net], + init_network_dirs=['network_dir_adam'], + pdes=[pde1, pde2, pde3, pde4, pde5, pde6], + max_iter=300, + opt_config = dict(optimizer='LBFGS', lr=1) + ) +``` + +The result is shown as follows: +![image](https://github.com/xiangzixuebit/picture/raw/3d73005f3642f10400975659479e856fb99f6518/NS11.png) + +## Unsteady 2D N-S equations with unknown parameters + +A two-dimensional incompressible flow and dynamic vortex shedding past a circular cylinder in a steady-state are numerically simulated. Respectively, the Reynolds number of the incompressible flow is $Re = 100$. The kinematic viscosity of the fluid is $\nu = 0.01$. The cylinder diameter D is 1. The simulation domain size is +$[-15,25] × [-8,8]$. The computational domain is much smaller: $[1,8] × [-2,2]× [0,20]$. + +![image](https://github.com/xiangzixuebit/picture/raw/3d73005f3642f10400975659479e856fb99f6518/NS2.png) + +$$ +\begin{equation} +\begin{aligned} +&u_t+\lambda_1\left(u u_x+v u_y\right)=-p_x+\lambda_2\left(u_{x x}+u_{y y}\right) \\ +&v_t+\lambda_1\left(u v_x+v v_y\right)=-p_y+\lambda_2\left(v_{x x}+v_{y y}\right) +\end{aligned} +\end{equation} +$$ + +where $\lambda_1$ and $\lambda_2$ are two unknown parameters to be recovered. We make the assumption that $u=\psi_y, \quad v=-\psi_x$ + +for some stream function $\psi(x, y)$. Under this assumption, the continuity equation will be automatically satisfied. The following architecture is used in this example, + +$$ +\begin{equation} +\psi, p=net\left(t, x, y, \lambda_1, \lambda_2\right) +\end{equation} +$$ + +### Define Symbols and Geometric Objects + +We define three coordinate symbols `x`, `y` and `t`, three variables $u,v,p$ are defined. + +```python +x = Symbol('x') +y = Symbol('y') +t = Symbol('t') +geo = sc.Rectangle((1., -2.), (8., 2.)) +u = sp.Function('u')(x, y, t) +v = sp.Function('v')(x, y, t) +p = sp.Function('p')(x, y, t) +time_range = {t: (0, 20)} +``` + +### Define Sampling Methods and Constraints + +This example has only two equation constraints, while the former has six equation constraints. We also use the LAD-PINN model. Then the PDE constrained optimization model is formulated as: + +$$ +\min _{\theta, \lambda} \frac{1}{\# \mathbf{D}_u} \sum_{\left(t_i, x_i, u_i\right) \in \mathbf{D}_u}\left|u_i-u_\theta\left(t_i, x_i ; \lambda\right)\right|+\omega \cdot L_{p d e} . +$$ + +```python +@sc.datanode(name='NS_domain', loss_fn='L1') +class NSExternal(sc.SampleDomain): + def __init__(self): + points = pd.read_csv('NSexternel_sample.csv') + self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns} + self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + +@sc.datanode(name='NS_external') +class NSEq(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_interior(density=2000, param_ranges=time_range) + constraints = {'continuity': 0, 'momentum_x': 0, 'momentum_y': 0} + return points, constraints +``` + +### Define Neural Networks and PDEs + +IDRLnet defines a network node to represent the unknown Parameters. + +```python +net = sc.MLP([3, 20, 20, 20, 20, 20, 20, 20, 20, 3], activation=sc.Activation.tanh) +net = sc.get_net_node(inputs=('x', 'y', 't'), outputs=('u', 'v', 'p'), name='net', arch=sc.Arch.mlp) +var_nr = sc.get_net_node(inputs=('x', 'y'), outputs=('nu', 'rho'), arch=sc.Arch.single_var) +pde = sc.NavierStokesNode(nu='nu', rho='rho', dim=2, time=True, u='u', v='v', p='p') +``` + +### Define A Solver + +Two nodes trained together + +```python +s = sc.Solver(sample_domains=(NSExternal(), NSEq()), + netnodes=[net, var_nr], + pdes=[pde], + network_dir='network_dir', + max_iter=10000) +``` + +Finally, the real velocity field and pressure field at t=10s are compared with the predicted results: +![image](https://github.com/xiangzixuebit/picture/raw/3d73005f3642f10400975659479e856fb99f6518/NS22.png) diff --git a/docs/user/get_started/tutorial.rst b/docs/user/get_started/tutorial.rst index 2acf4cf..02739eb 100644 --- a/docs/user/get_started/tutorial.rst +++ b/docs/user/get_started/tutorial.rst @@ -14,6 +14,8 @@ To make full use of IDRLnet. We strongly suggest following the following example 6. :ref:`Parameterized poisson equation `. The example introduces how to train a surrogate with parameters. 7. :ref:`Variational Minimization `. The example introduces how to solve variational minimization problems. 8. :ref:`Volterra integral differential equation `. The example introduces the way to solve IDEs. +9. :ref:`Navier-Stokes equation `. The example introduces how to use the LBFGS optimizer. +10. :ref:`Deepritz method `. The example introduces the way to solve PDEs with the Deepritz method. @@ -28,3 +30,5 @@ To make full use of IDRLnet. We strongly suggest following the following example 6_parameterized_poisson 7_minimal_surface 8_volterra_ide + 9_navier_stokes_equation + 10_deepritz diff --git a/examples/deepritz/deepritz.py b/examples/deepritz/deepritz.py new file mode 100644 index 0000000..23be1f3 --- /dev/null +++ b/examples/deepritz/deepritz.py @@ -0,0 +1,75 @@ +import matplotlib.pyplot as plt +import numpy as np +import sympy as sp +import matplotlib.tri as tri +import idrlnet.shortcut as sc + +x, y = sp.symbols("x y") +u = sp.Function("u")(x, y) +geo = sc.Rectangle((-1, -1), (1., 1.)) + + +@sc.datanode(sigma=1000.0) +class Boundary(sc.SampleDomain): + def __init__(self): + self.points = geo.sample_boundary(100,) + self.constraints = {"u": 0.} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + + +@sc.datanode(loss_fn="Identity") +class Interior(sc.SampleDomain): + def __init__(self): + self.points = geo.sample_interior(1000) + self.constraints = {"integral_dxdy": 0,} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + + +def f(x, y): + return 2 * sp.pi ** 2 * sp.sin(sp.pi * x) * sp.sin(sp.pi * y) + + +dx_exp = sc.ExpressionNode( + expression=0.5*(u.diff(x) ** 2 + u.diff(y) ** 2) - u * f(x, y), name="dxdy" +) +net = sc.get_net_node(inputs=("x", "y"), outputs=("u",), name="net", arch=sc.Arch.mlp) + +integral = sc.ICNode("dxdy", dim=2, time=False) + +s = sc.Solver( + sample_domains=(Boundary(), Interior()), + netnodes=[net], + pdes=[ + dx_exp, + integral, + ], + max_iter=10000, +) +s.solve() +coord = s.infer_step({"Interior": ["x", "y", "u"]}) +num_x = coord["Interior"]["x"].cpu().detach().numpy().ravel() +num_y = coord["Interior"]["y"].cpu().detach().numpy().ravel() +num_Up = coord["Interior"]["u"].cpu().detach().numpy().ravel() + +# Ground truth +num_U = np.sin(np.pi*num_x)*np.sin(np.pi*num_y) + +fig, ax = plt.subplots(1, 3, figsize=(10, 3)) +triang_total = tri.Triangulation(num_x, num_y) +ax[0].tricontourf(triang_total, num_Up, 100, cmap="bwr", vmin=-1, vmax=1) +ax[0].axis("off") +ax[0].set_title("prediction") +ax[1].tricontourf(triang_total, num_U, 100, cmap="bwr", vmin=-1, vmax=1) +ax[1].axis("off") +ax[1].set_title("ground truth") +ax[2].tricontourf( + triang_total, np.abs(num_U - num_Up), 100, cmap="bwr", vmin=0, vmax=0.5 +) +ax[2].axis("off") +ax[2].set_title("absolute error") + +plt.savefig("deepritz.png", dpi=300, bbox_inches="tight") diff --git a/examples/euler_beam/euler_beam.py b/examples/euler_beam/euler_beam.py index 8962298..7b1b75e 100644 --- a/examples/euler_beam/euler_beam.py +++ b/examples/euler_beam/euler_beam.py @@ -63,8 +63,7 @@ solver = sc.Solver( ), netnodes=[net], pdes=[pde1, pde2, pde3, pde4], - max_iter=2000, -) + max_iter=2000) solver.solve() diff --git a/examples/euler_beam/euler_beam_lbfgs.py b/examples/euler_beam/euler_beam_lbfgs.py new file mode 100644 index 0000000..4c54684 --- /dev/null +++ b/examples/euler_beam/euler_beam_lbfgs.py @@ -0,0 +1,79 @@ +import matplotlib.pyplot as plt +import sympy as sp +import numpy as np +import idrlnet.shortcut as sc + +x = sp.symbols('x') +Line = sc.Line1D(0, 1) +y = sp.Function('y')(x) + + +@sc.datanode(name='interior') +class Interior(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return Line.sample_interior(1000), {'dddd_y': 0} + + +@sc.datanode(name='left_boundary1') +class LeftBoundary1(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return Line.sample_boundary(100, sieve=(sp.Eq(x, 0))), {'y': 0} + + +@sc.datanode(name='left_boundary2') +class LeftBoundary2(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return Line.sample_boundary(100, sieve=(sp.Eq(x, 0))), {'d_y': 0} + + +@sc.datanode(name='right_boundary1') +class RightBoundary1(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return Line.sample_boundary(100, sieve=(sp.Eq(x, 1))), {'dd_y': 0} + + +@sc.datanode(name='right_boundary2') +class RightBoundary2(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return Line.sample_boundary(100, sieve=(sp.Eq(x, 1))), {'ddd_y': 0} + + +@sc.datanode(name='infer') +class Infer(sc.SampleDomain): + def sampling(self, *args, **kwargs): + return {'x': np.linspace(0, 1, 1000).reshape(-1, 1)}, {} + + +net = sc.get_net_node(inputs=('x',), outputs=('y',), name='net', arch=sc.Arch.mlp) + +pde1 = sc.ExpressionNode(name='dddd_y', expression=y.diff(x).diff(x).diff(x).diff(x) + 1) +pde2 = sc.ExpressionNode(name='d_y', expression=y.diff(x)) +pde3 = sc.ExpressionNode(name='dd_y', expression=y.diff(x).diff(x)) +pde4 = sc.ExpressionNode(name='ddd_y', expression=y.diff(x).diff(x).diff(x)) + +solver = sc.Solver( + sample_domains=(Interior(), LeftBoundary1(), LeftBoundary2(), RightBoundary1(), RightBoundary2()), + netnodes=[net], + pdes=[pde1, pde2, pde3, pde4], + max_iter=200, + opt_config=dict(optimizer='LBFGS', lr=1)) +solver.solve() + + +# inference +def exact(x): + return -(x ** 4) / 24 + x ** 3 / 6 - x ** 2 / 4 + + +solver.sample_domains = (Infer(),) +points = solver.infer_step({'infer': ['x', 'y']}) +xs = points['infer']['x'].detach().cpu().numpy().ravel() +y_pred = points['infer']['y'].detach().cpu().numpy().ravel() +plt.plot(xs, y_pred, label='Pred') +y_exact = exact(xs) +plt.plot(xs, y_exact, label='Exact', linestyle='--') +plt.legend() +plt.xlabel('x') +plt.ylabel('w') +plt.savefig('Euler_beam_LBFGS.png', dpi=300, bbox_inches='tight') +plt.show() diff --git a/examples/ns_steady/NS_steady.py b/examples/ns_steady/NS_steady.py new file mode 100644 index 0000000..0e1e4da --- /dev/null +++ b/examples/ns_steady/NS_steady.py @@ -0,0 +1,194 @@ +import matplotlib.pyplot as plt +import sympy as sp +import numpy as np +import idrlnet.shortcut as sc +from sympy import Symbol, sin +import pandas as pd +import torch +import matplotlib.tri as tri + +x = Symbol('x') +y = Symbol('y') +rec = sc.Rectangle((0., 0.), (1.1, 0.41)) +cir = sc.Circle((0.2, 0.2), 0.05) +geo = rec - cir +u = sp.Function('u')(x, y) +v = sp.Function('v')(x, y) +p = sp.Function('p')(x, y) +s11 = sp.Function('s11')(x, y) +s22 = sp.Function('s22')(x, y) +s12 = sp.Function('s12')(x, y) +nu=0.02 +rho=1 + +@sc.datanode +class Inlet(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = rec.sample_boundary(1000, sieve=(sp.Eq(x, 0.))) + constraints = sc.Variables({'u': 4 * (0.41 - y) * y / (0.41 * 0.41)}) + return points, constraints + +@sc.datanode +class Outlet(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_boundary(1000, sieve=(sp.Eq(x, 1.1))) + constraints = sc.Variables({'p': 0.}) + return points, constraints + +@sc.datanode +class Wall(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_boundary(1000, sieve=((x > 0.) & (x < 1.1))) + #print("points3", points) + constraints = sc.Variables({'u': 0., 'v': 0.}) + return points, constraints + +@sc.datanode(name='NS_external') +class Interior_domain(sc.SampleDomain): + def __init__(self): + self.density = 2000 + + def sampling(self, *args, **kwargs): + points = geo.sample_interior(2000) + constraints = {'f_s11': 0., 'f_s22': 0., 'f_s12': 0., 'f_u': 0., 'f_v': 0., 'f_p': 0.} + return points, constraints + +@sc.datanode(name='NS_domain', loss_fn='L1') +class NSExternal(sc.SampleDomain): + def __init__(self): + points = pd.read_csv('NSexternel_sample.csv') + self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns} + self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + +net = sc.MLP([2, 40, 40, 40, 40, 40, 40, 40, 40, 6], activation=sc.Activation.tanh) +net = sc.get_net_node(inputs=('x', 'y'), outputs=('u', 'v', 'p', 's11', 's22', 's12'), name='net', arch=sc.Arch.mlp) +pde1 = sc.ExpressionNode(name='f_s11', expression=-p + 2 * nu * u.diff(x) - s11) +pde2 = sc.ExpressionNode(name='f_s22', expression=-p + 2 * nu * v.diff(y) - s22) +pde3 = sc.ExpressionNode(name='f_s12', expression=nu * (u.diff(y) + v.diff(x)) - s12) +pde4 = sc.ExpressionNode(name='f_u', expression=u * u.diff(x) + v * u.diff(y) - nu * (s11.diff(x) + s12.diff(y))) +pde5 = sc.ExpressionNode(name='f_v', expression=u * v.diff(x) + v * v.diff(y) - nu * (s12.diff(x) + s22.diff(y))) +pde6 = sc.ExpressionNode(name='f_p', expression=p + (s11 + s22) / 2) +s = sc.Solver(sample_domains=(Inlet(), Outlet(), Wall(), Interior_domain(), NSExternal()), + netnodes=[net], + init_network_dirs=['network_dir_adam'], + pdes=[pde1, pde2, pde3, pde4, pde5, pde6], + max_iter=300, + opt_config = dict(optimizer='LBFGS', lr=1) + ) +#opt_config = dict(optimizer='LBFGS', lr=1) +#init_network_dirs=['network_dir_lbfgs'], +s.solve() + +points1 = pd.read_csv('NSexternel_test.csv') +points1 = {col: points1[col].to_numpy().reshape(-1, 1) for col in points1.columns} +x_test = torch.tensor(points1['x_test'].astype(np.float32)) +y_test = torch.tensor(points1['y_test'].astype(np.float32)) +u_test = torch.tensor(points1['u_test'].astype(np.float32)) +v_test = torch.tensor(points1['v_test'].astype(np.float32)) +p_test = torch.tensor(points1['p_test'].astype(np.float32)) + +U = s.netnodes[0].net(torch.cat([x_test, y_test], dim=1)) + +num_x = x_test.cpu().detach().numpy().ravel() +num_y = y_test.cpu().detach().numpy().ravel() +num_u = u_test.cpu().detach().numpy().ravel() +num_v = v_test.cpu().detach().numpy().ravel() +num_p = p_test.cpu().detach().numpy().ravel() + +num_up = U[:, 0:1].cpu().detach().numpy().ravel() +num_vp = U[:, 1:2].cpu().detach().numpy().ravel() +num_pp = U[:, 2:3].cpu().detach().numpy().ravel() + + +triang_total = tri.Triangulation(num_x, num_y) + +font2 = {'family': 'Times New Roman', + 'weight': 'normal', + 'size': 15, + } + + + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_u, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $u$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_up, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $u$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_u - num_up, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_u_Adam.png', dpi=300, bbox_inches='tight') +plt.show() + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_v, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $v$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_vp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $v$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_v - num_vp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_v_Adam.png', dpi=300, bbox_inches='tight') +plt.show() + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_p, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $p$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_pp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $p$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_p - num_pp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_p_Adam.png', dpi=300, bbox_inches='tight') +plt.show() \ No newline at end of file diff --git a/examples/ns_unsteady/NS_unsteady.py b/examples/ns_unsteady/NS_unsteady.py new file mode 100644 index 0000000..344b74f --- /dev/null +++ b/examples/ns_unsteady/NS_unsteady.py @@ -0,0 +1,188 @@ +import matplotlib.pyplot as plt +import sympy as sp +import numpy as np +import idrlnet.shortcut as sc +from sympy import Symbol, sin +import pandas as pd +import torch +import matplotlib.tri as tri + +x = Symbol('x') +y = Symbol('y') +t = Symbol('t') +geo = sc.Rectangle((1., -2.), (8., 2.)) +u = sp.Function('u')(x, y, t) +v = sp.Function('v')(x, y, t) +p = sp.Function('p')(x, y, t) +time_range = {t: (0, 20)} +nu=0.01 +rho=1 + +@sc.datanode(name='NS_domain', loss_fn='L1') +class NSExternal(sc.SampleDomain): + def __init__(self): + points = pd.read_csv('NSexternel_sample.csv') + self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns} + self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')} + + def sampling(self, *args, **kwargs): + return self.points, self.constraints + +@sc.datanode(name='NS_external') +class NSEq(sc.SampleDomain): + def sampling(self, *args, **kwargs): + points = geo.sample_interior(density=2000, param_ranges=time_range) + constraints = {'continuity': 0, 'momentum_x': 0, 'momentum_y': 0} + return points, constraints + +net = sc.MLP([3, 20, 20, 20, 20, 20, 20, 20, 20, 3], activation=sc.Activation.tanh) +net = sc.get_net_node(inputs=('x', 'y', 't'), outputs=('u', 'v', 'p'), name='net', arch=sc.Arch.mlp) +#var_nr = sc.get_net_node(inputs=('x', 'y'), outputs=('nu', 'rho'), arch=sc.Arch.single_var) +#pde = sc.NavierStokesNode(nu='nu', rho='rho', dim=2, time=True, u='u', v='v', p='p') +pde = sc.NavierStokesNode(nu=0.01, rho=1.0, dim=2, time=True) +s = sc.Solver(sample_domains=(NSExternal(), NSEq()), + netnodes=[net], + init_network_dirs=['network_dir_adam'], + pdes=[pde], + max_iter=100, + opt_config=dict(optimizer='LBFGS', lr=1) + ) +#opt_config=dict(optimizer='LBFGS', lr=1) +# s = sc.Solver(sample_domains=(NSExternal(), NSEq()), +# netnodes=[net, var_nr], +# pdes=[pde], +# network_dir='network_dir', +# max_iter=10) +s.solve() + + +coord = s.infer_step(domain_attr={'NS_domain': ['x', 'y', 'u', 'v', 'p']}) +num_xd = coord['NS_domain']['x'].cpu().detach().numpy().ravel() +num_yd = coord['NS_domain']['y'].cpu().detach().numpy().ravel() +num_ud = coord['NS_domain']['u'].cpu().detach().numpy().ravel() +num_vd = coord['NS_domain']['v'].cpu().detach().numpy().ravel() +num_pd = coord['NS_domain']['p'].cpu().detach().numpy().ravel() + +# print("true paratmeter rho: {:.4f}".format(rho)) +# predict_rho = var_nr.evaluate(torch.Tensor([[1.0]])).item() +# print("predicted parameter rho: {:.4f}".format(predict_rho)) + + +points1 = pd.read_csv('NSexternel_test.csv') +points1 = {col: points1[col].to_numpy().reshape(-1, 1) for col in points1.columns} +x_test = torch.tensor(points1['x_test'].astype(np.float32)) +y_test = torch.tensor(points1['y_test'].astype(np.float32)) +t_test = torch.tensor(points1['t_test'].astype(np.float32)) +u_test = torch.tensor(points1['u_test'].astype(np.float32)) +v_test = torch.tensor(points1['v_test'].astype(np.float32)) +p_test = torch.tensor(points1['p_test'].astype(np.float32)) + +U = s.netnodes[0].net(torch.cat([x_test, y_test, t_test], dim=1)) + +num_x = x_test.cpu().detach().numpy().ravel() +num_y = y_test.cpu().detach().numpy().ravel() +num_u = u_test.cpu().detach().numpy().ravel() +num_v = v_test.cpu().detach().numpy().ravel() +num_p = p_test.cpu().detach().numpy().ravel() + +num_up = U[:, 0:1].cpu().detach().numpy().ravel() +num_vp = U[:, 1:2].cpu().detach().numpy().ravel() +num_pp = U[:, 2:3].cpu().detach().numpy().ravel() + + +triang_total = tri.Triangulation(num_x, num_y) + +font2 = {'family': 'Times New Roman', + 'weight': 'normal', + 'size': 15, + } + +# fig = plt.figure(figsize=(6, 6)) +# ax = fig.add_subplot(111) +# ax.scatter(num_xi, num_yi, c='b', s=1, label='Domain') +# ax.set_xlabel('$x$', font2) +# ax.set_ylabel('$y$', font2) +# ax.set_title('collocation points', fontsize=18) +# plt.savefig('points.png', dpi=300, bbox_inches='tight') +# plt.show() + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_u, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $u$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_up, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $u$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_u - num_up, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_u_c.png', dpi=300, bbox_inches='tight') +plt.show() + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_v, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $v$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_v, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $v$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_v - num_vp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_v_c.png', dpi=300, bbox_inches='tight') +plt.show() + +fig = plt.figure(figsize=(20, 4)) +ax1 = fig.add_subplot(131) +tcf = ax1.tricontourf(triang_total, num_p, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax1.set_xlabel('$x$', font2) +ax1.set_ylabel('$y$', font2) +ax1.set_title('Exact $p$', fontsize=18) + +ax2 = fig.add_subplot(132) +tcf = ax2.tricontourf(triang_total, num_pp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax2.set_xlabel('$x$', font2) +ax2.set_ylabel('$y$', font2) +ax2.set_title('Predicted $p$', fontsize=18) + +ax3 = fig.add_subplot(133) +tcf = ax3.tricontourf(triang_total, num_p - num_pp, 100, cmap='jet') +tc_bar = plt.colorbar(tcf) +tc_bar.ax.tick_params(labelsize=12) +ax3.set_xlabel('$x$', font2) +ax3.set_ylabel('$y$', font2) +ax3.set_title('Point-wise Error', fontsize=18) +plt.savefig('test_NS_p_c.png', dpi=300, bbox_inches='tight') +plt.show() \ No newline at end of file diff --git a/idrlnet/solver.py b/idrlnet/solver.py index cd9723b..f4f9f2c 100644 --- a/idrlnet/solver.py +++ b/idrlnet/solver.py @@ -312,18 +312,37 @@ class Solver(Notifier, Optimizable): """ self.notify(self, message={Signal.TRAIN_PIPE_START: "defaults"}) for opt in self.optimizers: - opt.zero_grad() + # print('Running optimization with %s'%(self.optimizer_config['optimizer'])) + if self.optimizer_config['optimizer'] == 'LBFGS': + def closure(): + opt.zero_grad() + samples = self.sample_variables_from_domains() + in_var, true_out, lambda_out = self.generate_in_out_dict(samples) + pred_out_sample = self.forward_through_all_graph(in_var, self.outvar_dict_index) + loss = self.compute_loss(in_var, pred_out_sample, true_out, lambda_out) + self.notify(self, message={Signal.BEFORE_BACKWARD: 'defaults'}) + loss.backward() + return loss + opt.step(closure) + + else: + opt.zero_grad() + samples = self.sample_variables_from_domains() + in_var, true_out, lambda_out = self.generate_in_out_dict(samples) + pred_out_sample = self.forward_through_all_graph(in_var, self.outvar_dict_index) + try: + loss = self.compute_loss(in_var, pred_out_sample, true_out, lambda_out) + except RuntimeError: + raise + self.notify(self, message={Signal.BEFORE_BACKWARD: 'defaults'}) + loss.backward() + opt.step() + samples = self.sample_variables_from_domains() in_var, true_out, lambda_out = self.generate_in_out_dict(samples) pred_out_sample = self.forward_through_all_graph(in_var, self.outvar_dict_index) - try: - loss = self.compute_loss(in_var, pred_out_sample, true_out, lambda_out) - except RuntimeError: - raise - self.notify(self, message={Signal.BEFORE_BACKWARD: "defaults"}) - loss.backward() - for opt in self.optimizers: - opt.step() + loss = self.compute_loss(in_var, pred_out_sample, true_out, lambda_out) + self.global_step += 1 for scheduler in self.schedulers: diff --git a/idrlnet/variable.py b/idrlnet/variable.py index afc1755..717de20 100644 --- a/idrlnet/variable.py +++ b/idrlnet/variable.py @@ -21,6 +21,7 @@ class Loss(enum.Enum): L1 = "L1" square = "square" + Identity = "Identity" class LossFunction: @@ -32,6 +33,8 @@ class LossFunction: return LossFunction.weighted_L1_loss(variables, name=name) elif loss_function == Loss.square.name or loss_function == Loss.square: return LossFunction.weighted_square_loss(variables, name=name) + elif loss_function == Loss.Identity.name or loss_function == Loss.Identity: + return LossFunction.weighted_identity_loss(variables, name=name) raise NotImplementedError(f"loss function {loss_function} is not defined!") @staticmethod @@ -62,6 +65,20 @@ class LossFunction: loss += torch.sum((val ** 2) * variables["area"]) return Variables({name: loss}) + @staticmethod + def weighted_identity_loss(variables: "Variables", name: str) -> "Variables": + loss = 0.0 + for key, val in variables.items(): + if key.startswith("lambda_") or key == "area": + continue + elif "lambda_" + key in variables.keys(): + loss += torch.sum( + val * variables["lambda_" + key] * variables["area"] + ) + else: + loss += torch.sum(val * variables["area"]) + return Variables({name: loss}) + class Variables(dict): def __sub__(self, other: "Variables") -> "Variables": diff --git a/requirements.txt b/requirements.txt index b0cbd5a..89e1573 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,5 @@ transforms3d typing -numpy keras h5py pandas @@ -11,12 +10,13 @@ sphinx matplotlib myst_parser sphinx_markdown_parser +numpy==1.21.0 sphinx_rtd_theme==0.5.2 tensorboard==2.4.1 sympy==1.5.1 pyevtk==1.1.1 flask==1.1.2 requests==2.25.0 -torch>=1.7.1 +torch==1.7.1 networkx==2.5.1 protobuf~=3.20 \ No newline at end of file