diff --git a/.github/workflows/export_tutorials.yml b/.github/workflows/export_tutorials.yml new file mode 100644 index 000000000..ac2d10d3b --- /dev/null +++ b/.github/workflows/export_tutorials.yml @@ -0,0 +1,72 @@ +name: "Export Tutorials" + +on: + push: + branches: + - "**" # Run on push on all branches + paths: + - 'tutorials/**/*.ipynb' +jobs: + export_tutorials: + permissions: write-all + runs-on: ubuntu-latest + env: + TUTORIAL_TIMEOUT: 1200s + steps: + - uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: 3.8 + + - name: Install dependencies + run: | + # Dependencies for tutorials + python3 -m pip install --upgrade pip .[tutorial] black[jupyter] + - name: Setup FFmpeg + uses: FedericoCarboni/setup-ffmpeg@v2 + + - id: files + uses: jitterbit/get-changed-files@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + format: space-delimited + + - name: Configure git + run: | + git config user.name "github-actions[bot]" + git config user.email 41898282+github-actions[bot]@users.noreply.github.com + - name: Export tutorials to .py and .html + run: | + set -x + for file in ${{ steps.files.outputs.all }}; do + if [[ $file == *.ipynb ]]; then + filename=$(basename $file) + pyfilename=$(echo ${filename%?????})py + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert $file --to python --output $pyfilename --output-dir=$(dirname $file) + htmlfilename=$(echo ${filename%?????} | sed -e 's/-//g')html + htmldir="docs/source"/$(echo ${file%?????} | sed -e 's/-//g')html + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert --execute $file --to html --output $htmlfilename --output-dir=$htmldir + fi + done + set +x + + - name: Run formatter + run: black tutorials/ + + - uses: benjlevesque/short-sha@v2.1 + id: short-sha + + - name: Remove unwanted files + run: | + rm -rf build/ + - name: Create Pull Request + uses: peter-evans/create-pull-request@v5.0.2 + with: + labels: maintenance + title: Export tutorial changed in ${{ steps.short-sha.outputs.sha }} + branch: export-tutorial-${{ steps.short-sha.outputs.sha }} + base: ${{ github.head_ref }} + commit-message: export tutorials changed in ${{ steps.short-sha.outputs.sha }} + delete-branch: true diff --git a/docs/source/tutorials/tutorial1/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial1/tutorial.html/tutorial.html new file mode 100644 index 000000000..75cbcce57 --- /dev/null +++ b/docs/source/tutorials/tutorial1/tutorial.html/tutorial.html @@ -0,0 +1,8324 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+
+ +
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial10/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial10/tutorial.html/tutorial.html new file mode 100644 index 000000000..989b6f39a --- /dev/null +++ b/docs/source/tutorials/tutorial10/tutorial.html/tutorial.html @@ -0,0 +1,8071 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial11/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial11/tutorial.html/tutorial.html new file mode 100644 index 000000000..c535aa15c --- /dev/null +++ b/docs/source/tutorials/tutorial11/tutorial.html/tutorial.html @@ -0,0 +1,28770 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial12/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial12/tutorial.html/tutorial.html new file mode 100644 index 000000000..b34b5f397 --- /dev/null +++ b/docs/source/tutorials/tutorial12/tutorial.html/tutorial.html @@ -0,0 +1,7827 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+
+ +
+ +
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+
+ +
+
+ + diff --git a/docs/source/tutorials/tutorial13/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial13/tutorial.html/tutorial.html new file mode 100644 index 000000000..f44f1d671 --- /dev/null +++ b/docs/source/tutorials/tutorial13/tutorial.html/tutorial.html @@ -0,0 +1,8149 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial14/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial14/tutorial.html/tutorial.html new file mode 100644 index 000000000..a6b080680 --- /dev/null +++ b/docs/source/tutorials/tutorial14/tutorial.html/tutorial.html @@ -0,0 +1,8212 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + diff --git a/docs/source/tutorials/tutorial2/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial2/tutorial.html/tutorial.html new file mode 100644 index 000000000..ba0fac15c --- /dev/null +++ b/docs/source/tutorials/tutorial2/tutorial.html/tutorial.html @@ -0,0 +1,8429 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial3/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial3/tutorial.html/tutorial.html new file mode 100644 index 000000000..78ca8fd94 --- /dev/null +++ b/docs/source/tutorials/tutorial3/tutorial.html/tutorial.html @@ -0,0 +1,8256 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+
+ +
+ +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial4/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial4/tutorial.html/tutorial.html new file mode 100644 index 000000000..00c28b99f --- /dev/null +++ b/docs/source/tutorials/tutorial4/tutorial.html/tutorial.html @@ -0,0 +1,9046 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial5/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial5/tutorial.html/tutorial.html new file mode 100644 index 000000000..38cdb476b --- /dev/null +++ b/docs/source/tutorials/tutorial5/tutorial.html/tutorial.html @@ -0,0 +1,8039 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial6/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial6/tutorial.html/tutorial.html new file mode 100644 index 000000000..b72fb7e5b --- /dev/null +++ b/docs/source/tutorials/tutorial6/tutorial.html/tutorial.html @@ -0,0 +1,8222 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+
+ + diff --git a/docs/source/tutorials/tutorial7/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial7/tutorial.html/tutorial.html new file mode 100644 index 000000000..60fcfd984 --- /dev/null +++ b/docs/source/tutorials/tutorial7/tutorial.html/tutorial.html @@ -0,0 +1,8090 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial8/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial8/tutorial.html/tutorial.html new file mode 100644 index 000000000..fdf5a80a0 --- /dev/null +++ b/docs/source/tutorials/tutorial8/tutorial.html/tutorial.html @@ -0,0 +1,8227 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ +
+ +
+
+ +
+ +
+
+ +
+ + +
+
+ +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/docs/source/tutorials/tutorial9/tutorial.html/tutorial.html b/docs/source/tutorials/tutorial9/tutorial.html/tutorial.html new file mode 100644 index 000000000..4e1324891 --- /dev/null +++ b/docs/source/tutorials/tutorial9/tutorial.html/tutorial.html @@ -0,0 +1,8008 @@ + + + + + +tutorial + + + + + + + + + + + + +
+
+ +
+ +
+
+ +
+ +
+
+ +
+
+ +
+
+ +
+ +
+
+ +
+ + +
+ + +
+
+ +
+ + +
+
+ +
+ + +
+
+ +
+
+ + + diff --git a/tutorials/tutorial1/tutorial.ipynb b/tutorials/tutorial1/tutorial.ipynb index 4a0d205c2..b2b7fbd84 100644 --- a/tutorials/tutorial1/tutorial.ipynb +++ b/tutorials/tutorial1/tutorial.ipynb @@ -87,25 +87,27 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import warnings\n", "\n", "from pina.problem import SpatialProblem, TimeDependentProblem\n", "from pina.domain import CartesianDomain\n", "\n", - "warnings.filterwarnings('ignore')\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", "\n", "class TimeSpaceODE(SpatialProblem, TimeDependentProblem):\n", - " \n", - " output_variables = ['u']\n", - " spatial_domain = CartesianDomain({'x': [0, 1]})\n", - " temporal_domain = CartesianDomain({'t': [0, 1]})\n", + "\n", + " output_variables = [\"u\"]\n", + " spatial_domain = CartesianDomain({\"x\": [0, 1]})\n", + " temporal_domain = CartesianDomain({\"t\": [0, 1]})\n", "\n", " # other stuff ..." ] @@ -152,6 +154,7 @@ "from pina.domain import CartesianDomain\n", "from pina.equation import Equation, FixedValue\n", "\n", + "\n", "# defining the ode equation\n", "def ode_equation(input_, output_):\n", "\n", @@ -164,6 +167,7 @@ " # calculate the residual and return it\n", " return u_x - u\n", "\n", + "\n", "class SimpleODE(SpatialProblem):\n", "\n", " output_variables = [\"u\"]\n", diff --git a/tutorials/tutorial1/tutorial.py b/tutorials/tutorial1/tutorial.py new file mode 100644 index 000000000..b8da9c8cb --- /dev/null +++ b/tutorials/tutorial1/tutorial.py @@ -0,0 +1,340 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Physics Informed Neural Networks on PINA +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial1/tutorial.ipynb) +# + +# In this tutorial, we will demonstrate a typical use case of **PINA** on a toy problem, following the standard API procedure. +# +#

+# PINA API +#

+# +# Specifically, the tutorial aims to introduce the following topics: +# +# * Explaining how to build **PINA** Problems, +# * Showing how to generate data for `PINN` training +# +# These are the two main steps needed **before** starting the modelling optimization (choose model and solver, and train). We will show each step in detail, and at the end, we will solve a simple Ordinary Differential Equation (ODE) problem using the `PINN` solver. + +# ## Build a PINA problem + +# Problem definition in the **PINA** framework is done by building a python `class`, which inherits from one or more problem classes (`SpatialProblem`, `TimeDependentProblem`, `ParametricProblem`, ...) depending on the nature of the problem. Below is an example: +# ### Simple Ordinary Differential Equation +# Consider the following: +# +# $$ +# \begin{equation} +# \begin{cases} +# \frac{d}{dx}u(x) &= u(x) \quad x\in(0,1)\\ +# u(x=0) &= 1 \\ +# \end{cases} +# \end{equation} +# $$ +# +# with the analytical solution $u(x) = e^x$. In this case, our ODE depends only on the spatial variable $x\in(0,1)$ , meaning that our `Problem` class is going to be inherited from the `SpatialProblem` class: +# +# ```python +# from pina.problem import SpatialProblem +# from pina.domain import CartesianProblem +# +# class SimpleODE(SpatialProblem): +# +# output_variables = ['u'] +# spatial_domain = CartesianProblem({'x': [0, 1]}) +# +# # other stuff ... +# ``` +# +# Notice that we define `output_variables` as a list of symbols, indicating the output variables of our equation (in this case only $u$), this is done because in **PINA** the `torch.Tensor`s are labelled, allowing the user maximal flexibility for the manipulation of the tensor. The `spatial_domain` variable indicates where the sample points are going to be sampled in the domain, in this case $x\in[0,1]$. +# +# What if our equation is also time-dependent? In this case, our `class` will inherit from both `SpatialProblem` and `TimeDependentProblem`: +# + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import warnings + +from pina.problem import SpatialProblem, TimeDependentProblem +from pina.domain import CartesianDomain + +warnings.filterwarnings("ignore") + + +class TimeSpaceODE(SpatialProblem, TimeDependentProblem): + + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 1]}) + temporal_domain = CartesianDomain({"t": [0, 1]}) + + # other stuff ... + + +# where we have included the `temporal_domain` variable, indicating the time domain wanted for the solution. +# +# In summary, using **PINA**, we can initialize a problem with a class which inherits from different base classes: `SpatialProblem`, `TimeDependentProblem`, `ParametricProblem`, and so on depending on the type of problem we are considering. Here are some examples (more on the official documentation): +# * ``SpatialProblem`` $\rightarrow$ a differential equation with spatial variable(s) ``spatial_domain`` +# * ``TimeDependentProblem`` $\rightarrow$ a time-dependent differential equation with temporal variable(s) ``temporal_domain`` +# * ``ParametricProblem`` $\rightarrow$ a parametrized differential equation with parametric variable(s) ``parameter_domain`` +# * ``AbstractProblem`` $\rightarrow$ any **PINA** problem inherits from here + +# ### Write the problem class +# +# Once the `Problem` class is initialized, we need to represent the differential equation in **PINA**. In order to do this, we need to load the **PINA** operators from `pina.operator` module. Again, we'll consider Equation (1) and represent it in **PINA**: + +# In[ ]: + + +import torch +import matplotlib.pyplot as plt + +from pina.problem import SpatialProblem +from pina.operator import grad +from pina import Condition +from pina.domain import CartesianDomain +from pina.equation import Equation, FixedValue + + +# defining the ode equation +def ode_equation(input_, output_): + + # computing the derivative + u_x = grad(output_, input_, components=["u"], d=["x"]) + + # extracting the u input variable + u = output_.extract(["u"]) + + # calculate the residual and return it + return u_x - u + + +class SimpleODE(SpatialProblem): + + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 1]}) + + domains = { + "x0": CartesianDomain({"x": 0.0}), + "D": CartesianDomain({"x": [0, 1]}), + } + + # conditions to hold + conditions = { + "bound_cond": Condition(domain="x0", equation=FixedValue(1.0)), + "phys_cond": Condition(domain="D", equation=Equation(ode_equation)), + } + + # defining the true solution + def solution(self, pts): + return torch.exp(pts.extract(["x"])) + + +problem = SimpleODE() + + +# After we define the `Problem` class, we need to write different class methods, where each method is a function returning a residual. These functions are the ones minimized during PINN optimization, given the initial conditions. For example, in the domain $[0,1]$, the ODE equation (`ode_equation`) must be satisfied. We represent this by returning the difference between subtracting the variable `u` from its gradient (the residual), which we hope to minimize to 0. This is done for all conditions. Notice that we do not pass directly a `python` function, but an `Equation` object, which is initialized with the `python` function. This is done so that all the computations and internal checks are done inside **PINA**. +# +# Once we have defined the function, we need to tell the neural network where these methods are to be applied. To do so, we use the `Condition` class. In the `Condition` class, we pass the location points and the equation we want minimized on those points (other possibilities are allowed, see the documentation for reference). +# +# Finally, it's possible to define a `solution` function, which can be useful if we want to plot the results and see how the real solution compares to the expected (true) solution. Notice that the `solution` function is a method of the `PINN` class, but it is not mandatory for problem definition. +# + +# ## Generate data +# +# Data for training can come in form of direct numerical simulation results, or points in the domains. In case we perform unsupervised learning, we just need the collocation points for training, i.e. points where we want to evaluate the neural network. Sampling point in **PINA** is very easy, here we show three examples using the `.discretise_domain` method of the `AbstractProblem` class. + +# In[ ]: + + +# sampling 20 points in [0, 1] through discretization in all locations +problem.discretise_domain(n=20, mode="grid", domains="all") + +# sampling 20 points in (0, 1) through latin hypercube sampling in D, and 1 point in x0 +problem.discretise_domain(n=20, mode="latin", domains=["D"]) +problem.discretise_domain(n=1, mode="random", domains=["x0"]) + +# sampling 20 points in (0, 1) randomly +problem.discretise_domain(n=20, mode="random") + + +# We are going to use latin hypercube points for sampling. We need to sample in all the conditions domains. In our case we sample in `D` and `x0`. + +# In[ ]: + + +# sampling for training +problem.discretise_domain(1, "random", domains=["x0"]) +problem.discretise_domain(20, "lh", domains=["D"]) + + +# The points are saved in a python `dict`, and can be accessed by calling the attribute `input_pts` of the problem + +# In[ ]: + + +print("Input points:", problem.discretised_domains) +print("Input points labels:", problem.discretised_domains["D"].labels) + + +# To visualize the sampled points we can use `matplotlib.pyplot`: + +# In[ ]: + + +for location in problem.input_pts: + coords = ( + problem.input_pts[location].extract(problem.spatial_variables).flatten() + ) + plt.scatter(coords, torch.zeros_like(coords), s=10, label=location) +plt.legend() + + +# ## Perform a small training + +# Once we have defined the problem and generated the data we can start the modelling. Here we will choose a `FeedForward` neural network available in `pina.model`, and we will train using the `PINN` solver from `pina.solver`. We highlight that this training is fairly simple, for more advanced stuff consider the tutorials in the ***Physics Informed Neural Networks*** section of ***Tutorials***. For training we use the `Trainer` class from `pina.trainer`. Here we show a very short training and some method for plotting the results. Notice that by default all relevant metrics (e.g. MSE error during training) are going to be tracked using a `lightning` logger, by default `CSVLogger`. If you want to track the metric by yourself without a logger, use `pina.callback.MetricTracker`. + +# In[ ]: + + +from pina import Trainer +from pina.solver import PINN +from pina.model import FeedForward +from lightning.pytorch.loggers import TensorBoardLogger +from pina.optim import TorchOptimizer + + +# build the model +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) + +# create the PINN object +pinn = PINN(problem, model, TorchOptimizer(torch.optim.Adam, lr=0.005)) + +# create the trainer +trainer = Trainer( + solver=pinn, + max_epochs=1500, + logger=TensorBoardLogger("tutorial_logs"), + accelerator="cpu", + train_size=1.0, + test_size=0.0, + val_size=0.0, + enable_model_summary=False, +) # we train on CPU and avoid model summary at beginning of training (optional) + +# train +trainer.train() + + +# After the training we can inspect trainer logged metrics (by default **PINA** logs mean square error residual loss). The logged metrics can be accessed online using one of the `Lightning` loggers. The final loss can be accessed by `trainer.logged_metrics` + +# In[27]: + + +# inspecting final loss +trainer.logged_metrics + + +# By using `matplotlib` we can also do some qualitative plots of the solution. + +# In[ ]: + + +pts = pinn.problem.spatial_domain.sample(256, "grid", variables="x") +predicted_output = pinn.forward(pts).extract("u").tensor.detach() +true_output = pinn.problem.solution(pts).detach() +fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 8)) +ax.plot(pts.extract(["x"]), predicted_output, label="Neural Network solution") +ax.plot(pts.extract(["x"]), true_output, label="True solution") +plt.legend() + + +# The solution is overlapped with the actual one, and they are barely indistinguishable. We can also take a look at the loss using `TensorBoard`: + +# In[ ]: + + +print("\nTo load TensorBoard run load_ext tensorboard on your terminal") +print( + "To visualize the loss you can run tensorboard --logdir 'tutorial_logs' on your terminal\n" +) +# # uncomment for running tensorboard +# %load_ext tensorboard +# %tensorboard --logdir=tutorial_logs + + +# As we can see the loss has not reached a minimum, suggesting that we could train for longer! Alternatively, we can also take look at the loss using callbacks. Here we use `MetricTracker` from `pina.callback`: + +# In[ ]: + + +from pina.callback import MetricTracker + +# create the model +newmodel = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) + +# create the PINN object +newpinn = PINN( + problem, newmodel, optimizer=TorchOptimizer(torch.optim.Adam, lr=0.005) +) + +# create the trainer +newtrainer = Trainer( + solver=newpinn, + max_epochs=1500, + logger=True, # enable parameter logging + callbacks=[MetricTracker()], + accelerator="cpu", + train_size=1.0, + test_size=0.0, + val_size=0.0, + enable_model_summary=False, +) # we train on CPU and avoid model summary at beginning of training (optional) + +# train +newtrainer.train() + +# plot loss +trainer_metrics = newtrainer.callbacks[0].metrics +loss = trainer_metrics["train_loss"] +epochs = range(len(loss)) +plt.plot(epochs, loss.cpu()) +# plotting +plt.xlabel("epoch") +plt.ylabel("loss") +plt.yscale("log") + + +# ## What's next? +# +# Congratulations on completing the introductory tutorial of **PINA**! There are several directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the finaly accuracy +# +# 2. Train the network using other types of models (see `pina.model`) +# +# 3. GPU training and speed benchmarking +# +# 4. Many more... diff --git a/tutorials/tutorial10/tutorial.ipynb b/tutorials/tutorial10/tutorial.ipynb index 8f9baf797..748cb163a 100644 --- a/tutorials/tutorial10/tutorial.ipynb +++ b/tutorials/tutorial10/tutorial.ipynb @@ -25,16 +25,17 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", - " # get the data\n", - " !mkdir \"data\"\n", - " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS.mat\" -O \"data/Data_KS.mat\"\n", - " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS2.mat\" -O \"data/Data_KS2.mat\"\n", + " !pip install \"pina-mathlab\"\n", + " # get the data\n", + " !mkdir \"data\"\n", + " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS.mat\" -O \"data/Data_KS.mat\"\n", + " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS2.mat\" -O \"data/Data_KS2.mat\"\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", @@ -46,7 +47,7 @@ "from pina.solver import SupervisedSolver\n", "from pina.problem.zoo import SupervisedProblem\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -106,17 +107,24 @@ ], "source": [ "# load data\n", - "data=io.loadmat(\"data/Data_KS.mat\")\n", + "data = io.loadmat(\"data/Data_KS.mat\")\n", "\n", "# converting to label tensor\n", - "initial_cond_train = LabelTensor(torch.tensor(data['initial_cond_train'], dtype=torch.float), ['t','x','u0'])\n", - "initial_cond_test = LabelTensor(torch.tensor(data['initial_cond_test'], dtype=torch.float), ['t','x','u0'])\n", - "sol_train = LabelTensor(torch.tensor(data['sol_train'], dtype=torch.float), ['u'])\n", - "sol_test = LabelTensor(torch.tensor(data['sol_test'], dtype=torch.float), ['u'])\n", - "\n", - "print('Data Loaded')\n", - "print(f' shape initial condition: {initial_cond_train.shape}')\n", - "print(f' shape solution: {sol_train.shape}')" + "initial_cond_train = LabelTensor(\n", + " torch.tensor(data[\"initial_cond_train\"], dtype=torch.float),\n", + " [\"t\", \"x\", \"u0\"],\n", + ")\n", + "initial_cond_test = LabelTensor(\n", + " torch.tensor(data[\"initial_cond_test\"], dtype=torch.float), [\"t\", \"x\", \"u0\"]\n", + ")\n", + "sol_train = LabelTensor(\n", + " torch.tensor(data[\"sol_train\"], dtype=torch.float), [\"u\"]\n", + ")\n", + "sol_test = LabelTensor(torch.tensor(data[\"sol_test\"], dtype=torch.float), [\"u\"])\n", + "\n", + "print(\"Data Loaded\")\n", + "print(f\" shape initial condition: {initial_cond_train.shape}\")\n", + "print(f\" shape solution: {sol_train.shape}\")" ] }, { @@ -151,35 +159,58 @@ "# helper function\n", "def plot_trajectory(coords, real, no_sol=None):\n", " # find the x-t shapes\n", - " dim_x = len(torch.unique(coords.extract('x')))\n", - " dim_t = len(torch.unique(coords.extract('t')))\n", + " dim_x = len(torch.unique(coords.extract(\"x\")))\n", + " dim_t = len(torch.unique(coords.extract(\"t\")))\n", " # if we don't have the Neural Operator solution we simply plot the real one\n", " if no_sol is None:\n", " fig, axs = plt.subplots(1, 1, figsize=(15, 5), sharex=True, sharey=True)\n", - " c = axs.imshow(real.reshape(dim_t, dim_x).T.detach(),extent=[0, 50, 0, 64], cmap='PuOr_r', aspect='auto')\n", - " axs.set_title('Real solution')\n", + " c = axs.imshow(\n", + " real.reshape(dim_t, dim_x).T.detach(),\n", + " extent=[0, 50, 0, 64],\n", + " cmap=\"PuOr_r\",\n", + " aspect=\"auto\",\n", + " )\n", + " axs.set_title(\"Real solution\")\n", " fig.colorbar(c, ax=axs)\n", - " axs.set_xlabel('t')\n", - " axs.set_ylabel('x')\n", + " axs.set_xlabel(\"t\")\n", + " axs.set_ylabel(\"x\")\n", " # otherwise we plot the real one, the Neural Operator one, and their difference\n", " else:\n", " fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharex=True, sharey=True)\n", - " axs[0].imshow(real.reshape(dim_t, dim_x).T.detach(),extent=[0, 50, 0, 64], cmap='PuOr_r', aspect='auto')\n", - " axs[0].set_title('Real solution')\n", - " axs[1].imshow(no_sol.reshape(dim_t, dim_x).T.detach(),extent=[0, 50, 0, 64], cmap='PuOr_r', aspect='auto')\n", - " axs[1].set_title('NO solution')\n", - " c = axs[2].imshow((real - no_sol).abs().reshape(dim_t, dim_x).T.detach(),extent=[0, 50, 0, 64], cmap='PuOr_r', aspect='auto')\n", - " axs[2].set_title('Absolute difference')\n", + " axs[0].imshow(\n", + " real.reshape(dim_t, dim_x).T.detach(),\n", + " extent=[0, 50, 0, 64],\n", + " cmap=\"PuOr_r\",\n", + " aspect=\"auto\",\n", + " )\n", + " axs[0].set_title(\"Real solution\")\n", + " axs[1].imshow(\n", + " no_sol.reshape(dim_t, dim_x).T.detach(),\n", + " extent=[0, 50, 0, 64],\n", + " cmap=\"PuOr_r\",\n", + " aspect=\"auto\",\n", + " )\n", + " axs[1].set_title(\"NO solution\")\n", + " c = axs[2].imshow(\n", + " (real - no_sol).abs().reshape(dim_t, dim_x).T.detach(),\n", + " extent=[0, 50, 0, 64],\n", + " cmap=\"PuOr_r\",\n", + " aspect=\"auto\",\n", + " )\n", + " axs[2].set_title(\"Absolute difference\")\n", " fig.colorbar(c, ax=axs.ravel().tolist())\n", " for ax in axs:\n", - " ax.set_xlabel('t')\n", - " ax.set_ylabel('x')\n", + " ax.set_xlabel(\"t\")\n", + " ax.set_ylabel(\"x\")\n", " plt.show()\n", "\n", + "\n", "# a sample trajectory (we use the sample 5, feel free to change)\n", "sample_number = 20\n", - "plot_trajectory(coords=initial_cond_train[sample_number].extract(['x', 't']),\n", - " real=sol_train[sample_number].extract('u'))\n" + "plot_trajectory(\n", + " coords=initial_cond_train[sample_number].extract([\"x\", \"t\"]),\n", + " real=sol_train[sample_number].extract(\"u\"),\n", + ")" ] }, { @@ -300,7 +331,12 @@ ], "source": [ "# initialize problem\n", - "problem = SupervisedProblem(initial_cond_train, sol_train, input_variables=initial_cond_train.labels, output_variables=sol_train.labels)\n", + "problem = SupervisedProblem(\n", + " initial_cond_train,\n", + " sol_train,\n", + " input_variables=initial_cond_train.labels,\n", + " output_variables=sol_train.labels,\n", + ")\n", "# initialize solver\n", "solver = SupervisedSolver(problem=problem, model=model)\n", "# train, only CPU and avoid model summary at beginning of training (optional)\n", @@ -343,9 +379,11 @@ "source": [ "sample_number = 2\n", "no_sol = solver(initial_cond_test)\n", - "plot_trajectory(coords=initial_cond_test[sample_number].extract(['x', 't']),\n", - " real=sol_test[sample_number].extract('u'),\n", - " no_sol=no_sol[5])" + "plot_trajectory(\n", + " coords=initial_cond_test[sample_number].extract([\"x\", \"t\"]),\n", + " real=sol_test[sample_number].extract(\"u\"),\n", + " no_sol=no_sol[5],\n", + ")" ] }, { @@ -373,15 +411,19 @@ "source": [ "from pina.loss import PowerLoss\n", "\n", - "error_metric = PowerLoss(p=2) # we use the MSE loss\n", + "error_metric = PowerLoss(p=2) # we use the MSE loss\n", "\n", "with torch.no_grad():\n", " no_sol_train = solver(initial_cond_train)\n", - " err_train = error_metric(sol_train.extract('u'), no_sol_train).mean() # we average the error over trajectories\n", + " err_train = error_metric(\n", + " sol_train.extract(\"u\"), no_sol_train\n", + " ).mean() # we average the error over trajectories\n", " no_sol_test = solver(initial_cond_test)\n", - " err_test = error_metric(sol_test.extract('u'),no_sol_test).mean() # we average the error over trajectories\n", - " print(f'Training error: {float(err_train):.3f}')\n", - " print(f'Testing error: {float(err_test):.3f}')" + " err_test = error_metric(\n", + " sol_test.extract(\"u\"), no_sol_test\n", + " ).mean() # we average the error over trajectories\n", + " print(f\"Training error: {float(err_train):.3f}\")\n", + " print(f\"Testing error: {float(err_test):.3f}\")" ] }, { diff --git a/tutorials/tutorial10/tutorial.py b/tutorials/tutorial10/tutorial.py new file mode 100644 index 000000000..19f1ff7c4 --- /dev/null +++ b/tutorials/tutorial10/tutorial.py @@ -0,0 +1,318 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Averaging Neural Operator for solving Kuramoto Sivashinsky equation +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial10/tutorial.ipynb) +# +# In this tutorial we will build a Neural Operator using the +# `AveragingNeuralOperator` model and the `SupervisedSolver`. At the end of the +# tutorial you will be able to train a Neural Operator for learning +# the operator of time dependent PDEs. +# +# +# First of all, some useful imports. Note we use `scipy` for i/o operations. +# + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + # get the data + get_ipython().system('mkdir "data"') + get_ipython().system( + 'wget "https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS.mat" -O "data/Data_KS.mat"' + ) + get_ipython().system( + 'wget "https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial10/data/Data_KS2.mat" -O "data/Data_KS2.mat"' + ) + +import torch +import matplotlib.pyplot as plt +import warnings + +from scipy import io +from pina import Condition, Trainer, LabelTensor +from pina.model import AveragingNeuralOperator +from pina.solver import SupervisedSolver +from pina.problem.zoo import SupervisedProblem + +warnings.filterwarnings("ignore") + + +# ## Data Generation +# +# We will focus on solving a specific PDE, the **Kuramoto Sivashinsky** (KS) equation. +# The KS PDE is a fourth-order nonlinear PDE with the following form: +# +# $$ +# \frac{\partial u}{\partial t}(x,t) = -u(x,t)\frac{\partial u}{\partial x}(x,t)- \frac{\partial^{4}u}{\partial x^{4}}(x,t) - \frac{\partial^{2}u}{\partial x^{2}}(x,t). +# $$ +# +# In the above $x\in \Omega=[0, 64]$ represents a spatial location, $t\in\mathbb{T}=[0,50]$ the time and $u(x, t)$ is the value of the function $u:\Omega \times\mathbb{T}\in\mathbb{R}$. We indicate with $\mathbb{U}$ a suitable space for $u$, i.e. we have that the solution $u\in\mathbb{U}$. +# +# +# We impose Dirichlet boundary conditions on the derivative of $u$ on the border of the domain $\partial \Omega$ +# $$ +# \frac{\partial u}{\partial x}(x,t)=0 \quad \forall (x,t)\in \partial \Omega\times\mathbb{T}. +# $$ +# +# Initial conditions are sampled from a distribution over truncated Fourier series with random coefficients +# $\{A_k, \ell_k, \phi_k\}_k$ as +# $$ +# u(x,0) = \sum_{k=1}^N A_k \sin(2 \pi \ell_k x / L + \phi_k) \ , +# $$ +# +# where $A_k \in [-0.4, -0.3]$, $\ell_k = 2$, $\phi_k = 2\pi \quad \forall k=1,\dots,N$. +# +# +# We have already generated some data for differenti initial conditions, and our objective will +# be to build a Neural Operator that, given $u(x, t)$ will output $u(x, t+\delta)$, where +# $\delta$ is a fixed time step. We will come back on the Neural Operator architecture, for now +# we first need to import the data. +# +# **Note:** +# *The numerical integration is obtained by using pseudospectral method for spatial derivative discratization and +# implicit Runge Kutta 5 for temporal dynamics.* +# + +# In[2]: + + +# load data +data = io.loadmat("data/Data_KS.mat") + +# converting to label tensor +initial_cond_train = LabelTensor( + torch.tensor(data["initial_cond_train"], dtype=torch.float), + ["t", "x", "u0"], +) +initial_cond_test = LabelTensor( + torch.tensor(data["initial_cond_test"], dtype=torch.float), ["t", "x", "u0"] +) +sol_train = LabelTensor( + torch.tensor(data["sol_train"], dtype=torch.float), ["u"] +) +sol_test = LabelTensor(torch.tensor(data["sol_test"], dtype=torch.float), ["u"]) + +print("Data Loaded") +print(f" shape initial condition: {initial_cond_train.shape}") +print(f" shape solution: {sol_train.shape}") + + +# The data are saved in the form `B \times N \times D`, where `B` is the batch_size +# (basically how many initial conditions we sample), `N` the number of points in the mesh +# (which is the product of the discretization in `x` timese the one in `t`), and +# `D` the dimension of the problem (in this case we have three variables `[u, t, x]`). +# +# We are now going to plot some trajectories! + +# In[3]: + + +# helper function +def plot_trajectory(coords, real, no_sol=None): + # find the x-t shapes + dim_x = len(torch.unique(coords.extract("x"))) + dim_t = len(torch.unique(coords.extract("t"))) + # if we don't have the Neural Operator solution we simply plot the real one + if no_sol is None: + fig, axs = plt.subplots(1, 1, figsize=(15, 5), sharex=True, sharey=True) + c = axs.imshow( + real.reshape(dim_t, dim_x).T.detach(), + extent=[0, 50, 0, 64], + cmap="PuOr_r", + aspect="auto", + ) + axs.set_title("Real solution") + fig.colorbar(c, ax=axs) + axs.set_xlabel("t") + axs.set_ylabel("x") + # otherwise we plot the real one, the Neural Operator one, and their difference + else: + fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharex=True, sharey=True) + axs[0].imshow( + real.reshape(dim_t, dim_x).T.detach(), + extent=[0, 50, 0, 64], + cmap="PuOr_r", + aspect="auto", + ) + axs[0].set_title("Real solution") + axs[1].imshow( + no_sol.reshape(dim_t, dim_x).T.detach(), + extent=[0, 50, 0, 64], + cmap="PuOr_r", + aspect="auto", + ) + axs[1].set_title("NO solution") + c = axs[2].imshow( + (real - no_sol).abs().reshape(dim_t, dim_x).T.detach(), + extent=[0, 50, 0, 64], + cmap="PuOr_r", + aspect="auto", + ) + axs[2].set_title("Absolute difference") + fig.colorbar(c, ax=axs.ravel().tolist()) + for ax in axs: + ax.set_xlabel("t") + ax.set_ylabel("x") + plt.show() + + +# a sample trajectory (we use the sample 5, feel free to change) +sample_number = 20 +plot_trajectory( + coords=initial_cond_train[sample_number].extract(["x", "t"]), + real=sol_train[sample_number].extract("u"), +) + + +# As we can see, as the time progresses the solution becomes chaotic, which makes +# it really hard to learn! We will now focus on building a Neural Operator using the +# `SupervisedSolver` class to tackle the problem. +# +# ## Averaging Neural Operator +# +# We will build a neural operator $\texttt{NO}$ which takes the solution at time $t=0$ for any $x\in\Omega$, +# the time $(t)$ at which we want to compute the solution, and gives back the solution to the KS equation $u(x, t)$, mathematically: +# $$ +# \texttt{NO}_\theta : \mathbb{U} \rightarrow \mathbb{U}, +# $$ +# such that +# $$ +# \texttt{NO}_\theta[u(t=0)](x, t) \rightarrow u(x, t). +# $$ +# +# There are many ways on approximating the following operator, e.g. by 2D [FNO](https://mathlab.github.io/PINA/_rst/models/fno.html) (for regular meshes), +# a [DeepOnet](https://mathlab.github.io/PINA/_rst/models/deeponet.html), [Continuous Convolutional Neural Operator](https://mathlab.github.io/PINA/_rst/layers/convolution.html), +# [MIONet](https://mathlab.github.io/PINA/_rst/models/mionet.html). +# In this tutorial we will use the *Averaging Neural Operator* presented in [*The Nonlocal Neural Operator: Universal Approximation*](https://arxiv.org/abs/2304.13221) +# which is a [Kernel Neural Operator](https://mathlab.github.io/PINA/_rst/models/base_no.html) with integral kernel: +# +# $$ +# K(v) = \sigma\left(Wv(x) + b + \frac{1}{|\Omega|}\int_\Omega v(y)dy\right) +# $$ +# +# where: +# +# * $v(x)\in\mathbb{R}^{\rm{emb}}$ is the update for a function $v$ with $\mathbb{R}^{\rm{emb}}$ the embedding (hidden) size +# * $\sigma$ is a non-linear activation +# * $W\in\mathbb{R}^{\rm{emb}\times\rm{emb}}$ is a tunable matrix. +# * $b\in\mathbb{R}^{\rm{emb}}$ is a tunable bias. +# +# If PINA many Kernel Neural Operators are already implemented, and the modular componets of the [Kernel Neural Operator](https://mathlab.github.io/PINA/_rst/models/base_no.html) class permits to create new ones by composing base kernel layers. +# +# **Note:*** We will use the already built class* `AveragingNeuralOperator`, *as constructive excercise try to use the* [KernelNeuralOperator](https://mathlab.github.io/PINA/_rst/models/base_no.html) *class for building a kernel neural operator from scratch. You might employ the different layers that we have in pina, e.g.* [FeedForward](https://mathlab.github.io/PINA/_rst/models/fnn.html), *and* [AveragingNeuralOperator](https://mathlab.github.io/PINA/_rst/layers/avno_layer.html) *layers*. + +# In[4]: + + +class SIREN(torch.nn.Module): + def forward(self, x): + return torch.sin(x) + + +embedding_dimesion = 40 # hyperparameter embedding dimension +input_dimension = 3 # ['u', 'x', 't'] +number_of_coordinates = 2 # ['x', 't'] +lifting_net = torch.nn.Linear(input_dimension, embedding_dimesion) +projecting_net = torch.nn.Linear(embedding_dimesion + number_of_coordinates, 1) +model = AveragingNeuralOperator( + lifting_net=lifting_net, + projecting_net=projecting_net, + coordinates_indices=["x", "t"], + field_indices=["u0"], + n_layers=4, + func=SIREN, +) + + +# Super easy! Notice that we use the `SIREN` activation function, more on [Implicit Neural Representations with Periodic Activation Functions](https://arxiv.org/abs/2006.09661). +# +# ## Solving the KS problem +# +# We will now focus on solving the KS equation using the `SupervisedSolver` class +# and the `AveragingNeuralOperator` model. As done in the [FNO tutorial](https://github.com/mathLab/PINA/blob/master/tutorials/tutorial5/tutorial.ipynb) we now create the Neural Operator problem class with `SupervisedProblem`. + +# In[5]: + + +# initialize problem +problem = SupervisedProblem( + initial_cond_train, + sol_train, + input_variables=initial_cond_train.labels, + output_variables=sol_train.labels, +) +# initialize solver +solver = SupervisedSolver(problem=problem, model=model) +# train, only CPU and avoid model summary at beginning of training (optional) +trainer = Trainer( + solver=solver, + max_epochs=40, + accelerator="cpu", + enable_model_summary=False, + batch_size=5, # we train on CPU and avoid model summary at beginning of training (optional) + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# We can now see some plots for the solutions + +# In[6]: + + +sample_number = 2 +no_sol = solver(initial_cond_test) +plot_trajectory( + coords=initial_cond_test[sample_number].extract(["x", "t"]), + real=sol_test[sample_number].extract("u"), + no_sol=no_sol[5], +) + + +# As we can see we can obtain nice result considering the small training time and the difficulty of the problem! +# Let's take a look at the training and testing error: + +# In[7]: + + +from pina.loss import PowerLoss + +error_metric = PowerLoss(p=2) # we use the MSE loss + +with torch.no_grad(): + no_sol_train = solver(initial_cond_train) + err_train = error_metric( + sol_train.extract("u"), no_sol_train + ).mean() # we average the error over trajectories + no_sol_test = solver(initial_cond_test) + err_test = error_metric( + sol_test.extract("u"), no_sol_test + ).mean() # we average the error over trajectories + print(f"Training error: {float(err_train):.3f}") + print(f"Testing error: {float(err_test):.3f}") + + +# As we can see the error is pretty small, which agrees with what we can see from the previous plots. + +# ## What's next? +# +# Now you know how to solve a time dependent neural operator problem in **PINA**! There are multiple directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the final accuracy +# +# 2. We left a more challenging dataset [Data_KS2.mat](dat/Data_KS2.mat) where $A_k \in [-0.5, 0.5]$, $\ell_k \in [1, 2, 3]$, $\phi_k \in [0, 2\pi]$ for longer training +# +# 3. Compare the performance between the different neural operators (you can even try to implement your favourite one!) diff --git a/tutorials/tutorial11/tutorial.ipynb b/tutorials/tutorial11/tutorial.ipynb index 3506ca7ec..2a23eaecf 100644 --- a/tutorials/tutorial11/tutorial.ipynb +++ b/tutorials/tutorial11/tutorial.ipynb @@ -24,12 +24,13 @@ "outputs": [], "source": [ "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import torch\n", "import warnings\n", @@ -42,7 +43,7 @@ "from pina.domain import CartesianDomain\n", "from pina.equation import Equation, FixedValue\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -70,6 +71,7 @@ " # calculate the residual and return it\n", " return u_x - u\n", "\n", + "\n", "class SimpleODE(SpatialProblem):\n", "\n", " output_variables = [\"u\"]\n", @@ -101,7 +103,7 @@ " layers=[10, 10],\n", " func=torch.nn.Tanh,\n", " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", + " input_dimensions=len(problem.input_variables),\n", ")\n", "\n", "# create the PINN object\n", @@ -437,20 +439,22 @@ ], "source": [ "model = FeedForward(\n", - " layers=[10, 10],\n", - " func=torch.nn.Tanh,\n", - " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", - " )\n", + " layers=[10, 10],\n", + " func=torch.nn.Tanh,\n", + " output_dimensions=len(problem.output_variables),\n", + " input_dimensions=len(problem.input_variables),\n", + ")\n", "pinn = PINN(problem, model)\n", - "trainer = Trainer(solver=pinn,\n", - " accelerator='cpu',\n", - " logger=True,\n", - " callbacks=[NaiveMetricTracker()], # adding a callbacks\n", - " enable_model_summary=False,\n", - " train_size=1.0,\n", - " val_size=0.0,\n", - " test_size=0.0)\n", + "trainer = Trainer(\n", + " solver=pinn,\n", + " accelerator=\"cpu\",\n", + " logger=True,\n", + " callbacks=[NaiveMetricTracker()], # adding a callbacks\n", + " enable_model_summary=False,\n", + " train_size=1.0,\n", + " val_size=0.0,\n", + " test_size=0.0,\n", + ")\n", "trainer.train()" ] }, @@ -486,7 +490,7 @@ } ], "source": [ - "trainer.callbacks[0].saved_metrics[:3] # only the first three epochs" + "trainer.callbacks[0].saved_metrics[:3] # only the first three epochs" ] }, { @@ -615,19 +619,21 @@ "seed_everything(42, workers=True)\n", "\n", "model = FeedForward(\n", - " layers=[10, 10],\n", - " func=torch.nn.Tanh,\n", - " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", - " )\n", + " layers=[10, 10],\n", + " func=torch.nn.Tanh,\n", + " output_dimensions=len(problem.output_variables),\n", + " input_dimensions=len(problem.input_variables),\n", + ")\n", "\n", "pinn = PINN(problem, model)\n", - "trainer = Trainer(solver=pinn,\n", - " accelerator='cpu',\n", - " deterministic=True, # setting deterministic=True ensure reproducibility when a seed is imposed\n", - " max_epochs = 2000,\n", - " enable_model_summary=False,\n", - " callbacks=[Timer()]) # adding a callbacks\n", + "trainer = Trainer(\n", + " solver=pinn,\n", + " accelerator=\"cpu\",\n", + " deterministic=True, # setting deterministic=True ensure reproducibility when a seed is imposed\n", + " max_epochs=2000,\n", + " enable_model_summary=False,\n", + " callbacks=[Timer()],\n", + ") # adding a callbacks\n", "trainer.train()\n", "print(f'Total training time {trainer.callbacks[0].time_elapsed(\"train\"):.5f} s')" ] @@ -698,19 +704,20 @@ "seed_everything(42, workers=True)\n", "\n", "model = FeedForward(\n", - " layers=[10, 10],\n", - " func=torch.nn.Tanh,\n", - " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", - " )\n", + " layers=[10, 10],\n", + " func=torch.nn.Tanh,\n", + " output_dimensions=len(problem.output_variables),\n", + " input_dimensions=len(problem.input_variables),\n", + ")\n", "pinn = PINN(problem, model)\n", - "trainer = Trainer(solver=pinn,\n", - " accelerator='cpu',\n", - " deterministic=True,\n", - " max_epochs = 2000,\n", - " enable_model_summary=False,\n", - " callbacks=[Timer(),\n", - " StochasticWeightAveraging(swa_lrs=0.005)]) # adding StochasticWeightAveraging callbacks\n", + "trainer = Trainer(\n", + " solver=pinn,\n", + " accelerator=\"cpu\",\n", + " deterministic=True,\n", + " max_epochs=2000,\n", + " enable_model_summary=False,\n", + " callbacks=[Timer(), StochasticWeightAveraging(swa_lrs=0.005)],\n", + ") # adding StochasticWeightAveraging callbacks\n", "trainer.train()\n", "print(f'Total training time {trainer.callbacks[0].time_elapsed(\"train\"):.5f} s')" ] @@ -783,19 +790,20 @@ "seed_everything(42, workers=True)\n", "\n", "model = FeedForward(\n", - " layers=[10, 10],\n", - " func=torch.nn.Tanh,\n", - " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", - " )\n", + " layers=[10, 10],\n", + " func=torch.nn.Tanh,\n", + " output_dimensions=len(problem.output_variables),\n", + " input_dimensions=len(problem.input_variables),\n", + ")\n", "pinn = PINN(problem, model)\n", - "trainer = Trainer(solver=pinn,\n", - " accelerator='cpu',\n", - " max_epochs = 2000,\n", - " enable_model_summary=False,\n", - " gradient_clip_val=0.1, # clipping the gradient\n", - " callbacks=[Timer(),\n", - " StochasticWeightAveraging(swa_lrs=0.005)])\n", + "trainer = Trainer(\n", + " solver=pinn,\n", + " accelerator=\"cpu\",\n", + " max_epochs=2000,\n", + " enable_model_summary=False,\n", + " gradient_clip_val=0.1, # clipping the gradient\n", + " callbacks=[Timer(), StochasticWeightAveraging(swa_lrs=0.005)],\n", + ")\n", "trainer.train()\n", "print(f'Total training time {trainer.callbacks[0].time_elapsed(\"train\"):.5f} s')" ] diff --git a/tutorials/tutorial11/tutorial.py b/tutorials/tutorial11/tutorial.py new file mode 100644 index 000000000..0fd10295e --- /dev/null +++ b/tutorials/tutorial11/tutorial.py @@ -0,0 +1,387 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: PINA and PyTorch Lightning, training tips and visualizations +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial11/tutorial.ipynb) +# +# In this tutorial, we will delve deeper into the functionality of the `Trainer` class, which serves as the cornerstone for training **PINA** [Solvers](https://mathlab.github.io/PINA/_rst/_code.html#solvers). +# +# The `Trainer` class offers a plethora of features aimed at improving model accuracy, reducing training time and memory usage, facilitating logging visualization, and more thanks to the amazing job done by the PyTorch Lightning team! +# +# Our leading example will revolve around solving the `SimpleODE` problem, as outlined in the [*Introduction to PINA for Physics Informed Neural Networks training*](https://github.com/mathLab/PINA/blob/master/tutorials/tutorial1/tutorial.ipynb). If you haven't already explored it, we highly recommend doing so before diving into this tutorial. +# +# Let's start by importing useful modules, define the `SimpleODE` problem and the `PINN` solver. + +# In[1]: + + +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import warnings + +from pina import Condition, Trainer +from pina.solver import PINN +from pina.model import FeedForward +from pina.problem import SpatialProblem +from pina.operator import grad +from pina.domain import CartesianDomain +from pina.equation import Equation, FixedValue + +warnings.filterwarnings("ignore") + + +# Define problem and solver. + +# In[2]: + + +# defining the ode equation +def ode_equation(input_, output_): + + # computing the derivative + u_x = grad(output_, input_, components=["u"], d=["x"]) + + # extracting the u input variable + u = output_.extract(["u"]) + + # calculate the residual and return it + return u_x - u + + +class SimpleODE(SpatialProblem): + + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 1]}) + + domains = { + "x0": CartesianDomain({"x": 0.0}), + "D": CartesianDomain({"x": [0, 1]}), + } + + # conditions to hold + conditions = { + "bound_cond": Condition(domain="x0", equation=FixedValue(1.0)), + "phys_cond": Condition(domain="D", equation=Equation(ode_equation)), + } + + # defining the true solution + def solution(self, pts): + return torch.exp(pts.extract(["x"])) + + +# sampling for training +problem = SimpleODE() +problem.discretise_domain(1, "random", domains=["x0"]) +problem.discretise_domain(20, "lh", domains=["D"]) + +# build the model +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) + +# create the PINN object +pinn = PINN(problem, model) + + +# Till now we just followed the extact step of the previous tutorials. The `Trainer` object +# can be initialized by simiply passing the `PINN` solver + +# In[3]: + + +trainer = Trainer(solver=pinn) + + +# ## Trainer Accelerator +# +# When creating the trainer, **by defualt** the `Trainer` will choose the most performing `accelerator` for training which is available in your system, ranked as follow: +# 1. [TPU](https://cloud.google.com/tpu/docs/intro-to-tpu) +# 2. [IPU](https://www.graphcore.ai/products/ipu) +# 3. [HPU](https://habana.ai/) +# 4. [GPU](https://www.intel.com/content/www/us/en/products/docs/processors/what-is-a-gpu.html#:~:text=What%20does%20GPU%20stand%20for,video%20editing%2C%20and%20gaming%20applications) or [MPS](https://developer.apple.com/metal/pytorch/) +# 5. CPU + +# For setting manually the `accelerator` run: +# +# * `accelerator = {'gpu', 'cpu', 'hpu', 'mps', 'cpu', 'ipu'}` sets the accelerator to a specific one + +# In[4]: + + +trainer = Trainer(solver=pinn, accelerator="cpu") + + +# as you can see, even if in the used system `GPU` is available, it is not used since we set `accelerator='cpu'`. + +# ## Trainer Logging +# +# In **PINA** you can log metrics in different ways. The simplest approach is to use the `MetricTraker` class from `pina.callbacks` as seen in the [*Introduction to PINA for Physics Informed Neural Networks training*](https://github.com/mathLab/PINA/blob/master/tutorials/tutorial1/tutorial.ipynb) tutorial. +# +# However, expecially when we need to train multiple times to get an average of the loss across multiple runs, `pytorch_lightning.loggers` might be useful. Here we will use `TensorBoardLogger` (more on [logging](https://lightning.ai/docs/pytorch/stable/extensions/logging.html) here), but you can choose the one you prefer (or make your own one). +# +# We will now import `TensorBoardLogger`, do three runs of training and then visualize the results. Notice we set `enable_model_summary=False` to avoid model summary specifications (e.g. number of parameters), set it to true if needed. +# + +# In[5]: + + +from lightning.pytorch.loggers import TensorBoardLogger + +# three run of training, by default it trains for 1000 epochs +# we reinitialize the model each time otherwise the same parameters will be optimized +for _ in range(3): + model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), + ) + pinn = PINN(problem, model) + trainer = Trainer( + solver=pinn, + accelerator="cpu", + logger=TensorBoardLogger(save_dir="training_log"), + enable_model_summary=False, + train_size=1.0, + val_size=0.0, + test_size=0.0, + ) + trainer.train() + + +# We can now visualize the logs by simply running `tensorboard --logdir=training_log/` on terminal, you should obtain a webpage as the one shown below: + +#

+# \"Logging +#

+ +# as you can see, by default, **PINA** logs the losses which are shown in the progress bar, as well as the number of epochs. You can always insert more loggings by either defining a **callback** ([more on callbacks](https://lightning.ai/docs/pytorch/stable/extensions/callbacks.html)), or inheriting the solver and modify the programs with different **hooks** ([more on hooks](https://lightning.ai/docs/pytorch/stable/common/lightning_module.html#hooks)). + +# ## Trainer Callbacks + +# Whenever we need to access certain steps of the training for logging, do static modifications (i.e. not changing the `Solver`) or updating `Problem` hyperparameters (static variables), we can use `Callabacks`. Notice that `Callbacks` allow you to add arbitrary self-contained programs to your training. At specific points during the flow of execution (hooks), the Callback interface allows you to design programs that encapsulate a full set of functionality. It de-couples functionality that does not need to be in **PINA** `Solver`s. +# Lightning has a callback system to execute them when needed. Callbacks should capture NON-ESSENTIAL logic that is NOT required for your lightning module to run. +# +# The following are best practices when using/designing callbacks. +# +# * Callbacks should be isolated in their functionality. +# * Your callback should not rely on the behavior of other callbacks in order to work properly. +# * Do not manually call methods from the callback. +# * Directly calling methods (eg. on_validation_end) is strongly discouraged. +# * Whenever possible, your callbacks should not depend on the order in which they are executed. +# +# We will try now to implement a naive version of `MetricTraker` to show how callbacks work. Notice that this is a very easy application of callbacks, fortunately in **PINA** we already provide more advanced callbacks in `pina.callbacks`. +# +# + +# In[6]: + + +from lightning.pytorch.callbacks import Callback +from lightning.pytorch.callbacks import EarlyStopping +import torch + + +# define a simple callback +class NaiveMetricTracker(Callback): + def __init__(self): + self.saved_metrics = [] + + def on_train_epoch_end( + self, trainer, __ + ): # function called at the end of each epoch + self.saved_metrics.append( + {key: value for key, value in trainer.logged_metrics.items()} + ) + + +# Let's see the results when applyed to the `SimpleODE` problem. You can define callbacks when initializing the `Trainer` by the `callbacks` argument, which expects a list of callbacks. + +# In[7]: + + +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) +pinn = PINN(problem, model) +trainer = Trainer( + solver=pinn, + accelerator="cpu", + logger=True, + callbacks=[NaiveMetricTracker()], # adding a callbacks + enable_model_summary=False, + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# We can easily access the data by calling `trainer.callbacks[0].saved_metrics` (notice the zero representing the first callback in the list given at initialization). + +# In[8]: + + +trainer.callbacks[0].saved_metrics[:3] # only the first three epochs + + +# PyTorch Lightning also has some built in `Callbacks` which can be used in **PINA**, [here an extensive list](https://lightning.ai/docs/pytorch/stable/extensions/callbacks.html#built-in-callbacks). +# +# We can for example try the `EarlyStopping` routine, which automatically stops the training when a specific metric converged (here the `train_loss`). In order to let the training keep going forever set `max_epochs=-1`. + +# In[9]: + + +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) +pinn = PINN(problem, model) +trainer = Trainer( + solver=pinn, + accelerator="cpu", + max_epochs=-1, + enable_model_summary=False, + val_size=0.2, + train_size=0.8, + test_size=0.0, + callbacks=[EarlyStopping("val_loss")], +) # adding a callbacks +trainer.train() + + +# As we can see the model automatically stop when the logging metric stopped improving! + +# ## Trainer Tips to Boost Accuracy, Save Memory and Speed Up Training +# +# Untill now we have seen how to choose the right `accelerator`, how to log and visualize the results, and how to interface with the program in order to add specific parts of code at specific points by `callbacks`. +# Now, we well focus on how boost your training by saving memory and speeding it up, while mantaining the same or even better degree of accuracy! +# +# +# There are several built in methods developed in PyTorch Lightning which can be applied straight forward in **PINA**, here we report some: +# +# * [Stochastic Weight Averaging](https://pytorch.org/blog/pytorch-1.6-now-includes-stochastic-weight-averaging/) to boost accuracy +# * [Gradient Clippling](https://deepgram.com/ai-glossary/gradient-clipping) to reduce computational time (and improve accuracy) +# * [Gradient Accumulation](https://lightning.ai/docs/pytorch/stable/common/optimization.html#id3) to save memory consumption +# * [Mixed Precision Training](https://lightning.ai/docs/pytorch/stable/common/optimization.html#id3) to save memory consumption +# +# We will just demonstrate how to use the first two, and see the results compared to a standard training. +# We use the [`Timer`](https://lightning.ai/docs/pytorch/stable/api/lightning.pytorch.callbacks.Timer.html#lightning.pytorch.callbacks.Timer) callback from `pytorch_lightning.callbacks` to take the times. Let's start by training a simple model without any optimization (train for 2000 epochs). + +# In[10]: + + +from lightning.pytorch.callbacks import Timer +from lightning.pytorch import seed_everything + +# setting the seed for reproducibility +seed_everything(42, workers=True) + +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) + +pinn = PINN(problem, model) +trainer = Trainer( + solver=pinn, + accelerator="cpu", + deterministic=True, # setting deterministic=True ensure reproducibility when a seed is imposed + max_epochs=2000, + enable_model_summary=False, + callbacks=[Timer()], +) # adding a callbacks +trainer.train() +print(f'Total training time {trainer.callbacks[0].time_elapsed("train"):.5f} s') + + +# Now we do the same but with StochasticWeightAveraging + +# In[11]: + + +from lightning.pytorch.callbacks import StochasticWeightAveraging + +# setting the seed for reproducibility +seed_everything(42, workers=True) + +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) +pinn = PINN(problem, model) +trainer = Trainer( + solver=pinn, + accelerator="cpu", + deterministic=True, + max_epochs=2000, + enable_model_summary=False, + callbacks=[Timer(), StochasticWeightAveraging(swa_lrs=0.005)], +) # adding StochasticWeightAveraging callbacks +trainer.train() +print(f'Total training time {trainer.callbacks[0].time_elapsed("train"):.5f} s') + + +# As you can see, the training time does not change at all! Notice that around epoch `1600` +# the scheduler is switched from the defalut one `ConstantLR` to the Stochastic Weight Average Learning Rate (`SWALR`). +# This is because by default `StochasticWeightAveraging` will be activated after `int(swa_epoch_start * max_epochs)` with `swa_epoch_start=0.7` by default. Finally, the final `mean_loss` is lower when `StochasticWeightAveraging` is used. +# +# We will now now do the same but clippling the gradient to be relatively small. + +# In[12]: + + +# setting the seed for reproducibility +seed_everything(42, workers=True) + +model = FeedForward( + layers=[10, 10], + func=torch.nn.Tanh, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) +pinn = PINN(problem, model) +trainer = Trainer( + solver=pinn, + accelerator="cpu", + max_epochs=2000, + enable_model_summary=False, + gradient_clip_val=0.1, # clipping the gradient + callbacks=[Timer(), StochasticWeightAveraging(swa_lrs=0.005)], +) +trainer.train() +print(f'Total training time {trainer.callbacks[0].time_elapsed("train"):.5f} s') + + +# As we can see we by applying gradient clipping we were able to even obtain lower error! +# +# ## What's next? +# +# Now you know how to use efficiently the `Trainer` class **PINA**! There are multiple directions you can go now: +# +# 1. Explore training times on different devices (e.g.) `TPU` +# +# 2. Try to reduce memory cost by mixed precision training and gradient accumulation (especially useful when training Neural Operators) +# +# 3. Benchmark `Trainer` speed for different precisions. diff --git a/tutorials/tutorial12/tutorial.ipynb b/tutorials/tutorial12/tutorial.ipynb index 77538395a..0223da5ae 100644 --- a/tutorials/tutorial12/tutorial.ipynb +++ b/tutorials/tutorial12/tutorial.ipynb @@ -53,16 +53,17 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import torch\n", "\n", - "#useful imports\n", + "# useful imports\n", "from pina import Condition\n", "from pina.problem import SpatialProblem, TimeDependentProblem\n", "from pina.equation import Equation, FixedValue\n", @@ -166,22 +167,23 @@ "outputs": [], "source": [ "class Burgers1DEquation(Equation):\n", - " \n", - " def __init__(self, nu = 0.):\n", + "\n", + " def __init__(self, nu=0.0):\n", " \"\"\"\n", " Burgers1D class. This class can be\n", " used to enforce the solution u to solve the viscous Burgers 1D Equation.\n", - " \n", + "\n", " :param torch.float32 nu: the viscosity coefficient. Default value is set to 0.\n", " \"\"\"\n", - " self.nu = nu \n", - " \n", + " self.nu = nu\n", + "\n", " def equation(input_, output_):\n", - " return grad(output_, input_, d='t') +\\\n", - " output_*grad(output_, input_, d='x') -\\\n", - " self.nu*laplacian(output_, input_, d='x')\n", + " return (\n", + " grad(output_, input_, d=\"t\")\n", + " + output_ * grad(output_, input_, d=\"x\")\n", + " - self.nu * laplacian(output_, input_, d=\"x\")\n", + " )\n", "\n", - " \n", " super().__init__(equation)" ] }, diff --git a/tutorials/tutorial12/tutorial.py b/tutorials/tutorial12/tutorial.py new file mode 100644 index 000000000..300744081 --- /dev/null +++ b/tutorials/tutorial12/tutorial.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: The `Equation` Class +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial12/tutorial.ipynb) + +# In this tutorial, we will show how to use the `Equation` Class in PINA. Specifically, we will see how use the Class and its inherited classes to enforce residuals minimization in PINNs. + +# # Example: The Burgers 1D equation + +# We will start implementing the viscous Burgers 1D problem Class, described as follows: +# +# +# $$ +# \begin{equation} +# \begin{cases} +# \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} &= \nu \frac{\partial^2 u}{ \partial x^2}, \quad x\in(0,1), \quad t>0\\ +# u(x,0) &= -\sin (\pi x)\\ +# u(x,t) &= 0 \quad x = \pm 1\\ +# \end{cases} +# \end{equation} +# $$ +# +# where we set $ \nu = \frac{0.01}{\pi}$. +# +# In the class that models this problem we will see in action the `Equation` class and one of its inherited classes, the `FixedValue` class. + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch + +# useful imports +from pina import Condition +from pina.problem import SpatialProblem, TimeDependentProblem +from pina.equation import Equation, FixedValue +from pina.domain import CartesianDomain +from pina.operator import grad, laplacian + + +# In[2]: + + +# define the burger equation +def burger_equation(input_, output_): + du = grad(output_, input_) + ddu = grad(du, input_, components=["dudx"]) + return ( + du.extract(["dudt"]) + + output_.extract(["u"]) * du.extract(["dudx"]) + - (0.01 / torch.pi) * ddu.extract(["ddudxdx"]) + ) + + +# define initial condition +def initial_condition(input_, output_): + u_expected = -torch.sin(torch.pi * input_.extract(["x"])) + return output_.extract(["u"]) - u_expected + + +class Burgers1D(TimeDependentProblem, SpatialProblem): + + # assign output/ spatial and temporal variables + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [-1, 1]}) + temporal_domain = CartesianDomain({"t": [0, 1]}) + + domains = { + "bound_cond1": CartesianDomain({"x": -1, "t": [0, 1]}), + "bound_cond2": CartesianDomain({"x": 1, "t": [0, 1]}), + "time_cond": CartesianDomain({"x": [-1, 1], "t": 0}), + "phys_cond": CartesianDomain({"x": [-1, 1], "t": [0, 1]}), + } + # problem condition statement + conditions = { + "bound_cond1": Condition( + domain="bound_cond1", equation=FixedValue(0.0) + ), + "bound_cond2": Condition( + domain="bound_cond2", equation=FixedValue(0.0) + ), + "time_cond": Condition( + domain="time_cond", equation=Equation(initial_condition) + ), + "phys_cond": Condition( + domain="phys_cond", equation=Equation(burger_equation) + ), + } + + +# +# The `Equation` class takes as input a function (in this case it happens twice, with `initial_condition` and `burger_equation`) which computes a residual of an equation, such as a PDE. In a problem class such as the one above, the `Equation` class with such a given input is passed as a parameter in the specified `Condition`. +# +# The `FixedValue` class takes as input a value of same dimensions of the output functions; this class can be used to enforce a fixed value for a specific condition, e.g. Dirichlet boundary conditions, as it happens for instance in our example. +# +# Once the equations are set as above in the problem conditions, the PINN solver will aim to minimize the residuals described in each equation in the training phase. + +# Available classes of equations include also: +# - `FixedGradient` and `FixedFlux`: they work analogously to `FixedValue` class, where we can require a constant value to be enforced, respectively, on the gradient of the solution or the divergence of the solution; +# - `Laplace`: it can be used to enforce the laplacian of the solution to be zero; +# - `SystemEquation`: we can enforce multiple conditions on the same subdomain through this class, passing a list of residual equations defined in the problem. +# + +# # Defining a new Equation class + +# `Equation` classes can be also inherited to define a new class. As example, we can see how to rewrite the above problem introducing a new class `Burgers1D`; during the class call, we can pass the viscosity parameter $\nu$: + +# In[3]: + + +class Burgers1DEquation(Equation): + + def __init__(self, nu=0.0): + """ + Burgers1D class. This class can be + used to enforce the solution u to solve the viscous Burgers 1D Equation. + + :param torch.float32 nu: the viscosity coefficient. Default value is set to 0. + """ + self.nu = nu + + def equation(input_, output_): + return ( + grad(output_, input_, d="t") + + output_ * grad(output_, input_, d="x") + - self.nu * laplacian(output_, input_, d="x") + ) + + super().__init__(equation) + + +# Now we can just pass the above class as input for the last condition, setting $\nu= \frac{0.01}{\pi}$: + +# In[4]: + + +class Burgers1D(TimeDependentProblem, SpatialProblem): + + # define initial condition + def initial_condition(input_, output_): + u_expected = -torch.sin(torch.pi * input_.extract(["x"])) + return output_.extract(["u"]) - u_expected + + # assign output/ spatial and temporal variables + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [-1, 1]}) + temporal_domain = CartesianDomain({"t": [0, 1]}) + + domains = { + "bound_cond1": CartesianDomain({"x": -1, "t": [0, 1]}), + "bound_cond2": CartesianDomain({"x": 1, "t": [0, 1]}), + "time_cond": CartesianDomain({"x": [-1, 1], "t": 0}), + "phys_cond": CartesianDomain({"x": [-1, 1], "t": [0, 1]}), + } + # problem condition statement + conditions = { + "bound_cond1": Condition( + domain="bound_cond1", equation=FixedValue(0.0) + ), + "bound_cond2": Condition( + domain="bound_cond2", equation=FixedValue(0.0) + ), + "time_cond": Condition( + domain="time_cond", equation=Equation(initial_condition) + ), + "phys_cond": Condition( + domain="phys_cond", equation=Burgers1DEquation(nu=0.01 / torch.pi) + ), + } + + +# # What's next? + +# Congratulations on completing the `Equation` class tutorial of **PINA**! As we have seen, you can build new classes that inherit `Equation` to store more complex equations, as the Burgers 1D equation, only requiring to pass the characteristic coefficients of the problem. +# From now on, you can: +# - define additional complex equation classes (e.g. `SchrodingerEquation`, `NavierStokeEquation`..) +# - define more `FixedOperator` (e.g. `FixedCurl`) diff --git a/tutorials/tutorial13/tutorial.ipynb b/tutorials/tutorial13/tutorial.ipynb index 5bd059d56..9295835a8 100644 --- a/tutorials/tutorial13/tutorial.ipynb +++ b/tutorials/tutorial13/tutorial.ipynb @@ -25,12 +25,13 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", @@ -46,7 +47,7 @@ "from pina.model import FeedForward\n", "from pina.model.block import FourierFeatureEmbedding\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -312,9 +313,13 @@ "l2_loss = LpLoss(p=2, relative=False)\n", "\n", "# sample new test points\n", - "pts = pts = problem.spatial_domain.sample(100, 'grid')\n", - "print(f'Relative l2 error PINN {l2_loss(pinn(pts), problem.solution(pts)).item():.2%}')\n", - "print(f'Relative l2 error SAPINN {l2_loss(sapinn(pts), problem.solution(pts)).item():.2%}')" + "pts = pts = problem.spatial_domain.sample(100, \"grid\")\n", + "print(\n", + " f\"Relative l2 error PINN {l2_loss(pinn(pts), problem.solution(pts)).item():.2%}\"\n", + ")\n", + "print(\n", + " f\"Relative l2 error SAPINN {l2_loss(sapinn(pts), problem.solution(pts)).item():.2%}\"\n", + ")" ] }, { @@ -462,12 +467,14 @@ } ], "source": [ - "#plot solution obtained\n", - "plot_solution(multiscale_pinn, 'Multiscale PINN solution')\n", + "# plot solution obtained\n", + "plot_solution(multiscale_pinn, \"Multiscale PINN solution\")\n", "\n", "# sample new test points\n", - "pts = pts = problem.spatial_domain.sample(100, 'grid')\n", - "print(f'Relative l2 error PINN with MultiscaleFourierNet: {l2_loss(multiscale_pinn(pts), problem.solution(pts)).item():.2%}')" + "pts = pts = problem.spatial_domain.sample(100, \"grid\")\n", + "print(\n", + " f\"Relative l2 error PINN with MultiscaleFourierNet: {l2_loss(multiscale_pinn(pts), problem.solution(pts)).item():.2%}\"\n", + ")" ] }, { diff --git a/tutorials/tutorial13/tutorial.py b/tutorials/tutorial13/tutorial.py new file mode 100644 index 000000000..0ddd6204c --- /dev/null +++ b/tutorials/tutorial13/tutorial.py @@ -0,0 +1,283 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Multiscale PDE learning with Fourier Feature Network +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial13/tutorial.ipynb) +# +# This tutorial presents how to solve with Physics-Informed Neural Networks (PINNs) +# a PDE characterized by multiscale behaviour, as +# presented in [*On the eigenvector bias of Fourier feature networks: From regression to solving +# multi-scale PDEs with physics-informed neural networks*]( +# https://doi.org/10.1016/j.cma.2021.113938). +# +# First of all, some useful imports. + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import matplotlib.pyplot as plt +import warnings + +from pina import Condition, Trainer +from pina.problem import SpatialProblem +from pina.operator import laplacian +from pina.solver import PINN, SelfAdaptivePINN as SAPINN +from pina.loss import LpLoss +from pina.domain import CartesianDomain +from pina.equation import Equation, FixedValue +from pina.model import FeedForward +from pina.model.block import FourierFeatureEmbedding + +warnings.filterwarnings("ignore") + + +# ## Multiscale Problem +# +# We begin by presenting the problem which also can be found in Section 2 of [*On the eigenvector bias of Fourier feature networks: From regression to solving +# multi-scale PDEs with physics-informed neural networks*]( +# https://doi.org/10.1016/j.cma.2021.113938). The one-dimensional Poisson problem we aim to solve is mathematically written as: +# +# \begin{equation} +# \begin{cases} +# \Delta u (x) + f(x) = 0 \quad x \in [0,1], \\ +# u(x) = 0 \quad x \in \partial[0,1], \\ +# \end{cases} +# \end{equation} +# +# We impose the solution as $u(x) = \sin(2\pi x) + 0.1 \sin(50\pi x)$ and obtain the force term $f(x) = (2\pi)^2 \sin(2\pi x) + 0.1 (50 \pi)^2 \sin(50\pi x)$. +# Though this example is simple and pedagogical, it is worth noting that +# the solution exhibits low frequency in the macro-scale and high frequency in the micro-scale, which resembles many +# practical scenarios. +# +# +# In **PINA** this problem is written, as always, as a class [see here for a tutorial on the Problem class](https://mathlab.github.io/PINA/_rst/tutorials/tutorial1/tutorial.html). Below you can find the `Poisson` problem which is mathmatically described above. + +# In[2]: + + +class Poisson(SpatialProblem): + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 1]}) + + def poisson_equation(input_, output_): + x = input_.extract("x") + u_xx = laplacian(output_, input_, components=["u"], d=["x"]) + f = ((2 * torch.pi) ** 2) * torch.sin(2 * torch.pi * x) + 0.1 * ( + (50 * torch.pi) ** 2 + ) * torch.sin(50 * torch.pi * x) + return u_xx + f + + domains = { + "bound_cond0": CartesianDomain({"x": 0.0}), + "bound_cond1": CartesianDomain({"x": 1.0}), + "phys_cond": spatial_domain, + } + # here we write the problem conditions + conditions = { + "bound_cond0": Condition( + domain="bound_cond0", equation=FixedValue(0.0) + ), + "bound_cond1": Condition( + domain="bound_cond1", equation=FixedValue(0.0) + ), + "phys_cond": Condition( + domain="phys_cond", equation=Equation(poisson_equation) + ), + } + + def solution(self, x): + return torch.sin(2 * torch.pi * x) + 0.1 * torch.sin(50 * torch.pi * x) + + +problem = Poisson() + +# let's discretise the domain +problem.discretise_domain(128, "grid", domains=["phys_cond"]) +problem.discretise_domain(1, "grid", domains=["bound_cond0", "bound_cond1"]) + + +# A standard PINN approach would be to fit this model using a Feed Forward (fully connected) Neural Network. For a conventional fully-connected neural network is easy to +# approximate a function $u$, given sufficient data inside the computational domain. However solving high-frequency or multi-scale problems presents great challenges to PINNs especially when the number of data cannot capture the different scales. +# +# Below we run a simulation using the `PINN` solver and the self adaptive `SAPINN` solver, using a [`FeedForward`](https://mathlab.github.io/PINA/_modules/pina/model/feed_forward.html#FeedForward) model. + +# In[3]: + + +# training with PINN and visualize results +pinn = PINN( + problem=problem, + model=FeedForward( + input_dimensions=1, output_dimensions=1, layers=[100, 100, 100] + ), +) + +trainer = Trainer( + pinn, + max_epochs=1500, + accelerator="cpu", + enable_model_summary=False, + val_size=0.0, + train_size=1.0, + test_size=0.0, +) +trainer.train() + +# training with PINN and visualize results +sapinn = SAPINN( + problem=problem, + model=FeedForward( + input_dimensions=1, output_dimensions=1, layers=[100, 100, 100] + ), +) +trainer_sapinn = Trainer( + sapinn, + max_epochs=1500, + accelerator="cpu", + enable_model_summary=False, + val_size=0.0, + train_size=1.0, + test_size=0.0, +) +trainer_sapinn.train() + + +# In[4]: + + +# define the function to plot the solution obtained using matplotlib +def plot_solution(pinn_to_use, title): + pts = pinn_to_use.problem.spatial_domain.sample(256, "grid", variables="x") + predicted_output = pinn_to_use.forward(pts).extract("u").tensor.detach() + true_output = pinn_to_use.problem.solution(pts).detach() + plt.plot( + pts.extract(["x"]), predicted_output, label="Neural Network solution" + ) + plt.plot(pts.extract(["x"]), true_output, label="True solution") + plt.title(title) + plt.legend() + + +# plot the solution of the two PINNs +plot_solution(pinn, "PINN solution") +plt.figure() +plot_solution(sapinn, "Self Adaptive PINN solution") + + +# We can clearly see that the solution has not been learned by the two different solvers. Indeed the big problem is not in the optimization strategy (i.e. the solver), but in the model used to solve the problem. A simple `FeedForward` network can hardly handle multiscales if not enough collocation points are used! +# +# We can also compute the $l_2$ relative error for the `PINN` and `SAPINN` solutions: + +# In[5]: + + +# l2 loss from PINA losses +l2_loss = LpLoss(p=2, relative=False) + +# sample new test points +pts = pts = problem.spatial_domain.sample(100, "grid") +print( + f"Relative l2 error PINN {l2_loss(pinn(pts), problem.solution(pts)).item():.2%}" +) +print( + f"Relative l2 error SAPINN {l2_loss(sapinn(pts), problem.solution(pts)).item():.2%}" +) + + +# Which is indeed very high! + +# ## Fourier Feature Embedding in PINA + +# Fourier Feature Embedding is a way to transform the input features, to help the network in learning multiscale variations in the output. It was +# first introduced in [*On the eigenvector bias of Fourier feature networks: From regression to solving +# multi-scale PDEs with physics-informed neural networks*]( +# https://doi.org/10.1016/j.cma.2021.113938) showing great results for multiscale problems. The basic idea is to map the input $\mathbf{x}$ into an embedding $\tilde{\mathbf{x}}$ where: +# +# $$ \tilde{\mathbf{x}} =\left[\cos\left( \mathbf{B} \mathbf{x} \right), \sin\left( \mathbf{B} \mathbf{x} \right)\right] $$ +# +# and $\mathbf{B}_{ij} \sim \mathcal{N}(0, \sigma^2)$. This simple operation allow the network to learn on multiple scales! +# +# In PINA we already have implemented the feature as a `layer` called [`FourierFeatureEmbedding`](https://mathlab.github.io/PINA/_rst/layers/fourier_embedding.html). Below we will build the *Multi-scale Fourier Feature Architecture*. In this architecture multiple Fourier feature embeddings (initialized with different $\sigma$) +# are applied to input coordinates and then passed through the same fully-connected neural network, before the outputs are finally concatenated with a linear layer. + +# In[6]: + + +class MultiscaleFourierNet(torch.nn.Module): + def __init__(self): + super().__init__() + self.embedding1 = FourierFeatureEmbedding( + input_dimension=1, output_dimension=100, sigma=1 + ) + self.embedding2 = FourierFeatureEmbedding( + input_dimension=1, output_dimension=100, sigma=10 + ) + self.layers = FeedForward( + input_dimensions=100, output_dimensions=100, layers=[100] + ) + self.final_layer = torch.nn.Linear(2 * 100, 1) + + def forward(self, x): + e1 = self.layers(self.embedding1(x)) + e2 = self.layers(self.embedding2(x)) + return self.final_layer(torch.cat([e1, e2], dim=-1)) + + +# We will train the `MultiscaleFourierNet` with the `PINN` solver (and feel free to try also with our PINN variants (`SAPINN`, `GPINN`, `CompetitivePINN`, ...). + +# In[7]: + + +multiscale_pinn = PINN(problem=problem, model=MultiscaleFourierNet()) +trainer = Trainer( + multiscale_pinn, + max_epochs=1500, + accelerator="cpu", + enable_model_summary=False, + val_size=0.0, + train_size=1.0, + test_size=0.0, +) +trainer.train() + + +# Let us now plot the solution and compute the relative $l_2$ again! + +# In[8]: + + +# plot solution obtained +plot_solution(multiscale_pinn, "Multiscale PINN solution") + +# sample new test points +pts = pts = problem.spatial_domain.sample(100, "grid") +print( + f"Relative l2 error PINN with MultiscaleFourierNet: {l2_loss(multiscale_pinn(pts), problem.solution(pts)).item():.2%}" +) + + +# It is pretty clear that the network has learned the correct solution, with also a very low error. Obviously a longer training and a more expressive neural network could improve the results! +# +# ## What's next? +# +# Congratulations on completing the one dimensional Poisson tutorial of **PINA** using `FourierFeatureEmbedding`! There are multiple directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the finaly accuracy +# +# 2. Understand the role of `sigma` in `FourierFeatureEmbedding` (see original paper for a nice reference) +# +# 3. Code the *Spatio-temporal multi-scale Fourier feature architecture* for a more complex time dependent PDE (section 3 of the original reference) +# +# 4. Many more... diff --git a/tutorials/tutorial14/tutorial.ipynb b/tutorials/tutorial14/tutorial.ipynb index 312e725d7..6301722e5 100644 --- a/tutorials/tutorial14/tutorial.ipynb +++ b/tutorials/tutorial14/tutorial.ipynb @@ -33,12 +33,13 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "%matplotlib inline\n", "\n", @@ -50,7 +51,7 @@ "from pina.model.block import PODBlock, RBFBlock\n", "from pina import LabelTensor\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -70,6 +71,7 @@ "source": [ "import smithers\n", "from smithers.dataset import LidCavity\n", + "\n", "dataset = LidCavity()" ] }, @@ -108,13 +110,13 @@ ], "source": [ "fig, axs = plt.subplots(1, 3, figsize=(14, 3))\n", - "for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots['mag(v)'][:3]):\n", + "for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots[\"mag(v)\"][:3]):\n", " ax.tricontourf(dataset.triang, u, levels=16)\n", - " ax.set_title(f'$u$ field for $\\mu$ = {par[0]:.4f}')\n", + " ax.set_title(f\"$u$ field for $\\mu$ = {par[0]:.4f}\")\n", "fig, axs = plt.subplots(1, 3, figsize=(14, 3))\n", - "for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots['p'][:3]):\n", + "for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots[\"p\"][:3]):\n", " ax.tricontourf(dataset.triang, u, levels=16)\n", - " ax.set_title(f'$p$ field for $\\mu$ = {par[0]:.4f}')" + " ax.set_title(f\"$p$ field for $\\mu$ = {par[0]:.4f}\")" ] }, { @@ -130,15 +132,16 @@ "metadata": {}, "outputs": [], "source": [ - "'''velocity magnitude data, 5041 for each snapshot'''\n", - "u=torch.tensor(dataset.snapshots['mag(v)']).float() \n", - "u = LabelTensor(u, labels=[f's{i}' for i in range(u.shape[1])])\n", - "'''pressure data, 5041 for each snapshot'''\n", - "p=torch.tensor(dataset.snapshots['p']).float()\n", - "p = LabelTensor(p, labels=[f's{i}' for i in range(p.shape[1])])\n", - "'''mu corresponding to each snapshot'''\n", - "mu=torch.tensor(dataset.params).float()\n", - "mu = LabelTensor(mu, labels=['mu'])\n" + "\"\"\"velocity magnitude data, 5041 for each snapshot\"\"\"\n", + "\n", + "u = torch.tensor(dataset.snapshots[\"mag(v)\"]).float()\n", + "u = LabelTensor(u, labels=[f\"s{i}\" for i in range(u.shape[1])])\n", + "\"\"\"pressure data, 5041 for each snapshot\"\"\"\n", + "p = torch.tensor(dataset.snapshots[\"p\"]).float()\n", + "p = LabelTensor(p, labels=[f\"s{i}\" for i in range(p.shape[1])])\n", + "\"\"\"mu corresponding to each snapshot\"\"\"\n", + "mu = torch.tensor(dataset.params).float()\n", + "mu = LabelTensor(mu, labels=[\"mu\"])" ] }, { @@ -154,15 +157,16 @@ "metadata": {}, "outputs": [], "source": [ - "'''number of snapshots'''\n", + "\"\"\"number of snapshots\"\"\"\n", + "\n", "n = u.shape[0]\n", - "'''training over total snapshots ratio and number of training snapshots'''\n", - "ratio = 0.9 \n", - "n_train = int(n*ratio)\n", - "'''split u and p data'''\n", - "u_train, u_test = u[:n_train], u[n_train:] #for mag(v)\n", - "p_train, p_test = p[:n_train], p[n_train:] #for p\n", - "'''split snapshots'''\n", + "\"\"\"training over total snapshots ratio and number of training snapshots\"\"\"\n", + "ratio = 0.9\n", + "n_train = int(n * ratio)\n", + "\"\"\"split u and p data\"\"\"\n", + "u_train, u_test = u[:n_train], u[n_train:] # for mag(v)\n", + "p_train, p_test = p[:n_train], p[n_train:] # for p\n", + "\"\"\"split snapshots\"\"\"\n", "mu_train, mu_test = mu[:n_train], mu[n_train:]" ] }, @@ -183,8 +187,9 @@ " \"\"\"\n", " Proper orthogonal decomposition with Radial Basis Function interpolation model.\n", " \"\"\"\n", + "\n", " def __init__(self, pod_rank, rbf_kernel):\n", - " \n", + "\n", " super().__init__()\n", " self.pod = PODBlock(pod_rank)\n", " self.rbf = RBFBlock(kernel=rbf_kernel)" @@ -207,8 +212,9 @@ " \"\"\"\n", " Proper orthogonal decomposition with Radial Basis Function interpolation model.\n", " \"\"\"\n", + "\n", " def __init__(self, pod_rank, rbf_kernel):\n", - " \n", + "\n", " super().__init__()\n", " self.pod = PODBlock(pod_rank)\n", " self.rbf = RBFBlock(kernel=rbf_kernel)\n", @@ -223,6 +229,7 @@ " \"\"\"\n", " coefficients = self.rbf(x)\n", " return self.pod.expand(coefficients)\n", + "\n", " def fit(self, p, x):\n", " \"\"\"\n", " Call the :meth:`pina.model.layers.PODBlock.fit` method of the\n", @@ -231,8 +238,7 @@ " :attr:`pina.model.layers.RBFBlock` attribute to fit the interpolation.\n", " \"\"\"\n", " self.pod.fit(x)\n", - " self.rbf.fit(p, self.pod.reduce(x))\n", - " " + " self.rbf.fit(p, self.pod.reduce(x))" ] }, { @@ -248,15 +254,16 @@ "metadata": {}, "outputs": [], "source": [ - "'''create the model'''\n", - "pod_rbfu = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline')\n", + "\"\"\"create the model\"\"\"\n", "\n", - "'''fit the model to velocity training data'''\n", + "pod_rbfu = PODRBF(pod_rank=20, rbf_kernel=\"thin_plate_spline\")\n", + "\n", + "\"\"\"fit the model to velocity training data\"\"\"\n", "pod_rbfu.fit(mu_train, u_train)\n", "\n", - "'''predict the parameter using the fitted model'''\n", + "\"\"\"predict the parameter using the fitted model\"\"\"\n", "u_train_rbf = pod_rbfu(mu_train)\n", - "u_test_rbf = pod_rbfu(mu_test)\n" + "u_test_rbf = pod_rbfu(mu_test)" ] }, { @@ -282,12 +289,12 @@ } ], "source": [ - "relative_u_error_train = torch.norm(u_train_rbf - u_train)/torch.norm(u_train)\n", - "relative_u_error_test = torch.norm(u_test_rbf - u_test)/torch.norm(u_test)\n", + "relative_u_error_train = torch.norm(u_train_rbf - u_train) / torch.norm(u_train)\n", + "relative_u_error_test = torch.norm(u_test_rbf - u_test) / torch.norm(u_test)\n", "\n", - "print('Error summary for POD-RBF model:')\n", - "print(f' Train: {relative_u_error_train.item():e}')\n", - "print(f' Test: {relative_u_error_test.item():e}')" + "print(\"Error summary for POD-RBF model:\")\n", + "print(f\" Train: {relative_u_error_train.item():e}\")\n", + "print(f\" Test: {relative_u_error_test.item():e}\")" ] }, { @@ -323,23 +330,32 @@ "fig, axs = plt.subplots(3, 4, figsize=(14, 10))\n", "\n", "relative_u_error_rbf = np.abs(u_test[idx] - u_idx_rbf.detach())\n", - "relative_u_error_rbf = np.where(u_test[idx] < 1e-7, 1e-7, relative_u_error_rbf/u_test[idx])\n", - " \n", - "for i, (idx_, rbf_, rbf_err_) in enumerate(\n", - " zip(idx, u_idx_rbf, relative_u_error_rbf)):\n", - " axs[0, i].set_title('Prediction for ' f'$\\mu$ = {mu_test[idx_].item():.4f}')\n", - " axs[1, i].set_title('True snapshot for ' f'$\\mu$ = {mu_test[idx_].item():.4f}')\n", - " axs[2, i].set_title('Error for ' f'$\\mu$ = {mu_test[idx_].item():.4f}')\n", + "relative_u_error_rbf = np.where(\n", + " u_test[idx] < 1e-7, 1e-7, relative_u_error_rbf / u_test[idx]\n", + ")\n", "\n", - " cm = axs[0, i].tricontourf(dataset.triang, rbf_.detach()) # POD-RBF prediction\n", + "for i, (idx_, rbf_, rbf_err_) in enumerate(\n", + " zip(idx, u_idx_rbf, relative_u_error_rbf)\n", + "):\n", + " axs[0, i].set_title(\"Prediction for \" f\"$\\mu$ = {mu_test[idx_].item():.4f}\")\n", + " axs[1, i].set_title(\n", + " \"True snapshot for \" f\"$\\mu$ = {mu_test[idx_].item():.4f}\"\n", + " )\n", + " axs[2, i].set_title(\"Error for \" f\"$\\mu$ = {mu_test[idx_].item():.4f}\")\n", + "\n", + " cm = axs[0, i].tricontourf(\n", + " dataset.triang, rbf_.detach()\n", + " ) # POD-RBF prediction\n", " plt.colorbar(cm, ax=axs[0, i])\n", - " \n", - " cm = axs[1, i].tricontourf(dataset.triang, u_test[idx_].flatten()) # Truth\n", + "\n", + " cm = axs[1, i].tricontourf(dataset.triang, u_test[idx_].flatten()) # Truth\n", " plt.colorbar(cm, ax=axs[1, i])\n", "\n", - " cm = axs[2, i].tripcolor(dataset.triang, rbf_err_, norm=matplotlib.colors.LogNorm()) # Error for POD-RBF\n", + " cm = axs[2, i].tripcolor(\n", + " dataset.triang, rbf_err_, norm=matplotlib.colors.LogNorm()\n", + " ) # Error for POD-RBF\n", " plt.colorbar(cm, ax=axs[2, i])\n", - " \n", + "\n", "plt.show()" ] }, @@ -366,22 +382,23 @@ } ], "source": [ - "'''create the model'''\n", - "pod_rbfp = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline')\n", + "\"\"\"create the model\"\"\"\n", + "\n", + "pod_rbfp = PODRBF(pod_rank=20, rbf_kernel=\"thin_plate_spline\")\n", "\n", - "'''fit the model to pressure training data'''\n", + "\"\"\"fit the model to pressure training data\"\"\"\n", "pod_rbfp.fit(mu_train, p_train)\n", "\n", - "'''predict the parameter using the fitted model'''\n", + "\"\"\"predict the parameter using the fitted model\"\"\"\n", "p_train_rbf = pod_rbfp(mu_train)\n", "p_test_rbf = pod_rbfp(mu_test)\n", "\n", - "relative_p_error_train = torch.norm(p_train_rbf - p_train)/torch.norm(p_train)\n", - "relative_p_error_test = torch.norm(p_test_rbf - p_test)/torch.norm(p_test)\n", + "relative_p_error_train = torch.norm(p_train_rbf - p_train) / torch.norm(p_train)\n", + "relative_p_error_test = torch.norm(p_test_rbf - p_test) / torch.norm(p_test)\n", "\n", - "print('Error summary for POD-RBF model:')\n", - "print(f' Train: {relative_p_error_train.item():e}')\n", - "print(f' Test: {relative_p_error_test.item():e}')" + "print(\"Error summary for POD-RBF model:\")\n", + "print(f\" Train: {relative_p_error_train.item():e}\")\n", + "print(f\" Test: {relative_p_error_test.item():e}\")" ] }, { @@ -409,10 +426,12 @@ ], "source": [ "fig, axs = plt.subplots(2, 3, figsize=(14, 6))\n", - "for ax, par, u in zip(axs.ravel(), dataset.params[66:72], dataset.snapshots['p'][66:72]):\n", + "for ax, par, u in zip(\n", + " axs.ravel(), dataset.params[66:72], dataset.snapshots[\"p\"][66:72]\n", + "):\n", " cm = ax.tricontourf(dataset.triang, u, levels=16)\n", " plt.colorbar(cm, ax=ax)\n", - " ax.set_title(f'$p$ field for $\\mu$ = {par[0]:.4f}')\n", + " ax.set_title(f\"$p$ field for $\\mu$ = {par[0]:.4f}\")\n", "plt.tight_layout()\n", "plt.show()" ] @@ -442,11 +461,13 @@ ], "source": [ "fig, axs = plt.subplots(2, 3, figsize=(14, 6))\n", - "for ax, par, u in zip(axs.ravel(), dataset.params[98:104], dataset.snapshots['p'][98:104]):\n", + "for ax, par, u in zip(\n", + " axs.ravel(), dataset.params[98:104], dataset.snapshots[\"p\"][98:104]\n", + "):\n", " cm = ax.tricontourf(dataset.triang, u, levels=16)\n", " plt.colorbar(cm, ax=ax)\n", - " ax.set_title(f'$p$ field for $\\mu$ = {par[0]:.4f}')\n", - "plt.tight_layout() \n", + " ax.set_title(f\"$p$ field for $\\mu$ = {par[0]:.4f}\")\n", + "plt.tight_layout()\n", "plt.show()" ] }, @@ -473,37 +494,42 @@ } ], "source": [ - "'''excluding problematic snapshots'''\n", + "\"\"\"excluding problematic snapshots\"\"\"\n", + "\n", "data = list(range(300))\n", "data_to_consider = data[:67] + data[71:100] + data[102:]\n", - "'''proceed as before'''\n", - "newp=torch.tensor(dataset.snapshots['p'][data_to_consider]).float()\n", - "newp = LabelTensor(newp, labels=[f's{i}' for i in range(newp.shape[1])])\n", + "\"\"\"proceed as before\"\"\"\n", + "newp = torch.tensor(dataset.snapshots[\"p\"][data_to_consider]).float()\n", + "newp = LabelTensor(newp, labels=[f\"s{i}\" for i in range(newp.shape[1])])\n", "\n", - "newmu=torch.tensor(dataset.params[data_to_consider]).float()\n", - "newmu = LabelTensor(newmu, labels=['mu'])\n", + "newmu = torch.tensor(dataset.params[data_to_consider]).float()\n", + "newmu = LabelTensor(newmu, labels=[\"mu\"])\n", "\n", "newn = newp.shape[0]\n", - "ratio = 0.9 \n", - "new_train = int(newn*ratio)\n", + "ratio = 0.9\n", + "new_train = int(newn * ratio)\n", "\n", - "new_p_train, new_p_test = newp[:new_train], newp[new_train:] \n", + "new_p_train, new_p_test = newp[:new_train], newp[new_train:]\n", "\n", "new_mu_train, new_mu_test = newmu[:new_train], newmu[new_train:]\n", "\n", - "new_pod_rbfp = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline')\n", + "new_pod_rbfp = PODRBF(pod_rank=20, rbf_kernel=\"thin_plate_spline\")\n", "\n", "new_pod_rbfp.fit(new_mu_train, new_p_train)\n", "\n", "new_p_train_rbf = new_pod_rbfp(new_mu_train)\n", "new_p_test_rbf = new_pod_rbfp(new_mu_test)\n", "\n", - "new_relative_p_error_train = torch.norm(new_p_train_rbf - new_p_train)/torch.norm(new_p_train)\n", - "new_relative_p_error_test = torch.norm(new_p_test_rbf - new_p_test)/torch.norm(new_p_test)\n", + "new_relative_p_error_train = torch.norm(\n", + " new_p_train_rbf - new_p_train\n", + ") / torch.norm(new_p_train)\n", + "new_relative_p_error_test = torch.norm(\n", + " new_p_test_rbf - new_p_test\n", + ") / torch.norm(new_p_test)\n", "\n", - "print('Error summary for POD-RBF model:')\n", - "print(f' Train: {new_relative_p_error_train.item():e}')\n", - "print(f' Test: {new_relative_p_error_test.item():e}')" + "print(\"Error summary for POD-RBF model:\")\n", + "print(f\" Train: {new_relative_p_error_train.item():e}\")\n", + "print(f\" Test: {new_relative_p_error_test.item():e}\")" ] }, { diff --git a/tutorials/tutorial14/tutorial.py b/tutorials/tutorial14/tutorial.py index 97c0bccff..f60a598c3 100644 --- a/tutorials/tutorial14/tutorial.py +++ b/tutorials/tutorial14/tutorial.py @@ -2,11 +2,11 @@ # coding: utf-8 # # Tutorial: Predicting Lid-driven cavity problem parameters with POD-RBF -# +# # [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial14/tutorial.ipynb) -# In this tutorial we will show how to use the **PINA** library to predict the distributions of velocity and pressure the Lid-driven Cavity problem, a benchmark in Computational Fluid Dynamics. The problem consists of a square cavity with a lid on top moving with tangential velocity (by convention to the right), with the addition of no-slip conditions on the walls of the cavity and null static pressure on the lower left angle. -# +# In this tutorial we will show how to use the **PINA** library to predict the distributions of velocity and pressure the Lid-driven Cavity problem, a benchmark in Computational Fluid Dynamics. The problem consists of a square cavity with a lid on top moving with tangential velocity (by convention to the right), with the addition of no-slip conditions on the walls of the cavity and null static pressure on the lower left angle. +# # Our goal is to predict the distributions of velocity and pressure of the fluid inside the cavity as the Reynolds number of the inlet fluid varies. To do so we're using a Reduced Order Model (ROM) based on Proper Orthogonal Decomposition (POD). The parametric solution manifold is approximated here with Radial Basis Function (RBF) Interpolation, a common mesh-free interpolation method that doesn't require trainers or solvers as the found radial basis functions are used to interpolate new points. # Let's start with the necessary imports. We're particularly interested in the `PODBlock` and `RBFBlock` classes which will allow us to define the POD-RBF model. @@ -16,14 +16,15 @@ ## routine needed to run the notebook on Google Colab try: - import google.colab - IN_COLAB = True + import google.colab + + IN_COLAB = True except: - IN_COLAB = False + IN_COLAB = False if IN_COLAB: - get_ipython().system('pip install "pina-mathlab"') + get_ipython().system('pip install "pina-mathlab"') -get_ipython().run_line_magic('matplotlib', 'inline') +get_ipython().run_line_magic("matplotlib", "inline") import matplotlib.pyplot as plt import torch @@ -33,11 +34,11 @@ from pina.model.block import PODBlock, RBFBlock from pina import LabelTensor -warnings.filterwarnings('ignore') +warnings.filterwarnings("ignore") # In this tutorial we're gonna use the `LidCavity` class from the [Smithers](https://github.com/mathLab/Smithers) library, which contains a set of parametric solutions of the Lid-driven cavity problem in a square domain. The dataset consists of 300 snapshots of the parameter fields, which in this case are the magnitude of velocity and the pressure, and the corresponding parameter values $u$ and $p$. Each snapshot corresponds to a different value of the tangential velocity $\mu$ of the lid, which has been sampled uniformly between 0.01 m/s and 1 m/s. -# +# # Let's start by importing the dataset: # In[2]: @@ -45,6 +46,7 @@ import smithers from smithers.dataset import LidCavity + dataset = LidCavity() @@ -54,13 +56,13 @@ fig, axs = plt.subplots(1, 3, figsize=(14, 3)) -for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots['mag(v)'][:3]): +for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots["mag(v)"][:3]): ax.tricontourf(dataset.triang, u, levels=16) - ax.set_title(f'$u$ field for $\mu$ = {par[0]:.4f}') + ax.set_title(f"$u$ field for $\mu$ = {par[0]:.4f}") fig, axs = plt.subplots(1, 3, figsize=(14, 3)) -for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots['p'][:3]): +for ax, par, u in zip(axs, dataset.params[:3], dataset.snapshots["p"][:3]): ax.tricontourf(dataset.triang, u, levels=16) - ax.set_title(f'$p$ field for $\mu$ = {par[0]:.4f}') + ax.set_title(f"$p$ field for $\mu$ = {par[0]:.4f}") # To train the model we only need the snapshots for the two parameters. In order to be able to work with the snapshots in **PINA** we first need to assure they're in a compatible format, hence why we start by casting them into `LabelTensor` objects: @@ -68,15 +70,16 @@ # In[4]: -'''velocity magnitude data, 5041 for each snapshot''' -u=torch.tensor(dataset.snapshots['mag(v)']).float() -u = LabelTensor(u, labels=[f's{i}' for i in range(u.shape[1])]) -'''pressure data, 5041 for each snapshot''' -p=torch.tensor(dataset.snapshots['p']).float() -p = LabelTensor(p, labels=[f's{i}' for i in range(p.shape[1])]) -'''mu corresponding to each snapshot''' -mu=torch.tensor(dataset.params).float() -mu = LabelTensor(mu, labels=['mu']) +"""velocity magnitude data, 5041 for each snapshot""" + +u = torch.tensor(dataset.snapshots["mag(v)"]).float() +u = LabelTensor(u, labels=[f"s{i}" for i in range(u.shape[1])]) +"""pressure data, 5041 for each snapshot""" +p = torch.tensor(dataset.snapshots["p"]).float() +p = LabelTensor(p, labels=[f"s{i}" for i in range(p.shape[1])]) +"""mu corresponding to each snapshot""" +mu = torch.tensor(dataset.params).float() +mu = LabelTensor(mu, labels=["mu"]) # The goal of our training is to be able to predict the solution for new test parameters. The first thing we need to do is validate the accuracy of the model, and in order to do so we split the 300 snapshots in training and testing dataset. In the example we set the training `ratio` to 0.9, which means that 90% of the total snapshots is used for training and the remaining 10% for testing. @@ -84,15 +87,16 @@ # In[5]: -'''number of snapshots''' +"""number of snapshots""" + n = u.shape[0] -'''training over total snapshots ratio and number of training snapshots''' -ratio = 0.9 -n_train = int(n*ratio) -'''split u and p data''' -u_train, u_test = u[:n_train], u[n_train:] #for mag(v) -p_train, p_test = p[:n_train], p[n_train:] #for p -'''split snapshots''' +"""training over total snapshots ratio and number of training snapshots""" +ratio = 0.9 +n_train = int(n * ratio) +"""split u and p data""" +u_train, u_test = u[:n_train], u[n_train:] # for mag(v) +p_train, p_test = p[:n_train], p[n_train:] # for p +"""split snapshots""" mu_train, mu_test = mu[:n_train], mu[n_train:] @@ -105,8 +109,9 @@ class PODRBF(torch.nn.Module): """ Proper orthogonal decomposition with Radial Basis Function interpolation model. """ + def __init__(self, pod_rank, rbf_kernel): - + super().__init__() self.pod = PODBlock(pod_rank) self.rbf = RBFBlock(kernel=rbf_kernel) @@ -121,8 +126,9 @@ class PODRBF(torch.nn.Module): """ Proper orthogonal decomposition with Radial Basis Function interpolation model. """ + def __init__(self, pod_rank, rbf_kernel): - + super().__init__() self.pod = PODBlock(pod_rank) self.rbf = RBFBlock(kernel=rbf_kernel) @@ -137,6 +143,7 @@ def forward(self, x): """ coefficients = self.rbf(x) return self.pod.expand(coefficients) + def fit(self, p, x): """ Call the :meth:`pina.model.layers.PODBlock.fit` method of the @@ -146,7 +153,6 @@ def fit(self, p, x): """ self.pod.fit(x) self.rbf.fit(p, self.pod.reduce(x)) - # Now that we've built our class, we can fit the model and ask it to predict the parameters for the remaining snapshots. We remember that we don't need to train the model, as it doesn't involve any learnable parameter. The only things we have to set are the rank of the decomposition and the radial basis function (here we use thin plate). Here we focus on predicting the magnitude of velocity: @@ -154,13 +160,14 @@ def fit(self, p, x): # In[8]: -'''create the model''' -pod_rbfu = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline') +"""create the model""" -'''fit the model to velocity training data''' +pod_rbfu = PODRBF(pod_rank=20, rbf_kernel="thin_plate_spline") + +"""fit the model to velocity training data""" pod_rbfu.fit(mu_train, u_train) -'''predict the parameter using the fitted model''' +"""predict the parameter using the fitted model""" u_train_rbf = pod_rbfu(mu_train) u_test_rbf = pod_rbfu(mu_test) @@ -170,12 +177,12 @@ def fit(self, p, x): # In[9]: -relative_u_error_train = torch.norm(u_train_rbf - u_train)/torch.norm(u_train) -relative_u_error_test = torch.norm(u_test_rbf - u_test)/torch.norm(u_test) +relative_u_error_train = torch.norm(u_train_rbf - u_train) / torch.norm(u_train) +relative_u_error_test = torch.norm(u_test_rbf - u_test) / torch.norm(u_test) -print('Error summary for POD-RBF model:') -print(f' Train: {relative_u_error_train.item():e}') -print(f' Test: {relative_u_error_test.item():e}') +print("Error summary for POD-RBF model:") +print(f" Train: {relative_u_error_train.item():e}") +print(f" Test: {relative_u_error_test.item():e}") # The results are promising! Now let's visualise them, comparing four random predicted snapshots to the true ones: @@ -192,23 +199,32 @@ def fit(self, p, x): fig, axs = plt.subplots(3, 4, figsize=(14, 10)) relative_u_error_rbf = np.abs(u_test[idx] - u_idx_rbf.detach()) -relative_u_error_rbf = np.where(u_test[idx] < 1e-7, 1e-7, relative_u_error_rbf/u_test[idx]) - -for i, (idx_, rbf_, rbf_err_) in enumerate( - zip(idx, u_idx_rbf, relative_u_error_rbf)): - axs[0, i].set_title('Prediction for ' f'$\mu$ = {mu_test[idx_].item():.4f}') - axs[1, i].set_title('True snapshot for ' f'$\mu$ = {mu_test[idx_].item():.4f}') - axs[2, i].set_title('Error for ' f'$\mu$ = {mu_test[idx_].item():.4f}') +relative_u_error_rbf = np.where( + u_test[idx] < 1e-7, 1e-7, relative_u_error_rbf / u_test[idx] +) - cm = axs[0, i].tricontourf(dataset.triang, rbf_.detach()) # POD-RBF prediction +for i, (idx_, rbf_, rbf_err_) in enumerate( + zip(idx, u_idx_rbf, relative_u_error_rbf) +): + axs[0, i].set_title("Prediction for " f"$\mu$ = {mu_test[idx_].item():.4f}") + axs[1, i].set_title( + "True snapshot for " f"$\mu$ = {mu_test[idx_].item():.4f}" + ) + axs[2, i].set_title("Error for " f"$\mu$ = {mu_test[idx_].item():.4f}") + + cm = axs[0, i].tricontourf( + dataset.triang, rbf_.detach() + ) # POD-RBF prediction plt.colorbar(cm, ax=axs[0, i]) - - cm = axs[1, i].tricontourf(dataset.triang, u_test[idx_].flatten()) # Truth + + cm = axs[1, i].tricontourf(dataset.triang, u_test[idx_].flatten()) # Truth plt.colorbar(cm, ax=axs[1, i]) - cm = axs[2, i].tripcolor(dataset.triang, rbf_err_, norm=matplotlib.colors.LogNorm()) # Error for POD-RBF + cm = axs[2, i].tripcolor( + dataset.triang, rbf_err_, norm=matplotlib.colors.LogNorm() + ) # Error for POD-RBF plt.colorbar(cm, ax=axs[2, i]) - + plt.show() @@ -217,34 +233,37 @@ def fit(self, p, x): # In[11]: -'''create the model''' -pod_rbfp = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline') +"""create the model""" + +pod_rbfp = PODRBF(pod_rank=20, rbf_kernel="thin_plate_spline") -'''fit the model to pressure training data''' +"""fit the model to pressure training data""" pod_rbfp.fit(mu_train, p_train) -'''predict the parameter using the fitted model''' +"""predict the parameter using the fitted model""" p_train_rbf = pod_rbfp(mu_train) p_test_rbf = pod_rbfp(mu_test) -relative_p_error_train = torch.norm(p_train_rbf - p_train)/torch.norm(p_train) -relative_p_error_test = torch.norm(p_test_rbf - p_test)/torch.norm(p_test) +relative_p_error_train = torch.norm(p_train_rbf - p_train) / torch.norm(p_train) +relative_p_error_test = torch.norm(p_test_rbf - p_test) / torch.norm(p_test) -print('Error summary for POD-RBF model:') -print(f' Train: {relative_p_error_train.item():e}') -print(f' Test: {relative_p_error_test.item():e}') +print("Error summary for POD-RBF model:") +print(f" Train: {relative_p_error_train.item():e}") +print(f" Test: {relative_p_error_test.item():e}") -# Unfortunately here we obtain a very high relative test error, although this is likely due to the nature of the available data. Looking at the plots we can see that the pressure field is subject to high variations between subsequent snapshots, especially here: +# Unfortunately here we obtain a very high relative test error, although this is likely due to the nature of the available data. Looking at the plots we can see that the pressure field is subject to high variations between subsequent snapshots, especially here: # In[12]: fig, axs = plt.subplots(2, 3, figsize=(14, 6)) -for ax, par, u in zip(axs.ravel(), dataset.params[66:72], dataset.snapshots['p'][66:72]): +for ax, par, u in zip( + axs.ravel(), dataset.params[66:72], dataset.snapshots["p"][66:72] +): cm = ax.tricontourf(dataset.triang, u, levels=16) plt.colorbar(cm, ax=ax) - ax.set_title(f'$p$ field for $\mu$ = {par[0]:.4f}') + ax.set_title(f"$p$ field for $\mu$ = {par[0]:.4f}") plt.tight_layout() plt.show() @@ -255,11 +274,13 @@ def fit(self, p, x): fig, axs = plt.subplots(2, 3, figsize=(14, 6)) -for ax, par, u in zip(axs.ravel(), dataset.params[98:104], dataset.snapshots['p'][98:104]): +for ax, par, u in zip( + axs.ravel(), dataset.params[98:104], dataset.snapshots["p"][98:104] +): cm = ax.tricontourf(dataset.triang, u, levels=16) plt.colorbar(cm, ax=ax) - ax.set_title(f'$p$ field for $\mu$ = {par[0]:.4f}') -plt.tight_layout() + ax.set_title(f"$p$ field for $\mu$ = {par[0]:.4f}") +plt.tight_layout() plt.show() @@ -268,45 +289,50 @@ def fit(self, p, x): # In[14]: -'''excluding problematic snapshots''' +"""excluding problematic snapshots""" + data = list(range(300)) data_to_consider = data[:67] + data[71:100] + data[102:] -'''proceed as before''' -newp=torch.tensor(dataset.snapshots['p'][data_to_consider]).float() -newp = LabelTensor(newp, labels=[f's{i}' for i in range(newp.shape[1])]) +"""proceed as before""" +newp = torch.tensor(dataset.snapshots["p"][data_to_consider]).float() +newp = LabelTensor(newp, labels=[f"s{i}" for i in range(newp.shape[1])]) -newmu=torch.tensor(dataset.params[data_to_consider]).float() -newmu = LabelTensor(newmu, labels=['mu']) +newmu = torch.tensor(dataset.params[data_to_consider]).float() +newmu = LabelTensor(newmu, labels=["mu"]) newn = newp.shape[0] -ratio = 0.9 -new_train = int(newn*ratio) +ratio = 0.9 +new_train = int(newn * ratio) -new_p_train, new_p_test = newp[:new_train], newp[new_train:] +new_p_train, new_p_test = newp[:new_train], newp[new_train:] new_mu_train, new_mu_test = newmu[:new_train], newmu[new_train:] -new_pod_rbfp = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline') +new_pod_rbfp = PODRBF(pod_rank=20, rbf_kernel="thin_plate_spline") new_pod_rbfp.fit(new_mu_train, new_p_train) new_p_train_rbf = new_pod_rbfp(new_mu_train) new_p_test_rbf = new_pod_rbfp(new_mu_test) -new_relative_p_error_train = torch.norm(new_p_train_rbf - new_p_train)/torch.norm(new_p_train) -new_relative_p_error_test = torch.norm(new_p_test_rbf - new_p_test)/torch.norm(new_p_test) +new_relative_p_error_train = torch.norm( + new_p_train_rbf - new_p_train +) / torch.norm(new_p_train) +new_relative_p_error_test = torch.norm( + new_p_test_rbf - new_p_test +) / torch.norm(new_p_test) -print('Error summary for POD-RBF model:') -print(f' Train: {new_relative_p_error_train.item():e}') -print(f' Test: {new_relative_p_error_test.item():e}') +print("Error summary for POD-RBF model:") +print(f" Train: {new_relative_p_error_train.item():e}") +print(f" Test: {new_relative_p_error_test.item():e}") # ## What's next? -# +# # Congratulations on completing the **PINA** tutorial on building and using a custom POD class! Now you can try: -# +# # 1. Varying the inputs of the model (for a list of the supported RB functions look at the `rbf_layer.py` file in `pina.layers`) -# +# # 2. Changing the POD model, for example using Artificial Neural Networks. For a more in depth overview of POD-NN and a comparison with the POD-RBF model already shown, look at [Tutorial: Reduced order model (POD-RBF or POD-NN) for parametric problems](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial9/tutorial.ipynb) -# +# # 3. Building your own classes or adapt the one shown to other datasets/problems diff --git a/tutorials/tutorial2/tutorial.ipynb b/tutorials/tutorial2/tutorial.ipynb index 6a1a86282..d0d891c77 100644 --- a/tutorials/tutorial2/tutorial.ipynb +++ b/tutorials/tutorial2/tutorial.ipynb @@ -23,12 +23,13 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", @@ -39,7 +40,7 @@ "from pina.solver import PINN\n", "from torch.nn import Softplus\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -186,8 +187,7 @@ " train_size=0.8, # set train size\n", " val_size=0.0, # set validation size\n", " test_size=0.2, # set testing size\n", - " shuffle=True, # shuffle the data\n", - " \n", + " shuffle=True, # shuffle the data\n", ")\n", "\n", "# train\n", @@ -348,7 +348,7 @@ "\n", " def forward(self, pts):\n", " x, y = pts.extract([\"x\"]), pts.extract([\"y\"])\n", - " f = 2 *torch.pi**2 * torch.sin(x * torch.pi) * torch.sin(y * torch.pi)\n", + " f = 2 * torch.pi**2 * torch.sin(x * torch.pi) * torch.sin(y * torch.pi)\n", " return LabelTensor(f, [\"feat\"])\n", "\n", "\n", @@ -487,28 +487,36 @@ "source": [ "class SinSinAB(torch.nn.Module):\n", " \"\"\" \"\"\"\n", + "\n", " def __init__(self):\n", " super().__init__()\n", " self.alpha = torch.nn.Parameter(torch.tensor([1.0]))\n", " self.beta = torch.nn.Parameter(torch.tensor([1.0]))\n", "\n", " def forward(self, x):\n", - " t = (\n", - " self.beta*torch.sin(self.alpha*x.extract(['x'])*torch.pi)*\n", - " torch.sin(self.alpha*x.extract(['y'])*torch.pi)\n", + " t = (\n", + " self.beta\n", + " * torch.sin(self.alpha * x.extract([\"x\"]) * torch.pi)\n", + " * torch.sin(self.alpha * x.extract([\"y\"]) * torch.pi)\n", " )\n", - " return LabelTensor(t, ['b*sin(a*x)sin(a*y)'])\n", + " return LabelTensor(t, [\"b*sin(a*x)sin(a*y)\"])\n", "\n", "\n", "# make model + solver + trainer\n", "model_learn = FeedForwardWithExtraFeatures(\n", - " input_dimensions=len(problem.input_variables) + 1, #we add one as also we consider the extra feature dimension\n", + " input_dimensions=len(problem.input_variables)\n", + " + 1, # we add one as also we consider the extra feature dimension\n", " output_dimensions=len(problem.output_variables),\n", " func=Softplus,\n", " layers=[10, 10],\n", - " extra_features=SinSinAB())\n", + " extra_features=SinSinAB(),\n", + ")\n", "\n", - "pinn_learn = PINN(problem, model_learn, optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006,weight_decay=1e-8))\n", + "pinn_learn = PINN(\n", + " problem,\n", + " model_learn,\n", + " optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006, weight_decay=1e-8),\n", + ")\n", "trainer_learn = Trainer(\n", " solver=pinn_learn, # setting the solver, i.e. PINN\n", " max_epochs=1000, # setting max epochs in training\n", @@ -649,14 +657,14 @@ ], "source": [ "# test error base pinn\n", - "print('PINN')\n", + "print(\"PINN\")\n", "trainer_base.test()\n", "# test error extra features pinn\n", "print(\"PINN with extra features\")\n", "trainer_feat.test()\n", "# test error learnable extra features pinn\n", "print(\"PINN with learnable extra features\")\n", - "_=trainer_learn.test()" + "_ = trainer_learn.test()" ] }, { diff --git a/tutorials/tutorial2/tutorial.py b/tutorials/tutorial2/tutorial.py new file mode 100644 index 000000000..622783aaa --- /dev/null +++ b/tutorials/tutorial2/tutorial.py @@ -0,0 +1,360 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Two dimensional Poisson problem using Extra Features Learning +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial2/tutorial.ipynb) +# +# This tutorial presents how to solve with Physics-Informed Neural Networks (PINNs) a 2D Poisson problem with Dirichlet boundary conditions. We will train with standard PINN's training, and with extrafeatures. For more insights on extrafeature learning please read [*An extended physics informed neural network for preliminary analysis of parametric optimal control problems*](https://www.sciencedirect.com/science/article/abs/pii/S0898122123002018). +# +# First of all, some useful imports. + +# In[ ]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import matplotlib.pyplot as plt +import warnings + +from pina import LabelTensor, Trainer +from pina.model import FeedForward +from pina.solver import PINN +from torch.nn import Softplus + +warnings.filterwarnings("ignore") + + +# ## The problem definition + +# The two-dimensional Poisson problem is mathematically written as: +# \begin{equation} +# \begin{cases} +# \Delta u = 2\pi^2\sin{(\pi x)} \sin{(\pi y)} \text{ in } D, \\ +# u = 0 \text{ on } \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4, +# \end{cases} +# \end{equation} +# where $D$ is a square domain $[0,1]^2$, and $\Gamma_i$, with $i=1,...,4$, are the boundaries of the square. +# +# The Poisson problem is written in **PINA** code as a class. The equations are written as *conditions* that should be satisfied in the corresponding domains. The *solution* +# is the exact solution which will be compared with the predicted one. If interested in how to write problems see [this tutorial](https://mathlab.github.io/PINA/_rst/tutorials/tutorial1/tutorial.html). +# +# We will directly import the problem from `pina.problem.zoo`, which contains a vast list of PINN problems and more. + +# In[2]: + + +from pina.problem.zoo import Poisson2DSquareProblem as Poisson + +# initialize the problem +problem = Poisson() + +# print the conditions +print( + f"The problem is made of {len(problem.conditions.keys())} conditions: \n" + f"They are: {list(problem.conditions.keys())}" +) + +# let's discretise the domain +problem.discretise_domain(30, "grid", domains=["D"]) +problem.discretise_domain( + 100, + "grid", + domains=["g1", "g2", "g3", "g4"], +) + + +# ## Solving the problem with standard PINNs + +# After the problem, the feed-forward neural network is defined, through the class `FeedForward`. This neural network takes as input the coordinates (in this case $x$ and $y$) and provides the unkwown field of the Poisson problem. The residual of the equations are evaluated at several sampling points and the loss minimized by the neural network is the sum of the residuals. +# +# In this tutorial, the neural network is composed by two hidden layers of 10 neurons each, and it is trained for 1000 epochs with a learning rate of 0.006 and $l_2$ weight regularization set to $10^{-8}$. These parameters can be modified as desired. We set the `train_size` to 0.8 and `test_size` to 0.2, this mean that the discretised points will be divided in a 80%-20% fashion, where 80% will be used for training and the remaining 20% for testing. + +# In[3]: + + +# make model + solver + trainer +from pina.optim import TorchOptimizer + +model = FeedForward( + layers=[10, 10], + func=Softplus, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) +pinn = PINN( + problem, + model, + optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006, weight_decay=1e-8), +) +trainer_base = Trainer( + solver=pinn, # setting the solver, i.e. PINN + max_epochs=1000, # setting max epochs in training + accelerator="cpu", # we train on cpu, also other are available + enable_model_summary=False, # model summary statistics not printed + train_size=0.8, # set train size + val_size=0.0, # set validation size + test_size=0.2, # set testing size + shuffle=True, # shuffle the data +) + +# train +trainer_base.train() + + +# Now we plot the results using `matplotlib`. +# The solution predicted by the neural network is plotted on the left, the exact one is represented at the center and on the right the error between the exact and the predicted solutions is showed. + +# In[4]: + + +@torch.no_grad() +def plot_solution(solver): + # get the problem + problem = solver.problem + # get spatial points + spatial_samples = problem.spatial_domain.sample(30, "grid") + # compute pinn solution, true solution and absolute difference + data = { + "PINN solution": solver(spatial_samples), + "True solution": problem.solution(spatial_samples), + "Absolute Difference": torch.abs( + solver(spatial_samples) - problem.solution(spatial_samples) + ), + } + # plot the solution + for idx, (title, field) in enumerate(data.items()): + plt.subplot(1, 3, idx + 1) + plt.title(title) + plt.tricontourf( # convert to torch tensor + flatten + spatial_samples.extract("x").tensor.flatten(), + spatial_samples.extract("y").tensor.flatten(), + field.tensor.flatten(), + ) + plt.colorbar(), plt.tight_layout() + + +# Here the solution: + +# In[5]: + + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn) + + +# As you can see the solution is not very accurate, in what follows we will use **Extra Feature** as introduced in [*An extended physics informed neural network for preliminary analysis of parametric optimal control problems*](https://www.sciencedirect.com/science/article/abs/pii/S0898122123002018) to boost the training accuracy. Of course, even extra training will benefit, this tutorial is just to show that convergence using Extra Features is usally faster. + +# ## Solving the problem with extra-features PINNs + +# Now, the same problem is solved in a different way. +# A new neural network is now defined, with an additional input variable, named extra-feature, which coincides with the forcing term in the Laplace equation. +# The set of input variables to the neural network is: +# +# \begin{equation} +# [x, y, k(x, y)], \text{ with } k(x, y)= 2\pi^2\sin{(\pi x)}\sin{(\pi y)}, +# \end{equation} +# +# where $x$ and $y$ are the spatial coordinates and $k(x, y)$ is the added feature which is equal to the forcing term. +# +# This feature is initialized in the class `SinSin`, which is a simple `torch.nn.Module`. After declaring such feature, we can just adjust the `FeedForward` class by creating a subclass `FeedForwardWithExtraFeatures` with an adjusted forward method and the additional attribute `extra_features`. +# +# Finally, we perform the same training as before: the problem is `Poisson`, the network is composed by the same number of neurons and optimizer parameters are equal to previous test, the only change is the new extra feature. + +# In[6]: + + +class SinSin(torch.nn.Module): + """Feature: sin(x)*sin(y)""" + + def __init__(self): + super().__init__() + + def forward(self, pts): + x, y = pts.extract(["x"]), pts.extract(["y"]) + f = 2 * torch.pi**2 * torch.sin(x * torch.pi) * torch.sin(y * torch.pi) + return LabelTensor(f, ["feat"]) + + +class FeedForwardWithExtraFeatures(FeedForward): + def __init__(self, *args, extra_features, **kwargs): + super().__init__(*args, **kwargs) + self.extra_features = extra_features + + def forward(self, x): + extra_feature = self.extra_features(x) # we append extra features + x = x.append(extra_feature) + return super().forward(x) + + +model_feat = FeedForwardWithExtraFeatures( + input_dimensions=len(problem.input_variables) + 1, + output_dimensions=len(problem.output_variables), + func=Softplus, + layers=[10, 10], + extra_features=SinSin(), +) + +pinn_feat = PINN( + problem, + model_feat, + optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006, weight_decay=1e-8), +) +trainer_feat = Trainer( + solver=pinn_feat, # setting the solver, i.e. PINN + max_epochs=1000, # setting max epochs in training + accelerator="cpu", # we train on cpu, also other are available + enable_model_summary=False, # model summary statistics not printed + train_size=0.8, # set train size + val_size=0.0, # set validation size + test_size=0.2, # set testing size + shuffle=True, # shuffle the data +) + +trainer_feat.train() + + +# The predicted and exact solutions and the error between them are represented below. +# We can easily note that now our network, having almost the same condition as before, is able to reach additional order of magnitudes in accuracy. + +# In[7]: + + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn_feat) + + +# ## Solving the problem with learnable extra-features PINNs + +# We can still do better! +# +# Another way to exploit the extra features is the addition of learnable parameter inside them. +# In this way, the added parameters are learned during the training phase of the neural network. In this case, we use: +# +# \begin{equation} +# k(x, \mathbf{y}) = \beta \sin{(\alpha x)} \sin{(\alpha y)}, +# \end{equation} +# +# where $\alpha$ and $\beta$ are the abovementioned parameters. +# Their implementation is quite trivial: by using the class `torch.nn.Parameter` we cam define all the learnable parameters we need, and they are managed by `autograd` module! + +# In[8]: + + +class SinSinAB(torch.nn.Module): + """ """ + + def __init__(self): + super().__init__() + self.alpha = torch.nn.Parameter(torch.tensor([1.0])) + self.beta = torch.nn.Parameter(torch.tensor([1.0])) + + def forward(self, x): + t = ( + self.beta + * torch.sin(self.alpha * x.extract(["x"]) * torch.pi) + * torch.sin(self.alpha * x.extract(["y"]) * torch.pi) + ) + return LabelTensor(t, ["b*sin(a*x)sin(a*y)"]) + + +# make model + solver + trainer +model_learn = FeedForwardWithExtraFeatures( + input_dimensions=len(problem.input_variables) + + 1, # we add one as also we consider the extra feature dimension + output_dimensions=len(problem.output_variables), + func=Softplus, + layers=[10, 10], + extra_features=SinSinAB(), +) + +pinn_learn = PINN( + problem, + model_learn, + optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006, weight_decay=1e-8), +) +trainer_learn = Trainer( + solver=pinn_learn, # setting the solver, i.e. PINN + max_epochs=1000, # setting max epochs in training + accelerator="cpu", # we train on cpu, also other are available + enable_model_summary=False, # model summary statistics not printed + train_size=0.8, # set train size + val_size=0.0, # set validation size + test_size=0.2, # set testing size + shuffle=True, # shuffle the data +) +# train +trainer_learn.train() + + +# Umh, the final loss is not appreciabily better than previous model (with static extra features), despite the usage of learnable parameters. This is mainly due to the over-parametrization of the network: there are many parameter to optimize during the training, and the model in unable to understand automatically that only the parameters of the extra feature (and not the weights/bias of the FFN) should be tuned in order to fit our problem. A longer training can be helpful, but in this case the faster way to reach machine precision for solving the Poisson problem is removing all the hidden layers in the `FeedForward`, keeping only the $\alpha$ and $\beta$ parameters of the extra feature. + +# In[9]: + + +# make model + solver + trainer +model_learn = FeedForwardWithExtraFeatures( + layers=[], + func=Softplus, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables) + 1, + extra_features=SinSinAB(), +) +pinn_learn = PINN( + problem, + model_learn, + optimizer=TorchOptimizer(torch.optim.Adam, lr=0.006, weight_decay=1e-8), +) +trainer_learn = Trainer( + solver=pinn_learn, # setting the solver, i.e. PINN + max_epochs=1000, # setting max epochs in training + accelerator="cpu", # we train on cpu, also other are available + enable_model_summary=False, # model summary statistics not printed + train_size=0.8, # set train size + val_size=0.0, # set validation size + test_size=0.2, # set testing size + shuffle=True, # shuffle the data +) +# train +trainer_learn.train() + + +# In such a way, the model is able to reach a very high accuracy! +# Of course, this is a toy problem for understanding the usage of extra features: similar precision could be obtained if the extra features are very similar to the true solution. The analyzed Poisson problem shows a forcing term very close to the solution, resulting in a perfect problem to address with such an approach. + +# We conclude here by showing the test error for the analysed methodologies: the standard PINN, PINN with extra features, and PINN with learnable extra features. + +# In[12]: + + +# test error base pinn +print("PINN") +trainer_base.test() +# test error extra features pinn +print("PINN with extra features") +trainer_feat.test() +# test error learnable extra features pinn +print("PINN with learnable extra features") +_ = trainer_learn.test() + + +# ## What's next? +# +# Congratulations on completing the two dimensional Poisson tutorial of **PINA**! There are multiple directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the finaly accuracy +# +# 2. Propose new types of extrafeatures and see how they affect the learning +# +# 3. Exploit extrafeature training in more complex problems +# +# 4. Many more... diff --git a/tutorials/tutorial3/tutorial.ipynb b/tutorials/tutorial3/tutorial.ipynb index 5ffa55b9c..9962190d6 100644 --- a/tutorials/tutorial3/tutorial.ipynb +++ b/tutorials/tutorial3/tutorial.ipynb @@ -23,13 +23,14 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", - " \n", + " !pip install \"pina-mathlab\"\n", + "\n", "import torch\n", "import matplotlib.pyplot as plt\n", "import warnings\n", @@ -42,7 +43,7 @@ "from pina.equation import Equation, FixedValue\n", "from pina.callback import MetricTracker\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -261,7 +262,7 @@ " train_size=1.0,\n", " val_size=0.0,\n", " test_size=0.0,\n", - " callbacks=[MetricTracker(['train_loss', 'initial_loss','D_loss'])],\n", + " callbacks=[MetricTracker([\"train_loss\", \"initial_loss\", \"D_loss\"])],\n", ")\n", "trainer.train()" ] @@ -343,10 +344,10 @@ " \"True solution\": problem.solution(points),\n", " \"Absolute Difference\": torch.abs(\n", " solver(points) - problem.solution(points)\n", - " )\n", + " ),\n", " }\n", " # plot the solution\n", - " plt.suptitle(f'Solution for time {time.item()}')\n", + " plt.suptitle(f\"Solution for time {time.item()}\")\n", " for idx, (title, field) in enumerate(data.items()):\n", " plt.subplot(1, 3, idx + 1)\n", " plt.title(title)\n", diff --git a/tutorials/tutorial3/tutorial.py b/tutorials/tutorial3/tutorial.py new file mode 100644 index 000000000..82c574531 --- /dev/null +++ b/tutorials/tutorial3/tutorial.py @@ -0,0 +1,336 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Two dimensional Wave problem with hard constraint +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial3/tutorial.ipynb) +# +# In this tutorial we present how to solve the wave equation using hard constraint PINNs. For doing so we will build a costum `torch` model and pass it to the `PINN` solver. +# +# First of all, some useful imports. + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import matplotlib.pyplot as plt +import warnings + +from pina import Condition, LabelTensor, Trainer +from pina.problem import SpatialProblem, TimeDependentProblem +from pina.operator import laplacian, grad +from pina.domain import CartesianDomain +from pina.solver import PINN +from pina.equation import Equation, FixedValue +from pina.callback import MetricTracker + +warnings.filterwarnings("ignore") + + +# ## The problem definition + +# The problem is written in the following form: +# +# \begin{equation} +# \begin{cases} +# \Delta u(x,y,t) = \frac{\partial^2}{\partial t^2} u(x,y,t) \quad \text{in } D, \\\\ +# u(x, y, t=0) = \sin(\pi x)\sin(\pi y), \\\\ +# u(x, y, t) = 0 \quad \text{on } \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4, +# \end{cases} +# \end{equation} +# +# where $D$ is a squared domain $[0,1]^2$, and $\Gamma_i$, with $i=1,...,4$, are the boundaries of the square, and the velocity in the standard wave equation is fixed to one. + +# Now, the wave problem is written in PINA code as a class, inheriting from `SpatialProblem` and `TimeDependentProblem` since we deal with spatial, and time dependent variables. The equations are written as `conditions` that should be satisfied in the corresponding domains. `solution` is the exact solution which will be compared with the predicted one. + +# In[2]: + + +def wave_equation(input_, output_): + u_t = grad(output_, input_, components=["u"], d=["t"]) + u_tt = grad(u_t, input_, components=["dudt"], d=["t"]) + nabla_u = laplacian(output_, input_, components=["u"], d=["x", "y"]) + return nabla_u - u_tt + + +def initial_condition(input_, output_): + u_expected = torch.sin(torch.pi * input_.extract(["x"])) * torch.sin( + torch.pi * input_.extract(["y"]) + ) + return output_.extract(["u"]) - u_expected + + +class Wave(TimeDependentProblem, SpatialProblem): + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 1], "y": [0, 1]}) + temporal_domain = CartesianDomain({"t": [0, 1]}) + domains = { + "g1": CartesianDomain({"x": 1, "y": [0, 1], "t": [0, 1]}), + "g2": CartesianDomain({"x": 0, "y": [0, 1], "t": [0, 1]}), + "g3": CartesianDomain({"x": [0, 1], "y": 0, "t": [0, 1]}), + "g4": CartesianDomain({"x": [0, 1], "y": 1, "t": [0, 1]}), + "initial": CartesianDomain({"x": [0, 1], "y": [0, 1], "t": 0}), + "D": CartesianDomain({"x": [0, 1], "y": [0, 1], "t": [0, 1]}), + } + conditions = { + "g1": Condition(domain="g1", equation=FixedValue(0.0)), + "g2": Condition(domain="g2", equation=FixedValue(0.0)), + "g3": Condition(domain="g3", equation=FixedValue(0.0)), + "g4": Condition(domain="g4", equation=FixedValue(0.0)), + "initial": Condition( + domain="initial", equation=Equation(initial_condition) + ), + "D": Condition(domain="D", equation=Equation(wave_equation)), + } + + def solution(self, pts): + f = ( + torch.sin(torch.pi * pts.extract(["x"])) + * torch.sin(torch.pi * pts.extract(["y"])) + * torch.cos( + torch.sqrt(torch.tensor(2.0)) * torch.pi * pts.extract(["t"]) + ) + ) + return LabelTensor(f, self.output_variables) + + +# define problem +problem = Wave() + + +# ## Hard Constraint Model + +# After the problem, a **torch** model is needed to solve the PINN. Usually, many models are already implemented in **PINA**, but the user has the possibility to build his/her own model in `torch`. The hard constraint we impose is on the boundary of the spatial domain. Specifically, our solution is written as: +# +# $$ u_{\rm{pinn}} = xy(1-x)(1-y)\cdot NN(x, y, t), $$ +# +# where $NN$ is the neural net output. This neural network takes as input the coordinates (in this case $x$, $y$ and $t$) and provides the unknown field $u$. By construction, it is zero on the boundaries. The residuals of the equations are evaluated at several sampling points (which the user can manipulate using the method `discretise_domain`) and the loss minimized by the neural network is the sum of the residuals. + +# In[3]: + + +class HardMLP(torch.nn.Module): + + def __init__(self, input_dim, output_dim): + super().__init__() + + self.layers = torch.nn.Sequential( + torch.nn.Linear(input_dim, 40), + torch.nn.ReLU(), + torch.nn.Linear(40, 40), + torch.nn.ReLU(), + torch.nn.Linear(40, output_dim), + ) + + # here in the foward we implement the hard constraints + def forward(self, x): + hard = ( + x.extract(["x"]) + * (1 - x.extract(["x"])) + * x.extract(["y"]) + * (1 - x.extract(["y"])) + ) + return hard * self.layers(x) + + +# ## Train and Inference + +# In this tutorial, the neural network is trained for 1000 epochs with a learning rate of 0.001 (default in `PINN`). As always, we will log using `Tensorboard`. + +# In[4]: + + +# generate the data +problem.discretise_domain(1000, "random", domains="all") + +# define model +model = HardMLP(len(problem.input_variables), len(problem.output_variables)) + +# crete the solver +pinn = PINN(problem=problem, model=model) + +# create trainer and train +trainer = Trainer( + solver=pinn, + max_epochs=1000, + accelerator="cpu", + enable_model_summary=False, + train_size=1.0, + val_size=0.0, + test_size=0.0, + callbacks=[MetricTracker(["train_loss", "initial_loss", "D_loss"])], +) +trainer.train() + + +# Let's now plot the losses inside `MetricTracker` to see how they vary during training. + +# In[5]: + + +trainer_metrics = trainer.callbacks[0].metrics +for metric, loss in trainer_metrics.items(): + plt.plot(range(len(loss)), loss, label=metric) +# plotting +plt.xlabel("epoch") +plt.ylabel("loss") +plt.yscale("log") +plt.legend() + + +# Notice that the loss on the boundaries of the spatial domain is exactly zero, as expected! After the training is completed one can now plot some results using the `matplotlib`. We plot the predicted output on the left side, the true solution at the center and the difference on the right side using the `plot_solution` function. + +# In[6]: + + +@torch.no_grad() +def plot_solution(solver, time): + # get the problem + problem = solver.problem + # get spatial points + spatial_samples = problem.spatial_domain.sample(30, "grid") + # get temporal value + time = LabelTensor(torch.tensor([[time]]), "t") + # cross data + points = spatial_samples.append(time, mode="cross") + # compute pinn solution, true solution and absolute difference + data = { + "PINN solution": solver(points), + "True solution": problem.solution(points), + "Absolute Difference": torch.abs( + solver(points) - problem.solution(points) + ), + } + # plot the solution + plt.suptitle(f"Solution for time {time.item()}") + for idx, (title, field) in enumerate(data.items()): + plt.subplot(1, 3, idx + 1) + plt.title(title) + plt.tricontourf( # convert to torch tensor + flatten + points.extract("x").tensor.flatten(), + points.extract("y").tensor.flatten(), + field.tensor.flatten(), + ) + plt.colorbar(), plt.tight_layout() + + +# Let's take a look at the results at different times, for example `0.0`, `0.5` and `1.0`: + +# In[7]: + + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=0) + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=0.5) + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=1) + + +# The results are not so great, and we can clearly see that as time progresses the solution gets worse.... Can we do better? +# +# A valid option is to impose the initial condition as hard constraint as well. Specifically, our solution is written as: +# +# $$ u_{\rm{pinn}} = xy(1-x)(1-y)\cdot NN(x, y, t)\cdot t + \cos(\sqrt{2}\pi t)\sin(\pi x)\sin(\pi y), $$ +# +# Let us build the network first + +# In[8]: + + +class HardMLPtime(torch.nn.Module): + + def __init__(self, input_dim, output_dim): + super().__init__() + + self.layers = torch.nn.Sequential( + torch.nn.Linear(input_dim, 40), + torch.nn.ReLU(), + torch.nn.Linear(40, 40), + torch.nn.ReLU(), + torch.nn.Linear(40, output_dim), + ) + + # here in the foward we implement the hard constraints + def forward(self, x): + hard_space = ( + x.extract(["x"]) + * (1 - x.extract(["x"])) + * x.extract(["y"]) + * (1 - x.extract(["y"])) + ) + hard_t = ( + torch.sin(torch.pi * x.extract(["x"])) + * torch.sin(torch.pi * x.extract(["y"])) + * torch.cos( + torch.sqrt(torch.tensor(2.0)) * torch.pi * x.extract(["t"]) + ) + ) + return hard_space * self.layers(x) * x.extract(["t"]) + hard_t + + +# Now let's train with the same configuration as the previous test + +# In[9]: + + +# define model +model = HardMLPtime(len(problem.input_variables), len(problem.output_variables)) + +# crete the solver +pinn = PINN(problem=problem, model=model) + +# create trainer and train +trainer = Trainer( + solver=pinn, + max_epochs=1000, + accelerator="cpu", + enable_model_summary=False, + train_size=1.0, + val_size=0.0, + test_size=0.0, + callbacks=[MetricTracker(["train_loss", "initial_loss", "D_loss"])], +) +trainer.train() + + +# We can clearly see that the loss is way lower now. Let's plot the results + +# In[10]: + + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=0) + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=0.5) + +plt.figure(figsize=(12, 6)) +plot_solution(solver=pinn, time=1) + + +# We can see now that the results are way better! This is due to the fact that previously the network was not learning correctly the initial conditon, leading to a poor solution when time evolved. By imposing the initial condition the network is able to correctly solve the problem. + +# ## What's next? +# +# Congratulations on completing the two dimensional Wave tutorial of **PINA**! There are multiple directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the finaly accuracy +# +# 2. Propose new types of hard constraints in time, e.g. $$ u_{\rm{pinn}} = xy(1-x)(1-y)\cdot NN(x, y, t)(1-\exp(-t)) + \cos(\sqrt{2}\pi t)sin(\pi x)\sin(\pi y), $$ +# +# 3. Exploit extrafeature training for model 1 and 2 +# +# 4. Many more... diff --git a/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte b/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte new file mode 100644 index 000000000..1170b2cae Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte differ diff --git a/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte.gz b/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte.gz new file mode 100644 index 000000000..5ace8ea93 Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/t10k-images-idx3-ubyte.gz differ diff --git a/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte b/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte new file mode 100644 index 000000000..d1c3a9706 Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte differ diff --git a/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte.gz b/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte.gz new file mode 100644 index 000000000..a7e141541 Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/t10k-labels-idx1-ubyte.gz differ diff --git a/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte b/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte new file mode 100644 index 000000000..bbce27659 Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte differ diff --git a/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte.gz b/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte.gz new file mode 100644 index 000000000..b50e4b6bc Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/train-images-idx3-ubyte.gz differ diff --git a/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte b/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte new file mode 100644 index 000000000..d6b4c5db3 Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte differ diff --git a/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte.gz b/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte.gz new file mode 100644 index 000000000..707a576bb Binary files /dev/null and b/tutorials/tutorial4/data/MNIST/raw/train-labels-idx1-ubyte.gz differ diff --git a/tutorials/tutorial4/tutorial.ipynb b/tutorials/tutorial4/tutorial.ipynb index 9f3353801..f1df1b224 100644 --- a/tutorials/tutorial4/tutorial.ipynb +++ b/tutorials/tutorial4/tutorial.ipynb @@ -35,26 +35,27 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", - "import torch \n", - "import matplotlib.pyplot as plt \n", - "import torchvision # for MNIST dataset\n", + "import torch\n", + "import matplotlib.pyplot as plt\n", + "import torchvision # for MNIST dataset\n", "import warnings\n", "\n", "from pina import Trainer\n", "from pina.problem.zoo import SupervisedProblem\n", "from pina.solver import SupervisedSolver\n", "from pina.trainer import Trainer\n", - "from pina.model.block import ContinuousConvBlock \n", - "from pina.model import FeedForward # for building AE and MNIST classification\n", + "from pina.model.block import ContinuousConvBlock\n", + "from pina.model import FeedForward # for building AE and MNIST classification\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -517,7 +518,7 @@ "source": [ "# setting the problem\n", "problem = SupervisedProblem(\n", - " input_=train_data.train_data.unsqueeze(1), # adding channel dimension\n", + " input_=train_data.train_data.unsqueeze(1), # adding channel dimension\n", " output_=train_data.train_labels,\n", ")\n", "\n", @@ -568,7 +569,7 @@ "source": [ "correct = 0\n", "total = 0\n", - "trainer.data_module.setup('test')\n", + "trainer.data_module.setup(\"test\")\n", "with torch.no_grad():\n", " for data in trainer.data_module.test_dataloader():\n", " test_data = data[\"data\"]\n", @@ -580,9 +581,7 @@ " total += labels.size(0)\n", " correct += (predicted == labels).sum().item()\n", "\n", - "print(\n", - " f\"Accuracy of the network on the test images: {(correct / total):.3%}\"\n", - ")" + "print(f\"Accuracy of the network on the test images: {(correct / total):.3%}\")" ] }, { diff --git a/tutorials/tutorial4/tutorial.py b/tutorials/tutorial4/tutorial.py new file mode 100644 index 000000000..d4db53c89 --- /dev/null +++ b/tutorials/tutorial4/tutorial.py @@ -0,0 +1,676 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Unstructured convolutional autoencoder via continuous convolution +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial4/tutorial.ipynb) + +# In this tutorial, we will show how to use the Continuous Convolutional Filter, and how to build common Deep Learning architectures with it. The implementation of the filter follows the original work [*A Continuous Convolutional Trainable Filter for Modelling Unstructured Data*](https://arxiv.org/abs/2210.13416). + +# First of all we import the modules needed for the tutorial: + +# In[ ]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import matplotlib.pyplot as plt +import torchvision # for MNIST dataset +import warnings + +from pina import Trainer +from pina.problem.zoo import SupervisedProblem +from pina.solver import SupervisedSolver +from pina.trainer import Trainer +from pina.model.block import ContinuousConvBlock +from pina.model import FeedForward # for building AE and MNIST classification + +warnings.filterwarnings("ignore") + + +# The tutorial is structured as follow: +# * [Continuous filter background](#continuous-filter-background): understand how the convolutional filter works and how to use it. +# * [Building a MNIST Classifier](#building-a-mnist-classifier): show how to build a simple classifier using the MNIST dataset and how to combine a continuous convolutional layer with a feedforward neural network. +# * [Building a Continuous Convolutional Autoencoder](#building-a-continuous-convolutional-autoencoder): show how to use the continuous filter to work with unstructured data for autoencoding and up-sampling. + +# ## Continuous filter background + +# As reported by the authors in the original paper: in contrast to discrete convolution, continuous convolution is mathematically defined as: +# +# $$ +# \mathcal{I}_{\rm{out}}(\mathbf{x}) = \int_{\mathcal{X}} \mathcal{I}(\mathbf{x} + \mathbf{\tau}) \cdot \mathcal{K}(\mathbf{\tau}) d\mathbf{\tau}, +# $$ +# where $\mathcal{K} : \mathcal{X} \rightarrow \mathbb{R}$ is the *continuous filter* function, and $\mathcal{I} : \Omega \subset \mathbb{R}^N \rightarrow \mathbb{R}$ is the input function. The continuous filter function is approximated using a FeedForward Neural Network, thus trainable during the training phase. The way in which the integral is approximated can be different, currently on **PINA** we approximate it using a simple sum, as suggested by the authors. Thus, given $\{\mathbf{x}_i\}_{i=1}^{n}$ points in $\mathbb{R}^N$ of the input function mapped on the $\mathcal{X}$ filter domain, we approximate the above equation as: +# $$ +# \mathcal{I}_{\rm{out}}(\mathbf{\tilde{x}}_i) = \sum_{{\mathbf{x}_i}\in\mathcal{X}} \mathcal{I}(\mathbf{x}_i + \mathbf{\tau}) \cdot \mathcal{K}(\mathbf{x}_i), +# $$ +# where $\mathbf{\tau} \in \mathcal{S}$, with $\mathcal{S}$ the set of available strides, corresponds to the current stride position of the filter, and $\mathbf{\tilde{x}}_i$ points are obtained by taking the centroid of the filter position mapped on the $\Omega$ domain. + +# We will now try to pratically see how to work with the filter. From the above definition we see that what is needed is: +# 1. A domain and a function defined on that domain (the input) +# 2. A stride, corresponding to the positions where the filter needs to be $\rightarrow$ `stride` variable in `ContinuousConv` +# 3. The filter rectangular domain $\rightarrow$ `filter_dim` variable in `ContinuousConv` + +# ### Input function +# +# The input function for the continuous filter is defined as a tensor of shape: $$[B \times N_{in} \times N \times D]$$ where $B$ is the batch_size, $N_{in}$ is the number of input fields, $N$ the number of points in the mesh, $D$ the dimension of the problem. In particular: +# * $D$ is the number of spatial variables + 1. The last column must contain the field value. For example for 2D problems $D=3$ and the tensor will be something like `[first coordinate, second coordinate, field value]` +# * $N_{in}$ represents the number of vectorial function presented. For example a vectorial function $f = [f_1, f_2]$ will have $N_{in}=2$ +# +# Let's see an example to clear the ideas. We will be verbose to explain in details the input form. We wish to create the function: +# $$ +# f(x, y) = [\sin(\pi x) \sin(\pi y), -\sin(\pi x) \sin(\pi y)] \quad (x,y)\in[0,1]\times[0,1] +# $$ +# +# using a batch size equal to 1. + +# In[2]: + + +# batch size fixed to 1 +batch_size = 1 + +# points in the mesh fixed to 200 +N = 200 + +# vectorial 2 dimensional function, number_input_fields=2 +number_input_fields = 2 + +# 2 dimensional spatial variables, D = 2 + 1 = 3 +D = 3 + +# create the function f domain as random 2d points in [0, 1] +domain = torch.rand(size=(batch_size, number_input_fields, N, D - 1)) +print(f"Domain has shape: {domain.shape}") + +# create the functions +pi = torch.acos(torch.tensor([-1.0])) # pi value +f1 = torch.sin(pi * domain[:, 0, :, 0]) * torch.sin(pi * domain[:, 0, :, 1]) +f2 = -torch.sin(pi * domain[:, 1, :, 0]) * torch.sin(pi * domain[:, 1, :, 1]) + +# stacking the input domain and field values +data = torch.empty(size=(batch_size, number_input_fields, N, D)) +data[..., :-1] = domain # copy the domain +data[:, 0, :, -1] = f1 # copy first field value +data[:, 1, :, -1] = f1 # copy second field value +print(f"Filter input data has shape: {data.shape}") + + +# ### Stride +# +# The stride is passed as a dictionary `stride` which tells the filter where to go. Here is an example for the $[0,1]\times[0,5]$ domain: +# +# ```python +# # stride definition +# stride = {"domain": [1, 5], +# "start": [0, 0], +# "jump": [0.1, 0.3], +# "direction": [1, 1], +# } +# ``` +# This tells the filter: +# 1. `domain`: square domain (the only implemented) $[0,1]\times[0,5]$. The minimum value is always zero, while the maximum is specified by the user +# 2. `start`: start position of the filter, coordinate $(0, 0)$ +# 3. `jump`: the jumps of the centroid of the filter to the next position $(0.1, 0.3)$ +# 4. `direction`: the directions of the jump, with `1 = right`, `0 = no jump`, `-1 = left` with respect to the current position +# +# **Note** +# +# We are planning to release the possibility to directly pass a list of possible strides! + +# ### Filter definition +# +# Having defined all the previous blocks, we are now able to construct the continuous filter. +# +# Suppose we would like to get an output with only one field, and let us fix the filter dimension to be $[0.1, 0.1]$. + +# In[3]: + + +# filter dim +filter_dim = [0.1, 0.1] + +# stride +stride = { + "domain": [1, 1], + "start": [0, 0], + "jump": [0.08, 0.08], + "direction": [1, 1], +} + +# creating the filter +cConv = ContinuousConvBlock( + input_numb_field=number_input_fields, + output_numb_field=1, + filter_dim=filter_dim, + stride=stride, +) + + +# That's it! In just one line of code we have created the continuous convolutional filter. By default the `pina.model.FeedForward` neural network is intitialised, more on the [documentation](https://mathlab.github.io/PINA/_rst/fnn.html). In case the mesh doesn't change during training we can set the `optimize` flag equals to `True`, to exploit optimizations for finding the points to convolve. + +# In[4]: + + +# creating the filter + optimization +cConv = ContinuousConvBlock( + input_numb_field=number_input_fields, + output_numb_field=1, + filter_dim=filter_dim, + stride=stride, + optimize=True, +) + + +# Let's try to do a forward pass: + +# In[5]: + + +print(f"Filter input data has shape: {data.shape}") + +# input to the filter +output = cConv(data) + +print(f"Filter output data has shape: {output.shape}") + + +# If we don't want to use the default `FeedForward` neural network, we can pass a specified torch model in the `model` keyword as follow: +# + +# In[6]: + + +class SimpleKernel(torch.nn.Module): + def __init__(self) -> None: + super().__init__() + self.model = torch.nn.Sequential( + torch.nn.Linear(2, 20), + torch.nn.ReLU(), + torch.nn.Linear(20, 20), + torch.nn.ReLU(), + torch.nn.Linear(20, 1), + ) + + def forward(self, x): + return self.model(x) + + +cConv = ContinuousConvBlock( + input_numb_field=number_input_fields, + output_numb_field=1, + filter_dim=filter_dim, + stride=stride, + optimize=True, + model=SimpleKernel, +) + + +# Notice that we pass the class and not an already built object! + +# ## Building a MNIST Classifier +# +# Let's see how we can build a MNIST classifier using a continuous convolutional filter. We will use the MNIST dataset from PyTorch. In order to keep small training times we use only 6000 samples for training and 1000 samples for testing. + +# In[7]: + + +from torch.utils.data import DataLoader, SubsetRandomSampler + +numb_training = 6000 # get just 6000 images for training +numb_testing = 1000 # get just 1000 images for training +seed = 111 # for reproducibility +batch_size = 8 # setting batch size + +# setting the seed +torch.manual_seed(seed) + +# downloading the dataset +train_data = torchvision.datasets.MNIST( + "./data/", + download=True, + train=False, + transform=torchvision.transforms.Compose( + [ + torchvision.transforms.ToTensor(), + torchvision.transforms.Normalize((0.1307,), (0.3081,)), + ] + ), +) + + +# Let's now build a simple classifier. The MNIST dataset is composed by vectors of shape `[batch, 1, 28, 28]`, but we can image them as one field functions where the pixels $ij$ are the coordinate $x=i, y=j$ in a $[0, 27]\times[0,27]$ domain, and the pixels values are the field values. We just need a function to transform the regular tensor in a tensor compatible for the continuous filter: + +# In[8]: + + +def transform_input(x): + batch_size = x.shape[0] + dim_grid = tuple(x.shape[:-3:-1]) + + # creating the n dimensional mesh grid for a single channel image + values_mesh = [torch.arange(0, dim).float() for dim in dim_grid] + mesh = torch.meshgrid(values_mesh) + coordinates_mesh = [m.reshape(-1, 1).to(x.device) for m in mesh] + coordinates = ( + torch.cat(coordinates_mesh, dim=1) + .unsqueeze(0) + .repeat((batch_size, 1, 1)) + .unsqueeze(1) + ) + + return torch.cat((coordinates, x.flatten(2).unsqueeze(-1)), dim=-1) + + +# We can now build a simple classifier! We will use just one convolutional filter followed by a feedforward neural network + +# In[9]: + + +# setting the seed +torch.manual_seed(seed) + + +class ContinuousClassifier(torch.nn.Module): + def __init__(self): + super().__init__() + + # number of classes for classification + numb_class = 10 + + # convolutional block + self.convolution = ContinuousConvBlock( + input_numb_field=1, + output_numb_field=4, + stride={ + "domain": [27, 27], + "start": [0, 0], + "jumps": [4, 4], + "direction": [1, 1.0], + }, + filter_dim=[4, 4], + optimize=True, + ) + # feedforward net + self.nn = FeedForward( + input_dimensions=196, + output_dimensions=numb_class, + layers=[120, 64], + func=torch.nn.ReLU, + ) + + def forward(self, x): + # transform input + convolution + x = transform_input(x) + x = self.convolution(x) + # feed forward classification + return self.nn(x[..., -1].flatten(1)) + + +# We now aim to solve the classification problem. For this we will use the `SupervisedSolver` and the `SupervisedProblem`. The input of the supervised problems are the images, while the output the corresponding class. + +# In[10]: + + +# setting the problem +problem = SupervisedProblem( + input_=train_data.train_data.unsqueeze(1), # adding channel dimension + output_=train_data.train_labels, +) + +# setting the solver +solver = SupervisedSolver( + problem=problem, + model=ContinuousClassifier(), + loss=torch.nn.CrossEntropyLoss(), + use_lt=False, +) + +# setting the trainer +trainer = Trainer( + solver=solver, + max_epochs=1, + accelerator="cpu", + enable_model_summary=False, + train_size=0.7, + val_size=0.1, + test_size=0.2, + batch_size=64, +) +trainer.train() + + +# Let's see the performance on the test set! + +# In[11]: + + +correct = 0 +total = 0 +trainer.data_module.setup("test") +with torch.no_grad(): + for data in trainer.data_module.test_dataloader(): + test_data = data["data"] + images, labels = test_data["input"], test_data["target"] + # calculate outputs by running images through the network + outputs = solver(images) + # the class with the highest energy is what we choose as prediction + _, predicted = torch.max(outputs.data, 1) + total += labels.size(0) + correct += (predicted == labels).sum().item() + +print(f"Accuracy of the network on the test images: {(correct / total):.3%}") + + +# As we can see we have very good performance for having trained only for 1 epoch! Nevertheless, we are still using structured data... Let's see how we can build an autoencoder for unstructured data now. + +# ## Building a Continuous Convolutional Autoencoder +# +# Just as toy problem, we will now build an autoencoder for the following function $f(x,y)=\sin(\pi x)\sin(\pi y)$ on the unit circle domain centered in $(0.5, 0.5)$. We will also see the ability to up-sample (once trained) the results without retraining. Let's first create the input and visualize it, we will use firstly a mesh of $100$ points. + +# In[12]: + + +# create inputs +def circle_grid(N=100): + """Generate points withing a unit 2D circle centered in (0.5, 0.5) + + :param N: number of points + :type N: float + :return: [x, y] array of points + :rtype: torch.tensor + """ + + PI = torch.acos(torch.zeros(1)).item() * 2 + R = 0.5 + centerX = 0.5 + centerY = 0.5 + + r = R * torch.sqrt(torch.rand(N)) + theta = torch.rand(N) * 2 * PI + + x = centerX + r * torch.cos(theta) + y = centerY + r * torch.sin(theta) + + return torch.stack([x, y]).T + + +# create the grid +grid = circle_grid(500) + +# create input +input_data = torch.empty(size=(1, 1, grid.shape[0], 3)) +input_data[0, 0, :, :-1] = grid +input_data[0, 0, :, -1] = torch.sin(pi * grid[:, 0]) * torch.sin( + pi * grid[:, 1] +) + +# visualize data +plt.title("Training sample with 500 points") +plt.scatter(grid[:, 0], grid[:, 1], c=input_data[0, 0, :, -1]) +plt.colorbar() +plt.show() + + +# Let's now build a simple autoencoder using the continuous convolutional filter. The data is clearly unstructured and a simple convolutional filter might not work without projecting or interpolating first. Let's first build and `Encoder` and `Decoder` class, and then a `Autoencoder` class that contains both. + +# In[13]: + + +class Encoder(torch.nn.Module): + def __init__(self, hidden_dimension): + super().__init__() + + # convolutional block + self.convolution = ContinuousConvBlock( + input_numb_field=1, + output_numb_field=2, + stride={ + "domain": [1, 1], + "start": [0, 0], + "jumps": [0.05, 0.05], + "direction": [1, 1.0], + }, + filter_dim=[0.15, 0.15], + optimize=True, + ) + # feedforward net + self.nn = FeedForward( + input_dimensions=400, + output_dimensions=hidden_dimension, + layers=[240, 120], + ) + + def forward(self, x): + # convolution + x = self.convolution(x) + # feed forward pass + return self.nn(x[..., -1]) + + +class Decoder(torch.nn.Module): + def __init__(self, hidden_dimension): + super().__init__() + + # convolutional block + self.convolution = ContinuousConvBlock( + input_numb_field=2, + output_numb_field=1, + stride={ + "domain": [1, 1], + "start": [0, 0], + "jumps": [0.05, 0.05], + "direction": [1, 1.0], + }, + filter_dim=[0.15, 0.15], + optimize=True, + ) + # feedforward net + self.nn = FeedForward( + input_dimensions=hidden_dimension, + output_dimensions=400, + layers=[120, 240], + ) + + def forward(self, weights, grid): + # feed forward pass + x = self.nn(weights) + # transpose convolution + return torch.sigmoid(self.convolution.transpose(x, grid)) + + +# Very good! Notice that in the `Decoder` class in the `forward` pass we have used the `.transpose()` method of the `ContinuousConvolution` class. This method accepts the `weights` for upsampling and the `grid` on where to upsample. Let's now build the autoencoder! We set the hidden dimension in the `hidden_dimension` variable. We apply the sigmoid on the output since the field value is between $[0, 1]$. + +# In[14]: + + +class Autoencoder(torch.nn.Module): + def __init__(self, hidden_dimension=10): + super().__init__() + + self.encoder = Encoder(hidden_dimension) + self.decoder = Decoder(hidden_dimension) + + def forward(self, x): + # saving grid for later upsampling + grid = x.clone().detach() + # encoder + weights = self.encoder(x) + # decoder + out = self.decoder(weights, grid) + return out + + +# Let's now train the autoencoder, minimizing the mean square error loss and optimizing using Adam. We use the `SupervisedSolver` as solver, and the problem is a simple problem created by inheriting from `AbstractProblem`. It takes approximately two minutes to train on CPU. + +# In[15]: + + +# define the problem +problem = SupervisedProblem(input_data, input_data) + + +# define the solver +solver = SupervisedSolver( + problem=problem, + model=Autoencoder(), + loss=torch.nn.MSELoss(), + use_lt=False, +) + +# train +trainer = Trainer( + solver, + max_epochs=100, + accelerator="cpu", + enable_model_summary=False, # we train on CPU and avoid model summary at beginning of training (optional) + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# Let's visualize the two solutions side by side! + +# In[16]: + + +solver.eval() + +# get output and detach from computational graph for plotting +output = solver(input_data).detach() + +# visualize data +fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3)) +pic1 = axes[0].scatter(grid[:, 0], grid[:, 1], c=input_data[0, 0, :, -1]) +axes[0].set_title("Real") +fig.colorbar(pic1) +plt.subplot(1, 2, 2) +pic2 = axes[1].scatter(grid[:, 0], grid[:, 1], c=output[0, 0, :, -1]) +axes[1].set_title("Autoencoder") +fig.colorbar(pic2) +plt.tight_layout() +plt.show() + + +# As we can see, the two solutions are really similar! We can compute the $l_2$ error quite easily as well: + +# In[17]: + + +def l2_error(input_, target): + return torch.linalg.norm(input_ - target, ord=2) / torch.linalg.norm( + input_, ord=2 + ) + + +print(f"l2 error: {l2_error(input_data[0, 0, :, -1], output[0, 0, :, -1]):.2%}") + + +# More or less $4\%$ in $l_2$ error, which is really low considering the fact that we use just **one** convolutional layer and a simple feedforward to decrease the dimension. Let's see now some peculiarity of the filter. + +# ### Filter for upsampling +# +# Suppose we have already the hidden representation and we want to upsample on a differen grid with more points. Let's see how to do it: + +# In[18]: + + +# setting the seed +torch.manual_seed(seed) + +grid2 = circle_grid(1500) # triple number of points +input_data2 = torch.zeros(size=(1, 1, grid2.shape[0], 3)) +input_data2[0, 0, :, :-1] = grid2 +input_data2[0, 0, :, -1] = torch.sin(pi * grid2[:, 0]) * torch.sin( + pi * grid2[:, 1] +) + +# get the hidden representation from original input +latent = solver.model.encoder(input_data) + +# upsample on the second input_data2 +output = solver.model.decoder(latent, input_data2).detach() + +# show the picture +fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3)) +pic1 = axes[0].scatter(grid2[:, 0], grid2[:, 1], c=input_data2[0, 0, :, -1]) +axes[0].set_title("Real") +fig.colorbar(pic1) +plt.subplot(1, 2, 2) +pic2 = axes[1].scatter(grid2[:, 0], grid2[:, 1], c=output[0, 0, :, -1]) +axes[1].set_title("Up-sampling") +fig.colorbar(pic2) +plt.tight_layout() +plt.show() + + +# As we can see we have a very good approximation of the original function, even thought some noise is present. Let's calculate the error now: + +# In[19]: + + +print( + f"l2 error: {l2_error(input_data2[0, 0, :, -1], output[0, 0, :, -1]):.2%}" +) + + +# ### Autoencoding at different resolutions +# In the previous example we already had the hidden representation (of the original input) and we used it to upsample. Sometimes however we could have a finer mesh solution and we would simply want to encode it. This can be done without retraining! This procedure can be useful in case we have many points in the mesh and just a smaller part of them are needed for training. Let's see the results of this: + +# In[20]: + + +# setting the seed +torch.manual_seed(seed) + +grid2 = circle_grid(3500) # very fine mesh +input_data2 = torch.zeros(size=(1, 1, grid2.shape[0], 3)) +input_data2[0, 0, :, :-1] = grid2 +input_data2[0, 0, :, -1] = torch.sin(pi * grid2[:, 0]) * torch.sin( + pi * grid2[:, 1] +) + +# get the hidden representation from finer mesh input +latent = solver.model.encoder(input_data2) + +# upsample on the second input_data2 +output = solver.model.decoder(latent, input_data2).detach() + +# show the picture +fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3)) +pic1 = axes[0].scatter(grid2[:, 0], grid2[:, 1], c=input_data2[0, 0, :, -1]) +axes[0].set_title("Real") +fig.colorbar(pic1) +plt.subplot(1, 2, 2) +pic2 = axes[1].scatter(grid2[:, 0], grid2[:, 1], c=output[0, 0, :, -1]) +axes[1].set_title("Autoencoder not re-trained") +fig.colorbar(pic2) +plt.tight_layout() +plt.show() + +# calculate l2 error +print( + f"l2 error: {l2_error(input_data2[0, 0, :, -1], output[0, 0, :, -1]):.2%}" +) + + +# ## What's next? +# +# We have shown the basic usage of a convolutional filter. There are additional extensions possible: +# +# 1. Train using Physics Informed strategies +# +# 2. Use the filter to build an unstructured convolutional autoencoder for reduced order modelling +# +# 3. Many more... diff --git a/tutorials/tutorial5/tutorial.ipynb b/tutorials/tutorial5/tutorial.ipynb index 78d59afed..b02afd90c 100644 --- a/tutorials/tutorial5/tutorial.ipynb +++ b/tutorials/tutorial5/tutorial.ipynb @@ -33,15 +33,16 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", - " !pip install scipy\n", - " # get the data\n", - " !wget https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial5/Data_Darcy.mat\n", + " !pip install \"pina-mathlab\"\n", + " !pip install scipy\n", + " # get the data\n", + " !wget https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial5/Data_Darcy.mat\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", @@ -54,7 +55,7 @@ "from pina.solver import SupervisedSolver\n", "from pina.problem.zoo import SupervisedProblem\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -129,10 +130,10 @@ ], "source": [ "plt.subplot(1, 2, 1)\n", - "plt.title('permeability')\n", + "plt.title(\"permeability\")\n", "plt.imshow(k_train[0])\n", "plt.subplot(1, 2, 2)\n", - "plt.title('field solution')\n", + "plt.title(\"field solution\")\n", "plt.imshow(u_train[0])\n", "plt.show()" ] @@ -278,12 +279,10 @@ " )\n", " * 100\n", ")\n", - "print(f'Final error training {err:.2f}%')\n", + "print(f\"Final error training {err:.2f}%\")\n", "\n", "err = (\n", - " float(\n", - " metric_err(u_test.unsqueeze(-1), model(k_test.unsqueeze(-1))).mean()\n", - " )\n", + " float(metric_err(u_test.unsqueeze(-1), model(k_test.unsqueeze(-1))).mean())\n", " * 100\n", ")\n", "print(f\"Final error testing {err:.2f}%\")" diff --git a/tutorials/tutorial5/tutorial.py b/tutorials/tutorial5/tutorial.py new file mode 100644 index 000000000..6783fd557 --- /dev/null +++ b/tutorials/tutorial5/tutorial.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Two dimensional Darcy flow using the Fourier Neural Operator +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial5/tutorial.ipynb) +# + +# In this tutorial we are going to solve the Darcy flow problem in two dimensions, presented in [*Fourier Neural Operator for +# Parametric Partial Differential Equation*](https://openreview.net/pdf?id=c8P9NQVtmnO). First of all we import the modules needed for the tutorial. Importing `scipy` is needed for input-output operations. + +# In[1]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + get_ipython().system("pip install scipy") + # get the data + get_ipython().system( + "wget https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial5/Data_Darcy.mat" + ) + +import torch +import matplotlib.pyplot as plt +import warnings + +# !pip install scipy # install scipy +from scipy import io +from pina.model import FNO, FeedForward # let's import some models +from pina import Condition, Trainer +from pina.solver import SupervisedSolver +from pina.problem.zoo import SupervisedProblem + +warnings.filterwarnings("ignore") + + +# ## Data Generation +# +# We will focus on solving a specific PDE, the **Darcy Flow** equation. The Darcy PDE is a second-order elliptic PDE with the following form: +# +# $$ +# -\nabla\cdot(k(x, y)\nabla u(x, y)) = f(x) \quad (x, y) \in D. +# $$ +# +# Specifically, $u$ is the flow pressure, $k$ is the permeability field and $f$ is the forcing function. The Darcy flow can parameterize a variety of systems including flow through porous media, elastic materials and heat conduction. Here you will define the domain as a 2D unit square Dirichlet boundary conditions. The dataset is taken from the authors original reference. +# + +# In[2]: + + +# download the dataset +data = io.loadmat("Data_Darcy.mat") + +# extract data (we use only 100 data for train) +k_train = torch.tensor(data["k_train"], dtype=torch.float) +u_train = torch.tensor(data["u_train"], dtype=torch.float) +k_test = torch.tensor(data["k_test"], dtype=torch.float) +u_test = torch.tensor(data["u_test"], dtype=torch.float) +x = torch.tensor(data["x"], dtype=torch.float)[0] +y = torch.tensor(data["y"], dtype=torch.float)[0] + + +# Let's visualize some data + +# In[3]: + + +plt.subplot(1, 2, 1) +plt.title("permeability") +plt.imshow(k_train[0]) +plt.subplot(1, 2, 2) +plt.title("field solution") +plt.imshow(u_train[0]) +plt.show() + + +# We now create the Neural Operators problem class. Learning Neural Operators is similar as learning in a supervised manner, therefore we will use `SupervisedProblem`. + +# In[4]: + + +# make problem +problem = SupervisedProblem( + input_=k_train.unsqueeze(-1), output_=u_train.unsqueeze(-1) +) + + +# ## Solving the problem with a FeedForward Neural Network +# +# We will first solve the problem using a Feedforward neural network. We will use the `SupervisedSolver` for solving the problem, since we are training using supervised learning. + +# In[5]: + + +# make model +model = FeedForward(input_dimensions=1, output_dimensions=1) + + +# make solver +solver = SupervisedSolver(problem=problem, model=model, use_lt=False) + +# make the trainer and train +trainer = Trainer( + solver=solver, + max_epochs=10, + accelerator="cpu", + enable_model_summary=False, + batch_size=10, + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# The final loss is pretty high... We can calculate the error by importing `LpLoss`. + +# In[6]: + + +from pina.loss import LpLoss + +# make the metric +metric_err = LpLoss(relative=False) + +model = solver.model +err = ( + float( + metric_err(u_train.unsqueeze(-1), model(k_train.unsqueeze(-1))).mean() + ) + * 100 +) +print(f"Final error training {err:.2f}%") + +err = ( + float(metric_err(u_test.unsqueeze(-1), model(k_test.unsqueeze(-1))).mean()) + * 100 +) +print(f"Final error testing {err:.2f}%") + + +# ## Solving the problem with a Fourier Neural Operator (FNO) +# +# We will now move to solve the problem using a FNO. Since we are learning operator this approach is better suited, as we shall see. + +# In[7]: + + +# make model +lifting_net = torch.nn.Linear(1, 24) +projecting_net = torch.nn.Linear(24, 1) +model = FNO( + lifting_net=lifting_net, + projecting_net=projecting_net, + n_modes=8, + dimensions=2, + inner_size=24, + padding=8, +) + + +# make solver +solver = SupervisedSolver(problem=problem, model=model, use_lt=False) + +# make the trainer and train +trainer = Trainer( + solver=solver, + max_epochs=10, + accelerator="cpu", + enable_model_summary=False, + batch_size=10, + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# We can clearly see that the final loss is lower. Let's see in testing.. Notice that the number of parameters is way higher than a `FeedForward` network. We suggest to use GPU or TPU for a speed up in training, when many data samples are used. + +# In[8]: + + +model = solver.model +err = ( + float( + metric_err(u_train.unsqueeze(-1), model(k_train.unsqueeze(-1))).mean() + ) + * 100 +) +print(f"Final error training {err:.2f}%") + +err = ( + float(metric_err(u_test.unsqueeze(-1), model(k_test.unsqueeze(-1))).mean()) + * 100 +) +print(f"Final error testing {err:.2f}%") + + +# As we can see the loss is way lower! + +# ## What's next? +# +# We have made a very simple example on how to use the `FNO` for learning neural operator. Currently in **PINA** we implement 1D/2D/3D cases. We suggest to extend the tutorial using more complex problems and train for longer, to see the full potential of neural operators. diff --git a/tutorials/tutorial6/tutorial.ipynb b/tutorials/tutorial6/tutorial.ipynb index 563586178..522a9087f 100644 --- a/tutorials/tutorial6/tutorial.ipynb +++ b/tutorials/tutorial6/tutorial.ipynb @@ -26,21 +26,30 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import matplotlib.pyplot as plt\n", "\n", - "from pina.domain import EllipsoidDomain, Difference, CartesianDomain, Union, SimplexDomain, DomainInterface\n", + "from pina.domain import (\n", + " EllipsoidDomain,\n", + " Difference,\n", + " CartesianDomain,\n", + " Union,\n", + " SimplexDomain,\n", + " DomainInterface,\n", + ")\n", "from pina.label_tensor import LabelTensor\n", "\n", + "\n", "def plot_scatter(ax, pts, title):\n", " ax.title.set_text(title)\n", - " ax.scatter(pts.extract('x'), pts.extract('y'), color='blue', alpha=0.5)" + " ax.scatter(pts.extract(\"x\"), pts.extract(\"y\"), color=\"blue\", alpha=0.5)" ] }, { diff --git a/tutorials/tutorial6/tutorial.py b/tutorials/tutorial6/tutorial.py new file mode 100644 index 000000000..b35518434 --- /dev/null +++ b/tutorials/tutorial6/tutorial.py @@ -0,0 +1,293 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Building custom geometries with PINA `DomainInterface` class +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial6/tutorial.ipynb) +# +# In this tutorial we will show how to use geometries in PINA. Specifically, the tutorial will include how to create geometries and how to visualize them. The topics covered are: +# +# * Creating CartesianDomains and EllipsoidDomains +# * Getting the Union and Difference of Geometries +# * Sampling points in the domain (and visualize them) +# +# We import the relevant modules first. + +# In[ ]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import matplotlib.pyplot as plt + +from pina.domain import ( + EllipsoidDomain, + Difference, + CartesianDomain, + Union, + SimplexDomain, + DomainInterface, +) +from pina.label_tensor import LabelTensor + + +def plot_scatter(ax, pts, title): + ax.title.set_text(title) + ax.scatter(pts.extract("x"), pts.extract("y"), color="blue", alpha=0.5) + + +# ## Built-in Geometries + +# We will create one cartesian and two ellipsoids. For the sake of simplicity, we show here the 2-dimensional case, but the extension to 3D (and higher) cases is trivial. The geometries allow also the generation of samples belonging to the boundary. So, we will create one ellipsoid with the border and one without. + +# In[ ]: + + +cartesian = CartesianDomain({"x": [0, 2], "y": [0, 2]}) +ellipsoid_no_border = EllipsoidDomain({"x": [1, 3], "y": [1, 3]}) +ellipsoid_border = EllipsoidDomain( + {"x": [2, 4], "y": [2, 4]}, sample_surface=True +) + + +# The `{'x': [0, 2], 'y': [0, 2]}` are the bounds of the `CartesianDomain` being created. +# +# To visualize these shapes, we need to sample points on them. We will use the `sample` method of the `CartesianDomain` and `EllipsoidDomain` classes. This method takes a `n` argument which is the number of points to sample. It also takes different modes to sample, such as `'random'`. + +# In[ ]: + + +cartesian_samples = cartesian.sample(n=1000, mode="random") +ellipsoid_no_border_samples = ellipsoid_no_border.sample(n=1000, mode="random") +ellipsoid_border_samples = ellipsoid_border.sample(n=1000, mode="random") + + +# We can see the samples of each geometry to see what we are working with. + +# In[4]: + + +print(f"Cartesian Samples: {cartesian_samples}") +print(f"Ellipsoid No Border Samples: {ellipsoid_no_border_samples}") +print(f"Ellipsoid Border Samples: {ellipsoid_border_samples}") + + +# Notice how these are all `LabelTensor` objects. You can read more about these in the [documentation](https://mathlab.github.io/PINA/_rst/label_tensor.html). At a very high level, they are tensors where each element in a tensor has a label that we can access by doing `.labels`. We can also access the values of the tensor by doing `.extract(['x'])`. +# +# We are now ready to visualize the samples using matplotlib. + +# In[ ]: + + +fig, axs = plt.subplots(1, 3, figsize=(16, 4)) +pts_list = [ + cartesian_samples, + ellipsoid_no_border_samples, + ellipsoid_border_samples, +] +title_list = ["Cartesian Domain", "Ellipsoid Domain", "Ellipsoid Border Domain"] +for ax, pts, title in zip(axs, pts_list, title_list): + plot_scatter(ax, pts, title) + + +# We have now created, sampled, and visualized our first geometries! We can see that the `EllipsoidDomain` with the border has a border around it. We can also see that the `EllipsoidDomain` without the border is just the ellipse. We can also see that the `CartesianDomain` is just a square. + +# ### Simplex Domain +# +# Among the built-in shapes, we quickly show here the usage of `SimplexDomain`, which can be used for polygonal domains! + +# In[ ]: + + +import torch + +spatial_domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[1, 1]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "y"]), + ] +) + +spatial_domain2 = SimplexDomain( + [ + LabelTensor(torch.tensor([[0.0, -2.0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[-0.5, -0.5]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[-2.0, 0.0]]), labels=["x", "y"]), + ] +) + +pts = spatial_domain2.sample(100) +fig, axs = plt.subplots(1, 2, figsize=(16, 6)) +for domain, ax in zip([spatial_domain, spatial_domain2], axs): + pts = domain.sample(1000) + plot_scatter(ax, pts, "Simplex Domain") + + +# ## Boolean Operations + +# To create complex shapes we can use the boolean operations, for example to merge two default geometries. We need to simply use the `Union` class: it takes a list of geometries and returns the union of them. +# +# Let's create three unions. Firstly, it will be a union of `cartesian` and `ellipsoid_no_border`. Next, it will be a union of `ellipse_no_border` and `ellipse_border`. Lastly, it will be a union of all three geometries. + +# In[7]: + + +cart_ellipse_nb_union = Union([cartesian, ellipsoid_no_border]) +cart_ellipse_b_union = Union([cartesian, ellipsoid_border]) +three_domain_union = Union([cartesian, ellipsoid_no_border, ellipsoid_border]) + + +# We can of course sample points over the new geometries, by using the `sample` method as before. We highlight that the available sample strategy here is only *random*. + +# In[ ]: + + +c_e_nb_u_points = cart_ellipse_nb_union.sample(n=2000, mode="random") +c_e_b_u_points = cart_ellipse_b_union.sample(n=2000, mode="random") +three_domain_union_points = three_domain_union.sample(n=3000, mode="random") + + +# We can plot the samples of each of the unions to see what we are working with. + +# In[ ]: + + +fig, axs = plt.subplots(1, 3, figsize=(16, 4)) +pts_list = [c_e_nb_u_points, c_e_b_u_points, three_domain_union_points] +title_list = [ + "Cartesian with Ellipsoid No Border Union", + "Cartesian with Ellipsoid Border Union", + "Three Domain Union", +] +for ax, pts, title in zip(axs, pts_list, title_list): + plot_scatter(ax, pts, title) + + +# Now, we will find the differences of the geometries. We will find the difference of `cartesian` and `ellipsoid_no_border`. + +# In[ ]: + + +cart_ellipse_nb_difference = Difference([cartesian, ellipsoid_no_border]) +c_e_nb_d_points = cart_ellipse_nb_difference.sample(n=2000, mode="random") + +fig, ax = plt.subplots(1, 1, figsize=(8, 6)) +plot_scatter(ax, c_e_nb_d_points, "Difference") + + +# ## Create Custom DomainInterface + +# We will take a look on how to create our own geometry. The one we will try to make is a heart defined by the function $$(x^2+y^2-1)^3-x^2y^3 \le 0$$ + +# Let's start by importing what we will need to create our own geometry based on this equation. + +# In[11]: + + +import torch +from pina import LabelTensor + + +# Next, we will create the `Heart(DomainInterface)` class and initialize it. + +# In[ ]: + + +class Heart(DomainInterface): + """Implementation of the Heart Domain.""" + + def __init__(self, sample_border=False): + super().__init__() + + +# Because the `DomainInterface` class we are inheriting from requires both a `sample` method and `is_inside` method, we will create them and just add in "pass" for the moment. We also observe that the methods `sample_modes` and `variables` of the `DomainInterface` class are initialized as `abstractmethod`, so we need to redefine them both in the subclass `Heart` . + +# In[ ]: + + +class Heart(DomainInterface): + """Implementation of the Heart Domain.""" + + def __init__(self, sample_border=False): + super().__init__() + + def is_inside(self): + pass + + def sample(self): + pass + + @property + def sample_modes(self): + pass + + @property + def variables(self): + pass + + +# Now we have the skeleton for our `Heart` class. Also the `sample` method is where most of the work is done so let's fill it out. + +# In[ ]: + + +class Heart(DomainInterface): + """Implementation of the Heart Domain.""" + + def __init__(self, sample_border=False): + super().__init__() + + def is_inside(self): + pass + + def sample(self, n): + sampled_points = [] + + while len(sampled_points) < n: + x = torch.rand(1) * 3.0 - 1.5 + y = torch.rand(1) * 3.0 - 1.5 + if ((x**2 + y**2 - 1) ** 3 - (x**2) * (y**3)) <= 0: + sampled_points.append([x.item(), y.item()]) + + return LabelTensor(torch.tensor(sampled_points), labels=["x", "y"]) + + @property + def sample_modes(self): + pass + + @property + def variables(self): + pass + + +# To create the Heart geometry we simply run: + +# In[15]: + + +heart = Heart() + + +# To sample from the Heart geometry we simply run: + +# In[ ]: + + +pts_heart = heart.sample(1500) + +fig, ax = plt.subplots() +plot_scatter(ax, pts_heart, "Heart Domain") + + +# ## What's next? +# +# We have made a very simple tutorial on how to build custom geometries and use domain operation to compose base geometries. Now you can play around with different geometries and build your own! diff --git a/tutorials/tutorial7/tutorial.ipynb b/tutorials/tutorial7/tutorial.ipynb index 573f7954e..132578d11 100644 --- a/tutorials/tutorial7/tutorial.ipynb +++ b/tutorials/tutorial7/tutorial.ipynb @@ -75,17 +75,18 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", - " # get the data\n", - " !mkdir \"data\"\n", - " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pinn_solution_0.5_0.5\" -O \"data/pinn_solution_0.5_0.5\"\n", - " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pts_0.5_0.5\" -O \"data/pts_0.5_0.5\"\n", - " \n", + " !pip install \"pina-mathlab\"\n", + " # get the data\n", + " !mkdir \"data\"\n", + " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pinn_solution_0.5_0.5\" -O \"data/pinn_solution_0.5_0.5\"\n", + " !wget \"https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pts_0.5_0.5\" -O \"data/pts_0.5_0.5\"\n", + "\n", "import matplotlib.pyplot as plt\n", "import torch\n", "import warnings\n", @@ -101,7 +102,7 @@ "from lightning.pytorch import seed_everything\n", "from lightning.pytorch.callbacks import Callback\n", "\n", - "warnings.filterwarnings('ignore')\n", + "warnings.filterwarnings(\"ignore\")\n", "seed_everything(883)" ] }, @@ -152,11 +153,11 @@ } ], "source": [ - "points = data_input.extract(['x', 'y']).detach().numpy()\n", + "points = data_input.extract([\"x\", \"y\"]).detach().numpy()\n", "truth = data_output.detach().numpy()\n", "\n", "plt.scatter(points[:, 0], points[:, 1], c=truth, s=8)\n", - "plt.axis('equal')\n", + "plt.axis(\"equal\")\n", "plt.colorbar()\n", "plt.show()" ] @@ -255,8 +256,8 @@ " layers=[20, 20, 20],\n", " func=torch.nn.Softplus,\n", " output_dimensions=len(problem.output_variables),\n", - " input_dimensions=len(problem.input_variables)\n", - " )" + " input_dimensions=len(problem.input_variables),\n", + ")" ] }, { diff --git a/tutorials/tutorial7/tutorial.py b/tutorials/tutorial7/tutorial.py new file mode 100644 index 000000000..d5c000ee0 --- /dev/null +++ b/tutorials/tutorial7/tutorial.py @@ -0,0 +1,259 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Resolution of an inverse problem +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial7/tutorial.ipynb) + +# ### Introduction to the inverse problem + +# This tutorial shows how to solve an inverse Poisson problem with Physics-Informed Neural Networks. The problem definition is that of a Poisson problem with homogeneous boundary conditions and it reads: +# \begin{equation} +# \begin{cases} +# \Delta u = e^{-2(x-\mu_1)^2-2(y-\mu_2)^2} \text{ in } \Omega\, ,\\ +# u = 0 \text{ on }\partial \Omega,\\ +# u(\mu_1, \mu_2) = \text{ data} +# \end{cases} +# \end{equation} +# where $\Omega$ is a square domain $[-2, 2] \times [-2, 2]$, and $\partial \Omega=\Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$ is the union of the boundaries of the domain. +# +# This kind of problem, namely the "inverse problem", has two main goals: +# - find the solution $u$ that satisfies the Poisson equation; +# - find the unknown parameters ($\mu_1$, $\mu_2$) that better fit some given data (third equation in the system above). +# +# In order to achieve both goals we will need to define an `InverseProblem` in PINA. + +# Let's start with useful imports. + +# In[20]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + # get the data + get_ipython().system('mkdir "data"') + get_ipython().system( + 'wget "https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pinn_solution_0.5_0.5" -O "data/pinn_solution_0.5_0.5"' + ) + get_ipython().system( + 'wget "https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial7/data/pts_0.5_0.5" -O "data/pts_0.5_0.5"' + ) + +import matplotlib.pyplot as plt +import torch +import warnings + +from pina import Condition, Trainer +from pina.problem import SpatialProblem, InverseProblem +from pina.operator import laplacian +from pina.model import FeedForward +from pina.equation import Equation, FixedValue +from pina.solver import PINN +from pina.domain import CartesianDomain +from pina.optim import TorchOptimizer +from lightning.pytorch import seed_everything +from lightning.pytorch.callbacks import Callback + +warnings.filterwarnings("ignore") +seed_everything(883) + + +# Then, we import the pre-saved data, for ($\mu_1$, $\mu_2$)=($0.5$, $0.5$). These two values are the optimal parameters that we want to find through the neural network training. In particular, we import the `input` points (the spatial coordinates), and the `target` points (the corresponding $u$ values evaluated at the `input`). + +# In[21]: + + +data_output = torch.load( + "data/pinn_solution_0.5_0.5", weights_only=False +).detach() +data_input = torch.load("data/pts_0.5_0.5", weights_only=False) + + +# Moreover, let's plot also the data points and the reference solution: this is the expected output of the neural network. + +# In[22]: + + +points = data_input.extract(["x", "y"]).detach().numpy() +truth = data_output.detach().numpy() + +plt.scatter(points[:, 0], points[:, 1], c=truth, s=8) +plt.axis("equal") +plt.colorbar() +plt.show() + + +# ### Inverse problem definition in PINA + +# Then, we initialize the Poisson problem, that is inherited from the `SpatialProblem` and from the `InverseProblem` classes. We here have to define all the variables, and the domain where our unknown parameters ($\mu_1$, $\mu_2$) belong. Notice that the Laplace equation takes as inputs also the unknown variables, that will be treated as parameters that the neural network optimizes during the training process. + +# In[23]: + + +def laplace_equation(input_, output_, params_): + """ + Implementation of the laplace equation. + + :param LabelTensor input_: Input data of the problem. + :param LabelTensor output_: Output data of the problem. + :param dict params_: Parameters of the problem. + :return: The residual of the laplace equation. + :rtype: LabelTensor + """ + force_term = torch.exp( + -2 * (input_.extract(["x"]) - params_["mu1"]) ** 2 + - 2 * (input_.extract(["y"]) - params_["mu2"]) ** 2 + ) + delta_u = laplacian(output_, input_, components=["u"], d=["x", "y"]) + return delta_u - force_term + + +class Poisson(SpatialProblem, InverseProblem): + r""" + Implementation of the inverse 2-dimensional Poisson problem in the square + domain :math:`[0, 1] \times [0, 1]`, + with unknown parameter domain :math:`[-1, 1] \times [-1, 1]`. + """ + + output_variables = ["u"] + x_min, x_max = -2, 2 + y_min, y_max = -2, 2 + spatial_domain = CartesianDomain({"x": [x_min, x_max], "y": [y_min, y_max]}) + unknown_parameter_domain = CartesianDomain({"mu1": [-1, 1], "mu2": [-1, 1]}) + + domains = { + "g1": CartesianDomain({"x": [x_min, x_max], "y": y_max}), + "g2": CartesianDomain({"x": [x_min, x_max], "y": y_min}), + "g3": CartesianDomain({"x": x_max, "y": [y_min, y_max]}), + "g4": CartesianDomain({"x": x_min, "y": [y_min, y_max]}), + "D": CartesianDomain({"x": [x_min, x_max], "y": [y_min, y_max]}), + } + + conditions = { + "g1": Condition(domain="g1", equation=FixedValue(0.0)), + "g2": Condition(domain="g2", equation=FixedValue(0.0)), + "g3": Condition(domain="g3", equation=FixedValue(0.0)), + "g4": Condition(domain="g4", equation=FixedValue(0.0)), + "D": Condition(domain="D", equation=Equation(laplace_equation)), + "data": Condition(input=data_input, target=data_output), + } + + +problem = Poisson() + + +# Then, we define the neural network model we want to use. Here we used a model which imposes hard constrains on the boundary conditions, as also done in the Wave tutorial! + +# In[24]: + + +model = FeedForward( + layers=[20, 20, 20], + func=torch.nn.Softplus, + output_dimensions=len(problem.output_variables), + input_dimensions=len(problem.input_variables), +) + + +# After that, we discretize the spatial domain. + +# In[25]: + + +problem.discretise_domain(20, "grid", domains=["D"]) +problem.discretise_domain( + 1000, + "random", + domains=["g1", "g2", "g3", "g4"], +) + + +# Here, we define a simple callback for the trainer. We use this callback to save the parameters predicted by the neural network during the training. The parameters are saved every 100 epochs as `torch` tensors in a specified directory (`tmp_dir` in our case). +# The goal is to read the saved parameters after training and plot their trend across the epochs. + +# In[26]: + + +# temporary directory for saving logs of training +tmp_dir = "tmp_poisson_inverse" + + +class SaveParameters(Callback): + """ + Callback to save the parameters of the model every 100 epochs. + """ + + def on_train_epoch_end(self, trainer, __): + if trainer.current_epoch % 100 == 99: + torch.save( + trainer.solver.problem.unknown_parameters, + "{}/parameters_epoch{}".format(tmp_dir, trainer.current_epoch), + ) + + +# Then, we define the `PINN` object and train the solver using the `Trainer`. + +# In[27]: + + +max_epochs = 1500 +pinn = PINN( + problem, model, optimizer=TorchOptimizer(torch.optim.Adam, lr=0.005) +) +# define the trainer for the solver +trainer = Trainer( + solver=pinn, + accelerator="cpu", + max_epochs=max_epochs, + default_root_dir=tmp_dir, + enable_model_summary=False, + callbacks=[SaveParameters()], + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# One can now see how the parameters vary during the training by reading the saved solution and plotting them. The plot shows that the parameters stabilize to their true value before reaching the epoch $1000$! + +# In[28]: + + +epochs_saved = range(99, max_epochs, 100) +parameters = torch.empty((int(max_epochs / 100), 2)) +for i, epoch in enumerate(epochs_saved): + params_torch = torch.load( + "{}/parameters_epoch{}".format(tmp_dir, epoch), weights_only=False + ) + for e, var in enumerate(pinn.problem.unknown_variables): + parameters[i, e] = params_torch[var].data + +# Plot parameters +plt.close() +plt.plot(epochs_saved, parameters[:, 0], label="mu1", marker="o") +plt.plot(epochs_saved, parameters[:, 1], label="mu2", marker="s") +plt.ylim(-1, 1) +plt.grid() +plt.legend() +plt.xlabel("Epoch") +plt.ylabel("Parameter value") +plt.show() + + +# ## What's next? +# +# We have shown the basic usage PINNs in inverse problem modelling, further extensions include: +# +# 1. Train using different Physics Informed strategies +# +# 2. Try on more complex problems +# +# 3. Many more... diff --git a/tutorials/tutorial8/tutorial.ipynb b/tutorials/tutorial8/tutorial.ipynb index acc263d6c..b5273fcdb 100644 --- a/tutorials/tutorial8/tutorial.ipynb +++ b/tutorials/tutorial8/tutorial.ipynb @@ -38,12 +38,13 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "%matplotlib inline\n", "\n", @@ -60,7 +61,7 @@ "from pina.problem.zoo import SupervisedProblem\n", "from pina.model.block import PODBlock, RBFBlock\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -93,12 +94,13 @@ ], "source": [ "from smithers.dataset import NavierStokesDataset\n", + "\n", "dataset = NavierStokesDataset()\n", "\n", "fig, axs = plt.subplots(1, 4, figsize=(14, 3))\n", - "for ax, p, u in zip(axs, dataset.params[:4], dataset.snapshots['mag(v)'][:4]):\n", + "for ax, p, u in zip(axs, dataset.params[:4], dataset.snapshots[\"mag(v)\"][:4]):\n", " ax.tricontourf(dataset.triang, u, levels=16)\n", - " ax.set_title(f'$\\mu$ = {p[0]:.2f}')" + " ax.set_title(f\"$\\mu$ = {p[0]:.2f}\")" ] }, { @@ -118,7 +120,7 @@ "metadata": {}, "outputs": [], "source": [ - "u = torch.tensor(dataset.snapshots['mag(v)']).float()\n", + "u = torch.tensor(dataset.snapshots[\"mag(v)\"]).float()\n", "p = torch.tensor(dataset.params).float()\n", "problem = SupervisedProblem(input_=p, output_=u)" ] @@ -158,20 +160,18 @@ " \"\"\"\n", " Proper orthogonal decomposition with neural network model.\n", " \"\"\"\n", + "\n", " def __init__(self, pod_rank, layers, func):\n", - " \"\"\"\n", - " \n", - " \"\"\"\n", + " \"\"\" \"\"\"\n", " super().__init__()\n", - " \n", + "\n", " self.pod = PODBlock(pod_rank)\n", " self.nn = FeedForward(\n", " input_dimensions=1,\n", " output_dimensions=pod_rank,\n", " layers=layers,\n", - " func=func\n", + " func=func,\n", " )\n", - " \n", "\n", " def forward(self, x):\n", " \"\"\"\n", @@ -210,10 +210,11 @@ "source": [ "pod_nn = PODNN(pod_rank=20, layers=[10, 10, 10], func=torch.nn.Tanh)\n", "pod_nn_stokes = SupervisedSolver(\n", - " problem=problem, \n", - " model=pod_nn, \n", + " problem=problem,\n", + " model=pod_nn,\n", " optimizer=TorchOptimizer(torch.optim.Adam, lr=0.0001),\n", - " use_lt=False)" + " use_lt=False,\n", + ")" ] }, { @@ -278,14 +279,17 @@ " solver=pod_nn_stokes,\n", " max_epochs=1000,\n", " batch_size=None,\n", - " accelerator='cpu',\n", + " accelerator=\"cpu\",\n", " train_size=0.9,\n", " val_size=0.0,\n", - " test_size=0.1)\n", + " test_size=0.1,\n", + ")\n", "\n", "# fit the pod basis\n", - "trainer.data_module.setup(\"fit\") # set up the dataset\n", - "x_train = trainer.data_module.train_dataset.conditions_dict[\"data\"][\"target\"] # extract data for training\n", + "trainer.data_module.setup(\"fit\") # set up the dataset\n", + "x_train = trainer.data_module.train_dataset.conditions_dict[\"data\"][\n", + " \"target\"\n", + "] # extract data for training\n", "pod_nn.fit_pod(x=x_train)\n", "\n", "# now train\n", @@ -328,12 +332,12 @@ "u_test_nn = pod_nn_stokes(p_test)\n", "u_train_nn = pod_nn_stokes(p_train)\n", "\n", - "relative_error_train = torch.norm(u_train_nn - u_train)/torch.norm(u_train)\n", - "relative_error_test = torch.norm(u_test_nn - u_test)/torch.norm(u_test)\n", + "relative_error_train = torch.norm(u_train_nn - u_train) / torch.norm(u_train)\n", + "relative_error_test = torch.norm(u_test_nn - u_test) / torch.norm(u_test)\n", "\n", - "print('Error summary for POD-NN model:')\n", - "print(f' Train: {relative_error_train.item():e}')\n", - "print(f' Test: {relative_error_test.item():e}')" + "print(\"Error summary for POD-NN model:\")\n", + "print(f\" Train: {relative_error_train.item():e}\")\n", + "print(f\" Test: {relative_error_test.item():e}\")" ] }, { @@ -365,14 +369,11 @@ " \"\"\"\n", "\n", " def __init__(self, pod_rank, rbf_kernel):\n", - " \"\"\"\n", - " \n", - " \"\"\"\n", + " \"\"\" \"\"\"\n", " super().__init__()\n", - " \n", + "\n", " self.pod = PODBlock(pod_rank)\n", " self.rbf = RBFBlock(kernel=rbf_kernel)\n", - " \n", "\n", " def forward(self, x):\n", " \"\"\"\n", @@ -412,7 +413,7 @@ "metadata": {}, "outputs": [], "source": [ - "pod_rbf = PODRBF(pod_rank=20, rbf_kernel='thin_plate_spline')\n", + "pod_rbf = PODRBF(pod_rank=20, rbf_kernel=\"thin_plate_spline\")\n", "pod_rbf.fit(p_train, u_train)" ] }, @@ -444,12 +445,12 @@ "u_test_rbf = pod_rbf(p_test)\n", "u_train_rbf = pod_rbf(p_train)\n", "\n", - "relative_error_train = torch.norm(u_train_rbf - u_train)/torch.norm(u_train)\n", - "relative_error_test = torch.norm(u_test_rbf - u_test)/torch.norm(u_test)\n", + "relative_error_train = torch.norm(u_train_rbf - u_train) / torch.norm(u_train)\n", + "relative_error_test = torch.norm(u_test_rbf - u_test) / torch.norm(u_test)\n", "\n", - "print('Error summary for POD-RBF model:')\n", - "print(f' Train: {relative_error_train.item():e}')\n", - "print(f' Test: {relative_error_test.item():e}')" + "print(\"Error summary for POD-RBF model:\")\n", + "print(f\" Train: {relative_error_train.item():e}\")\n", + "print(f\" Test: {relative_error_test.item():e}\")" ] }, { diff --git a/tutorials/tutorial8/tutorial.py b/tutorials/tutorial8/tutorial.py new file mode 100644 index 000000000..63ca53172 --- /dev/null +++ b/tutorials/tutorial8/tutorial.py @@ -0,0 +1,315 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: Reduced order models (POD-NN and POD-RBF) for parametric problems +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial9/tutorial.ipynb) + +# The tutorial aims to show how to employ the **PINA** library in order to apply a reduced order modeling technique [1]. Such methodologies have several similarities with machine learning approaches, since the main goal consists in predicting the solution of differential equations (typically parametric PDEs) in a real-time fashion. +# +# In particular we are going to use the Proper Orthogonal Decomposition with either Radial Basis Function Interpolation (POD-RBF) or Neural Network (POD-NN) [2]. Here we basically perform a dimensional reduction using the POD approach, approximating the parametric solution manifold (at the reduced space) using a regression technique (NN) and comparing it to an RBF interpolation. In this example, we use a simple multilayer perceptron, but the plenty of different architectures can be plugged as well. + +# Let's start with the necessary imports. +# It's important to note the minimum PINA version to run this tutorial is the `0.1`. + +# In[82]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +get_ipython().run_line_magic("matplotlib", "inline") + +import matplotlib +import matplotlib.pyplot as plt +import torch +import numpy as np +import warnings + +from pina import Trainer +from pina.model import FeedForward +from pina.solver import SupervisedSolver +from pina.optim import TorchOptimizer +from pina.problem.zoo import SupervisedProblem +from pina.model.block import PODBlock, RBFBlock + +warnings.filterwarnings("ignore") + + +# We exploit the [Smithers](https://github.com/mathLab/Smithers) library to collect the parametric snapshots. In particular, we use the `NavierStokesDataset` class that contains a set of parametric solutions of the Navier-Stokes equations in a 2D L-shape domain. The parameter is the inflow velocity. +# The dataset is composed by 500 snapshots of the velocity (along $x$, $y$, and the magnitude) and pressure fields, and the corresponding parameter values. +# +# To visually check the snapshots, let's plot also the data points and the reference solution: this is the expected output of our model. + +# In[83]: + + +from smithers.dataset import NavierStokesDataset + +dataset = NavierStokesDataset() + +fig, axs = plt.subplots(1, 4, figsize=(14, 3)) +for ax, p, u in zip(axs, dataset.params[:4], dataset.snapshots["mag(v)"][:4]): + ax.tricontourf(dataset.triang, u, levels=16) + ax.set_title(f"$\mu$ = {p[0]:.2f}") + + +# The *snapshots* - aka the numerical solutions computed for several parameters - and the corresponding parameters are the only data we need to train the model, in order to predict the solution for any new test parameter. To properly validate the accuracy, we will split the 500 snapshots into the training dataset (90% of the original data) and the testing one (the reamining 10%) inside the `Trainer`. +# +# It is now time to define the problem! + +# In[84]: + + +u = torch.tensor(dataset.snapshots["mag(v)"]).float() +p = torch.tensor(dataset.params).float() +problem = SupervisedProblem(input_=p, output_=u) + + +# We can then build a `POD-NN` model (using an MLP architecture as approximation) and compare it with a `POD-RBF` model (using a Radial Basis Function interpolation as approximation). + +# ## POD-NN reduced order model + +# Let's build the `PODNN` class + +# In[85]: + + +class PODNN(torch.nn.Module): + """ + Proper orthogonal decomposition with neural network model. + """ + + def __init__(self, pod_rank, layers, func): + """ """ + super().__init__() + + self.pod = PODBlock(pod_rank) + self.nn = FeedForward( + input_dimensions=1, + output_dimensions=pod_rank, + layers=layers, + func=func, + ) + + def forward(self, x): + """ + Defines the computation performed at every call. + + :param x: The tensor to apply the forward pass. + :type x: torch.Tensor + :return: the output computed by the model. + :rtype: torch.Tensor + """ + coefficents = self.nn(x) + return self.pod.expand(coefficents) + + def fit_pod(self, x): + """ + Just call the :meth:`pina.model.layers.PODBlock.fit` method of the + :attr:`pina.model.layers.PODBlock` attribute. + """ + self.pod.fit(x) + + +# We highlight that the POD modes are directly computed by means of the singular value decomposition (computed over the input data), and not trained using the backpropagation approach. Only the weights of the MLP are actually trained during the optimization loop. + +# In[86]: + + +pod_nn = PODNN(pod_rank=20, layers=[10, 10, 10], func=torch.nn.Tanh) +pod_nn_stokes = SupervisedSolver( + problem=problem, + model=pod_nn, + optimizer=TorchOptimizer(torch.optim.Adam, lr=0.0001), + use_lt=False, +) + + +# Before starting we need to fit the POD basis on the training dataset, this can be easily done in PINA as well: + +# In[87]: + + +trainer = Trainer( + solver=pod_nn_stokes, + max_epochs=1000, + batch_size=None, + accelerator="cpu", + train_size=0.9, + val_size=0.0, + test_size=0.1, +) + +# fit the pod basis +trainer.data_module.setup("fit") # set up the dataset +x_train = trainer.data_module.train_dataset.conditions_dict["data"][ + "target" +] # extract data for training +pod_nn.fit_pod(x=x_train) + +# now train +trainer.train() + + +# Done! Now that the computational expensive part is over, we can load in future the model to infer new parameters (simply loading the checkpoint file automatically created by `Lightning`) or test its performances. We measure the relative error for the training and test datasets, printing the mean one. + +# In[ ]: + + +# extract train and test data +trainer.data_module.setup("test") # set up the dataset +p_train = trainer.data_module.train_dataset.conditions_dict["data"]["input"] +u_train = trainer.data_module.train_dataset.conditions_dict["data"]["target"] +p_test = trainer.data_module.test_dataset.conditions_dict["data"]["input"] +u_test = trainer.data_module.test_dataset.conditions_dict["data"]["target"] + +# compute statistics +u_test_nn = pod_nn_stokes(p_test) +u_train_nn = pod_nn_stokes(p_train) + +relative_error_train = torch.norm(u_train_nn - u_train) / torch.norm(u_train) +relative_error_test = torch.norm(u_test_nn - u_test) / torch.norm(u_test) + +print("Error summary for POD-NN model:") +print(f" Train: {relative_error_train.item():e}") +print(f" Test: {relative_error_test.item():e}") + + +# ## POD-RBF reduced order model + +# Then, we define the model we want to use, with the POD (`PODBlock`) and the RBF (`RBFBlock`) objects. + +# In[89]: + + +class PODRBF(torch.nn.Module): + """ + Proper orthogonal decomposition with Radial Basis Function interpolation model. + """ + + def __init__(self, pod_rank, rbf_kernel): + """ """ + super().__init__() + + self.pod = PODBlock(pod_rank) + self.rbf = RBFBlock(kernel=rbf_kernel) + + def forward(self, x): + """ + Defines the computation performed at every call. + + :param x: The tensor to apply the forward pass. + :type x: torch.Tensor + :return: the output computed by the model. + :rtype: torch.Tensor + """ + coefficents = self.rbf(x) + return self.pod.expand(coefficents) + + def fit(self, p, x): + """ + Call the :meth:`pina.model.layers.PODBlock.fit` method of the + :attr:`pina.model.layers.PODBlock` attribute to perform the POD, + and the :meth:`pina.model.layers.RBFBlock.fit` method of the + :attr:`pina.model.layers.RBFBlock` attribute to fit the interpolation. + """ + self.pod.fit(x) + self.rbf.fit(p, self.pod.reduce(x)) + + +# We can then fit the model and ask it to predict the required field for unseen values of the parameters. Note that this model does not need a `Trainer` since it does not include any neural network or learnable parameters. + +# In[90]: + + +pod_rbf = PODRBF(pod_rank=20, rbf_kernel="thin_plate_spline") +pod_rbf.fit(p_train, u_train) + + +# Compute errors + +# In[91]: + + +u_test_rbf = pod_rbf(p_test) +u_train_rbf = pod_rbf(p_train) + +relative_error_train = torch.norm(u_train_rbf - u_train) / torch.norm(u_train) +relative_error_test = torch.norm(u_test_rbf - u_test) / torch.norm(u_test) + +print("Error summary for POD-RBF model:") +print(f" Train: {relative_error_train.item():e}") +print(f" Test: {relative_error_test.item():e}") + + +# ## POD-RBF vs POD-NN + +# We can of course also plot the solutions predicted by the `PODRBF` and by the `PODNN` model, comparing them to the original ones. We can note here, in the `PODNN` model and for low velocities, some differences, but improvements can be accomplished thanks to longer training. + +# In[92]: + + +idx = torch.randint(0, len(u_test), (4,)) +u_idx_rbf = pod_rbf(p_test[idx]) +u_idx_nn = pod_nn_stokes(p_test[idx]) + + +fig, axs = plt.subplots(4, 5, figsize=(14, 9)) + +relative_error_rbf = np.abs(u_test[idx] - u_idx_rbf.detach()) +relative_error_rbf = np.where( + u_test[idx] < 1e-7, 1e-7, relative_error_rbf / u_test[idx] +) + +relative_error_nn = np.abs(u_test[idx] - u_idx_nn.detach()) +relative_error_nn = np.where( + u_test[idx] < 1e-7, 1e-7, relative_error_nn / u_test[idx] +) + +for i, (idx_, rbf_, nn_, rbf_err_, nn_err_) in enumerate( + zip(idx, u_idx_rbf, u_idx_nn, relative_error_rbf, relative_error_nn) +): + + axs[0, 0].set_title(f"Real Snapshots") + axs[0, 1].set_title(f"POD-RBF") + axs[0, 2].set_title(f"POD-NN") + axs[0, 3].set_title(f"Error POD-RBF") + axs[0, 4].set_title(f"Error POD-NN") + + cm = axs[i, 0].tricontourf( + dataset.triang, rbf_.detach() + ) # POD-RBF prediction + plt.colorbar(cm, ax=axs[i, 0]) + + cm = axs[i, 1].tricontourf( + dataset.triang, nn_.detach() + ) # POD-NN prediction + plt.colorbar(cm, ax=axs[i, 1]) + + cm = axs[i, 2].tricontourf(dataset.triang, u_test[idx_].flatten()) # Truth + plt.colorbar(cm, ax=axs[i, 2]) + + cm = axs[i, 3].tripcolor( + dataset.triang, rbf_err_, norm=matplotlib.colors.LogNorm() + ) # Error for POD-RBF + plt.colorbar(cm, ax=axs[i, 3]) + + cm = axs[i, 4].tripcolor( + dataset.triang, nn_err_, norm=matplotlib.colors.LogNorm() + ) # Error for POD-NN + plt.colorbar(cm, ax=axs[i, 4]) + +plt.show() + + +# #### References +# 1. Rozza G., Stabile G., Ballarin F. (2022). Advanced Reduced Order Methods and Applications in Computational Fluid Dynamics, Society for Industrial and Applied Mathematics. +# 2. Hesthaven, J. S., & Ubbiali, S. (2018). Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics, 363, 55-78. diff --git a/tutorials/tutorial9/tutorial.ipynb b/tutorials/tutorial9/tutorial.ipynb index 3f19c6c0e..88e1aee51 100644 --- a/tutorials/tutorial9/tutorial.ipynb +++ b/tutorials/tutorial9/tutorial.ipynb @@ -26,12 +26,13 @@ "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", - " import google.colab\n", - " IN_COLAB = True\n", + " import google.colab\n", + "\n", + " IN_COLAB = True\n", "except:\n", - " IN_COLAB = False\n", + " IN_COLAB = False\n", "if IN_COLAB:\n", - " !pip install \"pina-mathlab\"\n", + " !pip install \"pina-mathlab\"\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", @@ -47,7 +48,7 @@ "from pina.equation import Equation\n", "from pina.callback import MetricTracker\n", "\n", - "warnings.filterwarnings('ignore')" + "warnings.filterwarnings(\"ignore\")" ] }, { @@ -111,6 +112,7 @@ " def solution(self, pts):\n", " return torch.sin(torch.pi * pts) * torch.cos(3.0 * torch.pi * pts)\n", "\n", + "\n", "problem = Helmholtz()\n", "\n", "# let's discretise the domain\n", @@ -234,7 +236,7 @@ " pinn,\n", " max_epochs=5000,\n", " accelerator=\"cpu\",\n", - " enable_model_summary=False, \n", + " enable_model_summary=False,\n", " callbacks=[MetricTracker()],\n", " train_size=1.0,\n", " val_size=0.0,\n", @@ -266,9 +268,9 @@ " range(len(trainer_metrics[\"train_loss\"])), trainer_metrics[\"train_loss\"]\n", ")\n", "# plotting\n", - "plt.xlabel('epoch')\n", - "plt.ylabel('loss')\n", - "plt.yscale('log') " + "plt.xlabel(\"epoch\")\n", + "plt.ylabel(\"loss\")\n", + "plt.yscale(\"log\")" ] }, { @@ -305,11 +307,11 @@ } ], "source": [ - "pts = pinn.problem.spatial_domain.sample(256, 'grid', variables='x')\n", - "predicted_output = pinn.forward(pts).extract('u').tensor.detach()\n", + "pts = pinn.problem.spatial_domain.sample(256, \"grid\", variables=\"x\")\n", + "predicted_output = pinn.forward(pts).extract(\"u\").tensor.detach()\n", "true_output = pinn.problem.solution(pts)\n", - "plt.plot(pts.extract(['x']), predicted_output, label='Neural Network solution')\n", - "plt.plot(pts.extract(['x']), true_output, label='True solution')\n", + "plt.plot(pts.extract([\"x\"]), predicted_output, label=\"Neural Network solution\")\n", + "plt.plot(pts.extract([\"x\"]), true_output, label=\"True solution\")\n", "plt.legend()" ] }, @@ -340,21 +342,21 @@ "# plotting solution\n", "with torch.no_grad():\n", " # Notice here we put [-4, 4]!!!\n", - " new_domain = CartesianDomain({'x' : [0, 4]})\n", - " x = new_domain.sample(1000, mode='grid')\n", + " new_domain = CartesianDomain({\"x\": [0, 4]})\n", + " x = new_domain.sample(1000, mode=\"grid\")\n", " fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n", " # Plot 1\n", - " axes[0].plot(x, problem.solution(x), label=r'$u(x)$', color='blue')\n", - " axes[0].set_title(r'True solution $u(x)$')\n", + " axes[0].plot(x, problem.solution(x), label=r\"$u(x)$\", color=\"blue\")\n", + " axes[0].set_title(r\"True solution $u(x)$\")\n", " axes[0].legend(loc=\"upper right\")\n", " # Plot 2\n", - " axes[1].plot(x, pinn(x), label=r'$u_{\\theta}(x)$', color='green')\n", - " axes[1].set_title(r'PINN solution $u_{\\theta}(x)$')\n", + " axes[1].plot(x, pinn(x), label=r\"$u_{\\theta}(x)$\", color=\"green\")\n", + " axes[1].set_title(r\"PINN solution $u_{\\theta}(x)$\")\n", " axes[1].legend(loc=\"upper right\")\n", " # Plot 3\n", " diff = torch.abs(problem.solution(x) - pinn(x))\n", - " axes[2].plot(x, diff, label=r'$|u(x) - u_{\\theta}(x)|$', color='red')\n", - " axes[2].set_title(r'Absolute difference $|u(x) - u_{\\theta}(x)|$')\n", + " axes[2].plot(x, diff, label=r\"$|u(x) - u_{\\theta}(x)|$\", color=\"red\")\n", + " axes[2].set_title(r\"Absolute difference $|u(x) - u_{\\theta}(x)|$\")\n", " axes[2].legend(loc=\"upper right\")\n", " # Adjust layout\n", " plt.tight_layout()\n", diff --git a/tutorials/tutorial9/tutorial.py b/tutorials/tutorial9/tutorial.py new file mode 100644 index 000000000..00073d4e0 --- /dev/null +++ b/tutorials/tutorial9/tutorial.py @@ -0,0 +1,246 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial: One dimensional Helmholtz equation using Periodic Boundary Conditions +# +# [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial9/tutorial.ipynb) +# +# This tutorial presents how to solve with Physics-Informed Neural Networks (PINNs) +# a one dimensional Helmholtz equation with periodic boundary conditions (PBC). +# We will train with standard PINN's training by augmenting the input with +# periodic expansion as presented in [*An expert’s guide to training +# physics-informed neural networks*]( +# https://arxiv.org/abs/2308.08468). +# +# First of all, some useful imports. + +# In[14]: + + +## routine needed to run the notebook on Google Colab +try: + import google.colab + + IN_COLAB = True +except: + IN_COLAB = False +if IN_COLAB: + get_ipython().system('pip install "pina-mathlab"') + +import torch +import matplotlib.pyplot as plt +import warnings + +from pina import Condition, Trainer +from pina.problem import SpatialProblem +from pina.operator import laplacian +from pina.model import FeedForward +from pina.model.block import PeriodicBoundaryEmbedding # The PBC module +from pina.solver import PINN +from pina.domain import CartesianDomain +from pina.equation import Equation +from pina.callback import MetricTracker + +warnings.filterwarnings("ignore") + + +# ## The problem definition +# +# The one-dimensional Helmholtz problem is mathematically written as: +# $$ +# \begin{cases} +# \frac{d^2}{dx^2}u(x) - \lambda u(x) -f(x) &= 0 \quad x\in(0,2)\\ +# u^{(m)}(x=0) - u^{(m)}(x=2) &= 0 \quad m\in[0, 1, \cdots]\\ +# \end{cases} +# $$ +# In this case we are asking the solution to be $C^{\infty}$ periodic with +# period $2$, on the infinite domain $x\in(-\infty, \infty)$. Notice that the +# classical PINN would need infinite conditions to evaluate the PBC loss function, +# one for each derivative, which is of course infeasible... +# A possible solution, diverging from the original PINN formulation, +# is to use *coordinates augmentation*. In coordinates augmentation you seek for +# a coordinates transformation $v$ such that $x\rightarrow v(x)$ such that +# the periodicity condition $ u^{(m)}(x=0) - u^{(m)}(x=2) = 0 \quad m\in[0, 1, \cdots] $ is +# satisfied. +# +# For demonstration purposes, the problem specifics are $\lambda=-10\pi^2$, +# and $f(x)=-6\pi^2\sin(3\pi x)\cos(\pi x)$ which give a solution that can be +# computed analytically $u(x) = \sin(\pi x)\cos(3\pi x)$. + +# In[15]: + + +def helmholtz_equation(input_, output_): + x = input_.extract("x") + u_xx = laplacian(output_, input_, components=["u"], d=["x"]) + f = ( + -6.0 + * torch.pi**2 + * torch.sin(3 * torch.pi * x) + * torch.cos(torch.pi * x) + ) + lambda_ = -10.0 * torch.pi**2 + return u_xx - lambda_ * output_ - f + + +class Helmholtz(SpatialProblem): + output_variables = ["u"] + spatial_domain = CartesianDomain({"x": [0, 2]}) + + # here we write the problem conditions + conditions = { + "phys_cond": Condition( + domain=spatial_domain, equation=Equation(helmholtz_equation) + ), + } + + def solution(self, pts): + return torch.sin(torch.pi * pts) * torch.cos(3.0 * torch.pi * pts) + + +problem = Helmholtz() + +# let's discretise the domain +problem.discretise_domain(200, "grid", domains=["phys_cond"]) + + +# As usual, the Helmholtz problem is written in **PINA** code as a class. +# The equations are written as `conditions` that should be satisfied in the +# corresponding domains. The `solution` +# is the exact solution which will be compared with the predicted one. We used +# Latin Hypercube Sampling for choosing the collocation points. + +# ## Solving the problem with a Periodic Network + +# Any $\mathcal{C}^{\infty}$ periodic function +# $u : \mathbb{R} \rightarrow \mathbb{R}$ with period +# $L\in\mathbb{N}$ can be constructed by composition of an +# arbitrary smooth function $f : \mathbb{R}^n \rightarrow \mathbb{R}$ and a +# given smooth periodic function $v : \mathbb{R} \rightarrow \mathbb{R}^n$ with +# period $L$, that is $u(x) = f(v(x))$. The formulation is generalizable for +# arbitrary dimension, see [*A method for representing periodic functions and +# enforcing exactly periodic boundary conditions with +# deep neural networks*](https://arxiv.org/pdf/2007.07442). +# +# In our case, we rewrite +# $v(x) = \left[1, \cos\left(\frac{2\pi}{L} x\right), +# \sin\left(\frac{2\pi}{L} x\right)\right]$, i.e +# the coordinates augmentation, and $f(\cdot) = NN_{\theta}(\cdot)$ i.e. a neural +# network. The resulting neural network obtained by composing $f$ with $v$ gives +# the PINN approximate solution, that is +# $u(x) \approx u_{\theta}(x)=NN_{\theta}(v(x))$. +# +# In **PINA** this translates in using the `PeriodicBoundaryEmbedding` layer for $v$, and any +# `pina.model` for $NN_{\theta}$. Let's see it in action! +# + +# In[16]: + + +# we encapsulate all modules in a torch.nn.Sequential container +model = torch.nn.Sequential( + PeriodicBoundaryEmbedding(input_dimension=1, periods=2), + FeedForward( + input_dimensions=3, # output of PeriodicBoundaryEmbedding = 3 * input_dimension + output_dimensions=1, + layers=[10, 10], + ), +) + + +# As simple as that! Notice that in higher dimension you can specify different periods +# for all dimensions using a dictionary, e.g. `periods={'x':2, 'y':3, ...}` +# would indicate a periodicity of $2$ in $x$, $3$ in $y$, and so on... +# +# We will now solve the problem as usually with the `PINN` and `Trainer` class, then we will look at the losses using the `MetricTracker` callback from `pina.callback`. + +# In[17]: + + +pinn = PINN( + problem=problem, + model=model, +) +trainer = Trainer( + pinn, + max_epochs=5000, + accelerator="cpu", + enable_model_summary=False, + callbacks=[MetricTracker()], + train_size=1.0, + val_size=0.0, + test_size=0.0, +) +trainer.train() + + +# In[18]: + + +# plot loss +trainer_metrics = trainer.callbacks[0].metrics +plt.plot( + range(len(trainer_metrics["train_loss"])), trainer_metrics["train_loss"] +) +# plotting +plt.xlabel("epoch") +plt.ylabel("loss") +plt.yscale("log") + + +# We are going to plot the solution now! + +# In[19]: + + +pts = pinn.problem.spatial_domain.sample(256, "grid", variables="x") +predicted_output = pinn.forward(pts).extract("u").tensor.detach() +true_output = pinn.problem.solution(pts) +plt.plot(pts.extract(["x"]), predicted_output, label="Neural Network solution") +plt.plot(pts.extract(["x"]), true_output, label="True solution") +plt.legend() + + +# Great, they overlap perfectly! This seems a good result, considering the simple neural network used to some this (complex) problem. We will now test the neural network on the domain $[-4, 4]$ without retraining. In principle the periodicity should be present since the $v$ function ensures the periodicity in $(-\infty, \infty)$. + +# In[20]: + + +# plotting solution +with torch.no_grad(): + # Notice here we put [-4, 4]!!! + new_domain = CartesianDomain({"x": [0, 4]}) + x = new_domain.sample(1000, mode="grid") + fig, axes = plt.subplots(1, 3, figsize=(15, 5)) + # Plot 1 + axes[0].plot(x, problem.solution(x), label=r"$u(x)$", color="blue") + axes[0].set_title(r"True solution $u(x)$") + axes[0].legend(loc="upper right") + # Plot 2 + axes[1].plot(x, pinn(x), label=r"$u_{\theta}(x)$", color="green") + axes[1].set_title(r"PINN solution $u_{\theta}(x)$") + axes[1].legend(loc="upper right") + # Plot 3 + diff = torch.abs(problem.solution(x) - pinn(x)) + axes[2].plot(x, diff, label=r"$|u(x) - u_{\theta}(x)|$", color="red") + axes[2].set_title(r"Absolute difference $|u(x) - u_{\theta}(x)|$") + axes[2].legend(loc="upper right") + # Adjust layout + plt.tight_layout() + # Show the plots + plt.show() + + +# It is pretty clear that the network is periodic, with also the error following a periodic pattern. Obviously a longer training and a more expressive neural network could improve the results! +# +# ## What's next? +# +# Congratulations on completing the one dimensional Helmholtz tutorial of **PINA**! There are multiple directions you can go now: +# +# 1. Train the network for longer or with different layer sizes and assert the finaly accuracy +# +# 2. Apply the `PeriodicBoundaryEmbedding` layer for a time-dependent problem (see reference in the documentation) +# +# 3. Exploit extrafeature training ? +# +# 4. Many more...