From fb16fc7f3ac5c896b691825f66ca941938598948 Mon Sep 17 00:00:00 2001 From: Your Name Date: Mon, 29 Nov 2021 15:29:00 +0100 Subject: [PATCH] tmp commit - toward 0.0.1 --- .gitignore | 139 +++++ README.md | 3 +- examples/run_burgers.py | 115 +++++ ...rametric_elliptic_optimal_control_alpha.py | 207 ++++++++ examples/run_parametric_poisson.py | 157 ++++++ examples/run_poisson.py | 70 +++ pina/__init__.py | 3 + pina/adaptive_functions/__init__.py | 6 + pina/adaptive_functions/adaptive_cos.py | 50 ++ pina/adaptive_functions/adaptive_exp.py | 45 ++ pina/adaptive_functions/adaptive_linear.py | 42 ++ pina/adaptive_functions/adaptive_relu.py | 43 ++ pina/adaptive_functions/adaptive_sin.py | 46 ++ pina/adaptive_functions/adaptive_softplus.py | 42 ++ pina/adaptive_functions/adaptive_square.py | 42 ++ pina/adaptive_functions/adaptive_tanh.py | 52 ++ pina/chebyshev.py | 6 + pina/cube.py | 29 ++ pina/deep_feed_forward.py | 83 +++ pina/label_tensor.py | 49 ++ pina/multi_deep_feed_forward.py | 27 + pina/parametricproblem2d.py | 9 + pina/pinn.py | 488 ++++++++++++++++++ pina/ppinn.py | 465 +++++++++++++++++ pina/problem.py | 49 ++ pina/problem1d.py | 11 + pina/problem2d.py | 16 + pina/segment.py | 27 + pina/tdproblem1d.py | 27 + problems/burgers.py | 78 +++ problems/elliptic_optimal_control.py | 49 ++ .../laplacian_optimal_control_parametric.py | 55 ++ problems/laplacian_parametric.py | 29 ++ .../parametric_elliptic_optimal_control.py | 53 ++ ...elliptic_optimal_control_alpha_variable.py | 52 ++ problems/parametric_poisson.py | 43 ++ problems/poisson2D.py | 41 ++ problems/poisson_2.py | 43 ++ 38 files changed, 2790 insertions(+), 1 deletion(-) create mode 100644 .gitignore create mode 100644 examples/run_burgers.py create mode 100644 examples/run_parametric_elliptic_optimal_control_alpha.py create mode 100644 examples/run_parametric_poisson.py create mode 100644 examples/run_poisson.py create mode 100644 pina/__init__.py create mode 100644 pina/adaptive_functions/__init__.py create mode 100644 pina/adaptive_functions/adaptive_cos.py create mode 100644 pina/adaptive_functions/adaptive_exp.py create mode 100644 pina/adaptive_functions/adaptive_linear.py create mode 100644 pina/adaptive_functions/adaptive_relu.py create mode 100644 pina/adaptive_functions/adaptive_sin.py create mode 100644 pina/adaptive_functions/adaptive_softplus.py create mode 100644 pina/adaptive_functions/adaptive_square.py create mode 100644 pina/adaptive_functions/adaptive_tanh.py create mode 100644 pina/chebyshev.py create mode 100644 pina/cube.py create mode 100644 pina/deep_feed_forward.py create mode 100644 pina/label_tensor.py create mode 100644 pina/multi_deep_feed_forward.py create mode 100644 pina/parametricproblem2d.py create mode 100644 pina/pinn.py create mode 100644 pina/ppinn.py create mode 100644 pina/problem.py create mode 100644 pina/problem1d.py create mode 100644 pina/problem2d.py create mode 100644 pina/segment.py create mode 100644 pina/tdproblem1d.py create mode 100644 problems/burgers.py create mode 100644 problems/elliptic_optimal_control.py create mode 100644 problems/laplacian_optimal_control_parametric.py create mode 100644 problems/laplacian_parametric.py create mode 100644 problems/parametric_elliptic_optimal_control.py create mode 100644 problems/parametric_elliptic_optimal_control_alpha_variable.py create mode 100644 problems/parametric_poisson.py create mode 100644 problems/poisson2D.py create mode 100644 problems/poisson_2.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..73358ad --- /dev/null +++ b/.gitignore @@ -0,0 +1,139 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + diff --git a/README.md b/README.md index cbab1f6..be7a139 100644 --- a/README.md +++ b/README.md @@ -1 +1,2 @@ -void +# PINA +# Physic Informed Neural networks for Advance modeling diff --git a/examples/run_burgers.py b/examples/run_burgers.py new file mode 100644 index 0000000..059243c --- /dev/null +++ b/examples/run_burgers.py @@ -0,0 +1,115 @@ +import sys +import numpy as np +import torch +import argparse +from pina.pinn import PINN +from pina.ppinn import ParametricPINN as pPINN +from pina.label_tensor import LabelTensor +from torch.nn import ReLU, Tanh, Softplus +from problems.burgers import Burgers1D +from pina.deep_feed_forward import DeepFeedForward + +from pina.adaptive_functions import AdaptiveSin, AdaptiveCos, AdaptiveTanh + + +class myFeature(torch.nn.Module): + """ + Feature: sin(x) + """ + + def __init__(self, idx): + + super(myFeature, self).__init__() + self.idx = idx + + def forward(self, x): + + return torch.sin(torch.pi * x[:, self.idx]) + +class myExp(torch.nn.Module): + def __init__(self, idx): + + super().__init__() + self.idx = idx + + def forward(self, x): + + return torch.exp(x[:, self.idx]) + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Run PINA") + group = parser.add_mutually_exclusive_group(required=True) + group.add_argument("-s", "-save", action="store_true") + group.add_argument("-l", "-load", action="store_true") + parser.add_argument("id_run", help="number of run", type=int) + parser.add_argument("features", help="extra features", type=int) + args = parser.parse_args() + + feat = [myFeature(0)] if args.features else [] + + burgers_problem = Burgers1D() + model = DeepFeedForward( + layers=[20, 10, 5], + #layers=[8, 4, 2], + #layers=[16, 8, 4, 4], + #layers=[20, 4, 4, 4], + output_variables=burgers_problem.output_variables, + input_variables=burgers_problem.input_variables, + func=Tanh, + extra_features=feat + ) + + pinn = PINN( + burgers_problem, + model, + lr=0.006, + error_norm='mse', + regularizer=0, + lr_accelerate=None) + + if args.s: + + pinn.span_pts(8000, 'latin', ['D']) + pinn.span_pts(50, 'random', ['gamma1', 'gamma2', 'initia']) + #pinn.plot_pts() + pinn.train(10000, 1000) + #with open('burgers_history_{}_{}.txt'.format(args.id_run, args.features), 'w') as file_: + # for i, losses in enumerate(pinn.history): + # file_.write('{} {}\n'.format(i, sum(losses).item())) + pinn.save_state('pina.burger.{}.{}'.format(args.id_run, args.features)) + + else: + pinn.load_state('pina.burger.{}.{}'.format(args.id_run, args.features)) + #pinn.plot(256,filename='pina.burger.{}.{}.jpg'.format(args.id_run, args.features)) + + + print(pinn.history) + with open('burgers_history_{}_{}.txt'.format(args.id_run, args.features), 'w') as file_: + for i, losses in enumerate(pinn.history): + print(losses) + file_.write('{} {}\n'.format(i, sum(losses))) + import scipy + data = scipy.io.loadmat('Data/burgers_shock.mat') + data_solution = {'grid': np.meshgrid(data['x'], data['t']), 'grid_solution': data['usol'].T} + import matplotlib + matplotlib.use('Qt5Agg') + import matplotlib.pyplot as plt + + t =75 + for t in [25, 50, 75]: + input = torch.cat([ + torch.linspace(-1, 1, 256).reshape(-1, 1), + torch.ones(size=(256, 1)) * t /100], + axis=1).double() + output = pinn.model(input) + fout = 'pina.burgers.{}.{}.t{}.dat'.format(args.id_run, args.features, t) + with open(fout, 'w') as f_: + f_.write('x utruth upinn\n') + for x, utruth, upinn in zip(data['x'], data['usol'][:, t], output.tensor.detach()): + f_.write('{} {} {}\n'.format(x[0], utruth, upinn.item())) + plt.plot(data['usol'][:, t], label='truth') + plt.plot(output.tensor.detach(), 'x', label='pinn') + plt.legend() + plt.show() diff --git a/examples/run_parametric_elliptic_optimal_control_alpha.py b/examples/run_parametric_elliptic_optimal_control_alpha.py new file mode 100644 index 0000000..b3b73b5 --- /dev/null +++ b/examples/run_parametric_elliptic_optimal_control_alpha.py @@ -0,0 +1,207 @@ +import numpy as np +import torch +import argparse +from pina.pinn import PINN +from pina.ppinn import ParametricPINN as pPINN +from pina.label_tensor import LabelTensor +from torch.nn import ReLU, Tanh, Softplus +from pina.adaptive_functions.adaptive_softplus import AdaptiveSoftplus +from problems.parametric_elliptic_optimal_control_alpha_variable import ParametricEllipticOptimalControl +from pina.multi_deep_feed_forward import MultiDeepFeedForward +from pina.deep_feed_forward import DeepFeedForward + +alpha = 1 + +class myFeature(torch.nn.Module): + """ + Feature: sin(x) + """ + + def __init__(self): + super(myFeature, self).__init__() + + def forward(self, x): + return (-x[:, 0]**2+1) * (-x[:, 1]**2+1) + + +class CustomMultiDFF(MultiDeepFeedForward): + + def __init__(self, dff_dict): + super().__init__(dff_dict) + + def forward(self, x): + out = self.uu(x) + p = LabelTensor((out['u_param'] * x[:, 3]).reshape(-1, 1), ['p']) + a = LabelTensor.hstack([out, p]) + return a + +model = CustomMultiDFF( + { + 'uu': { + 'input_variables': ['x1', 'x2', 'mu', 'alpha'], + 'output_variables': ['u_param', 'y'], + 'layers': [40, 40, 20], + 'func': Softplus, + 'extra_features': [myFeature()], + }, + # 'u_param': { + # 'input_variables': ['u', 'mu'], + # 'output_variables': ['u_param'], + # 'layers': [], + # 'func': None + # }, + # 'p': { + # 'input_variables': ['u'], + # 'output_variables': ['p'], + # 'layers': [10], + # 'func': None + # }, + } +) + + +opc = ParametricEllipticOptimalControl(alpha) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Run PINA") + group = parser.add_mutually_exclusive_group(required=True) + group.add_argument("-s", "-save", action="store_true") + group.add_argument("-l", "-load", action="store_true") + args = parser.parse_args() + + # model = DeepFeedForward( + # layers=[40, 40, 20], + # output_variables=['u_param', 'y', 'p'], + # input_variables=opc.input_variables+['mu', 'alpha'], + # func=Softplus, + # extra_features=[myFeature()] + # ) + + + pinn = pPINN( + opc, + model, + lr=0.002, + error_norm='mse', + regularizer=1e-8, + lr_accelerate=None) + + if args.s: + + pinn.span_pts(30, 'grid', ['D1']) + pinn.span_pts(50, 'grid', ['gamma1', 'gamma2', 'gamma3', 'gamma4']) + pinn.train(10000, 20) + # with open('ocp_wrong_history.txt', 'w') as file_: + # for i, losses in enumerate(pinn.history): + # file_.write('{} {}\n'.format(i, sum(losses).item())) + + pinn.save_state('pina.ocp') + + else: + pinn.load_state('working.pina.ocp') + pinn.load_state('pina.ocp') + + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + # res = 64 + # param = torch.tensor([[3., 1]]) + # pts_container = [] + # for mn, mx in [[-1, 1], [-1, 1]]: + # pts_container.append(np.linspace(mn, mx, res)) + # grids_container = np.meshgrid(*pts_container) + # unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + # unrolled_pts = torch.cat([unrolled_pts, param.double().repeat(unrolled_pts.shape[0], 1).reshape(-1, 2)], axis=1) + + # unrolled_pts = LabelTensor(unrolled_pts, ['x1', 'x2', 'mu', 'alpha']) + # Z_pred = pinn.model(unrolled_pts.tensor) + # print(Z_pred.tensor.shape) + + # plt.subplot(2, 3, 1) + # plt.pcolor(Z_pred['y'].reshape(res, res).detach()) + # plt.colorbar() + # plt.subplot(2, 3, 2) + # plt.pcolor(Z_pred['u_param'].reshape(res, res).detach()) + # plt.colorbar() + # plt.subplot(2, 3, 3) + # plt.pcolor(Z_pred['p'].reshape(res, res).detach()) + # plt.colorbar() + # with open('ocp_mu3_a1_plot.txt', 'w') as f_: + # f_.write('x y u p ys\n') + # for (x, y), tru, pre, e in zip(unrolled_pts[:, :2], + # Z_pred['u_param'].reshape(-1, 1), + # Z_pred['p'].reshape(-1, 1), + # Z_pred['y'].reshape(-1, 1), + # ): + # f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + + + # param = torch.tensor([[3.0, 0.01]]) + # unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + # unrolled_pts = torch.cat([unrolled_pts, param.double().repeat(unrolled_pts.shape[0], 1).reshape(-1, 2)], axis=1) + # unrolled_pts = LabelTensor(unrolled_pts, ['x1', 'x2', 'mu', 'alpha']) + # Z_pred = pinn.model(unrolled_pts.tensor) + + # plt.subplot(2, 3, 4) + # plt.pcolor(Z_pred['y'].reshape(res, res).detach()) + # plt.colorbar() + # plt.subplot(2, 3, 5) + # plt.pcolor(Z_pred['u_param'].reshape(res, res).detach()) + # plt.colorbar() + # plt.subplot(2, 3, 6) + # plt.pcolor(Z_pred['p'].reshape(res, res).detach()) + # plt.colorbar() + + # plt.show() + # with open('ocp_mu3_a0.01_plot.txt', 'w') as f_: + # f_.write('x y u p ys\n') + # for (x, y), tru, pre, e in zip(unrolled_pts[:, :2], + # Z_pred['u_param'].reshape(-1, 1), + # Z_pred['p'].reshape(-1, 1), + # Z_pred['y'].reshape(-1, 1), + # ): + # f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + + + + + y = {} + u = {} + for alpha in [0.01, 0.1, 1]: + y[alpha] = [] + u[alpha] = [] + for p in np.linspace(0.5, 3, 32): + a = pinn.model(LabelTensor(torch.tensor([[0, 0, p, alpha]]).double(), ['x1', 'x2', 'mu', 'alpha']).tensor) + y[alpha].append(a['y'].detach().numpy()[0]) + u[alpha].append(a['u_param'].detach().numpy()[0]) + + + + plt.plot(np.linspace(0.5, 3, 32), u[1], label='u') + plt.plot(np.linspace(0.5, 3, 32), u[0.01], label='u') + plt.plot(np.linspace(0.5, 3, 32), u[0.1], label='u') + plt.plot([1, 2, 3], [0.28, 0.56, 0.85], 'o', label='Truth values') + plt.legend() + plt.show() + print(y[1]) + print(y[0.1]) + print(y[0.01]) + with open('elliptic_param_y.txt', 'w') as f_: + f_.write('mu 1 01 001\n') + for mu, y1, y01, y001 in zip(np.linspace(0.5, 3, 32), y[1], y[0.1], y[0.01]): + f_.write('{} {} {} {}\n'.format(mu, y1, y01, y001)) + + with open('elliptic_param_u.txt', 'w') as f_: + f_.write('mu 1 01 001\n') + for mu, y1, y01, y001 in zip(np.linspace(0.5, 3, 32), u[1], u[0.1], u[0.01]): + f_.write('{} {} {} {}\n'.format(mu, y1, y01, y001)) + + + plt.plot(np.linspace(0.5, 3, 32), y, label='y') + plt.plot([1, 2, 3], [0.062, 0.12, 0.19], 'o', label='Truth values') + plt.legend() + plt.show() + + diff --git a/examples/run_parametric_poisson.py b/examples/run_parametric_poisson.py new file mode 100644 index 0000000..f25cb87 --- /dev/null +++ b/examples/run_parametric_poisson.py @@ -0,0 +1,157 @@ +import sys +import numpy as np +import torch +import argparse +from pina.pinn import PINN +from pina.ppinn import ParametricPINN as pPINN +from pina.label_tensor import LabelTensor +from torch.nn import ReLU, Tanh, Softplus +from problems.parametric_poisson import ParametricPoisson2DProblem as Poisson2D +from pina.deep_feed_forward import DeepFeedForward + +from pina.adaptive_functions import AdaptiveSin, AdaptiveCos, AdaptiveTanh + + +class myFeature(torch.nn.Module): + """ + Feature: sin(x) + """ + + def __init__(self): + super(myFeature, self).__init__() + + def forward(self, x): + return torch.exp(- 2*(x[:, 0] - x[:, 2])**2 - 2*(x[:, 1] - x[:, 3])**2) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Run PINA") + group = parser.add_mutually_exclusive_group(required=True) + group.add_argument("-s", "-save", action="store_true") + group.add_argument("-l", "-load", action="store_true") + parser.add_argument("id_run", help="number of run", type=int) + parser.add_argument("features", help="extra features", type=int) + args = parser.parse_args() + + feat = [myFeature()] if args.features else [] + + poisson_problem = Poisson2D() + model = DeepFeedForward( + layers=[200, 40, 10], + output_variables=poisson_problem.output_variables, + input_variables=poisson_problem.input_variables+['mu1', 'mu2'], + func=Softplus, + extra_features=feat + ) + + pinn = pPINN( + poisson_problem, + model, + lr=0.0006, + regularizer=1e-6, + lr_accelerate=None) + + if args.s: + + pinn.span_pts(30, 'chebyshev', ['D']) + pinn.span_pts(50, 'grid', ['gamma1', 'gamma2', 'gamma3', 'gamma4']) + #pinn.plot_pts() + pinn.train(10000, 10) + pinn.save_state('pina.poisson_param') + + else: + pinn.load_state('pina.poisson_param') + #pinn.plot(40, torch.tensor([-0.8, -0.8])) + #pinn.plot(40, torch.tensor([ 0.8, 0.8])) + + from smithers.io import VTUHandler + from scipy.interpolate import griddata + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + res = 64 + fname_minus = 'Poisson_param_08minus000000.vtu' + param = torch.tensor([-0.8, -0.8]) + pts_container = [] + for mn, mx in [[-1, 1], [-1, 1]]: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + unrolled_pts = torch.cat([unrolled_pts, param.double().repeat(unrolled_pts.shape[0]).reshape(-1, 2)], axis=1) + + #unrolled_pts.to(dtype=self.dtype) + unrolled_pts = LabelTensor(unrolled_pts, ['x1', 'x2', 'mu1', 'mu2']) + + Z_pred = pinn.model(unrolled_pts.tensor) + data = VTUHandler.read(fname_minus) + + + print(data['points'][:, :2].shape) + print(data['point_data']['f_16'].shape) + print(grids_container[0].shape) + print(grids_container[1].shape) + Z_truth = griddata(data['points'][:, :2], data['point_data']['f_16'], (grids_container[0], grids_container[1])) + + + err = np.abs(Z_truth + Z_pred.tensor.reshape(res, res).detach().numpy()) + + plt.subplot(1, 3, 1) + plt.pcolor(-Z_pred.tensor.reshape(res, res).detach()) + plt.colorbar() + plt.subplot(1, 3, 2) + plt.pcolor(Z_truth) + plt.colorbar() + plt.subplot(1, 3, 3) + plt.pcolor(err, vmin=0, vmax=0.009) + plt.colorbar() + plt.show() + + print(unrolled_pts.tensor.shape) + with open('parpoisson_minus_plot.txt', 'w') as f_: + f_.write('x y truth pred e\n') + for (x, y), tru, pre, e in zip(unrolled_pts[:, :2], Z_truth.reshape(-1, 1), -Z_pred.tensor.reshape(-1, 1), err.reshape(-1, 1)): + f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + + fname_plus = 'Poisson_param_08plus000000.vtu' + param = torch.tensor([0.8, 0.8]) + pts_container = [] + for mn, mx in [[-1, 1], [-1, 1]]: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + unrolled_pts = torch.cat([unrolled_pts, param.double().repeat(unrolled_pts.shape[0]).reshape(-1, 2)], axis=1) + + #unrolled_pts.to(dtype=self.dtype) + unrolled_pts = LabelTensor(unrolled_pts, ['x1', 'x2', 'mu1', 'mu2']) + + Z_pred = pinn.model(unrolled_pts.tensor) + data = VTUHandler.read(fname_plus) + + + print(data['points'][:, :2].shape) + print(data['point_data']['f_16'].shape) + print(grids_container[0].shape) + print(grids_container[1].shape) + Z_truth = griddata(data['points'][:, :2], data['point_data']['f_16'], (grids_container[0], grids_container[1])) + + + err = np.abs(Z_truth + Z_pred.tensor.reshape(res, res).detach().numpy()) + + plt.subplot(1, 3, 1) + plt.pcolor(-Z_pred.tensor.reshape(res, res).detach()) + plt.colorbar() + plt.subplot(1, 3, 2) + plt.pcolor(Z_truth) + plt.colorbar() + plt.subplot(1, 3, 3) + print('gggggg') + plt.pcolor(err, vmin=0, vmax=0.001) + plt.colorbar() + plt.show() + with open('parpoisson_plus_plot.txt', 'w') as f_: + f_.write('x y truth pred e\n') + for (x, y), tru, pre, e in zip(unrolled_pts[:, :2], Z_truth.reshape(-1, 1), -Z_pred.tensor.reshape(-1, 1), err.reshape(-1, 1)): + f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + + diff --git a/examples/run_poisson.py b/examples/run_poisson.py new file mode 100644 index 0000000..5917ffd --- /dev/null +++ b/examples/run_poisson.py @@ -0,0 +1,70 @@ +import sys +import numpy as np +import torch +import argparse +from pina.pinn import PINN +from pina.ppinn import ParametricPINN as pPINN +from pina.label_tensor import LabelTensor +from torch.nn import ReLU, Tanh, Softplus +from problems.poisson2D import Poisson2DProblem as Poisson2D +from pina.deep_feed_forward import DeepFeedForward + +from pina.adaptive_functions import AdaptiveSin, AdaptiveCos, AdaptiveTanh + + +class myFeature(torch.nn.Module): + """ + Feature: sin(x) + """ + + def __init__(self): + super(myFeature, self).__init__() + + def forward(self, x): + return torch.sin(x[:, 0]*torch.pi) * torch.sin(x[:, 1]*torch.pi) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Run PINA") + group = parser.add_mutually_exclusive_group(required=True) + group.add_argument("-s", "-save", action="store_true") + group.add_argument("-l", "-load", action="store_true") + parser.add_argument("id_run", help="number of run", type=int) + parser.add_argument("features", help="extra features", type=int) + args = parser.parse_args() + + feat = [myFeature()] if args.features else [] + + poisson_problem = Poisson2D() + model = DeepFeedForward( + layers=[10, 10], + output_variables=poisson_problem.output_variables, + input_variables=poisson_problem.input_variables, + func=Softplus, + extra_features=feat + ) + + pinn = PINN( + poisson_problem, + model, + lr=0.003, + error_norm='mse', + regularizer=1e-8, + lr_accelerate=None) + + if args.s: + + pinn.span_pts(10, 'grid', ['D']) + pinn.span_pts(10, 'grid', ['gamma1', 'gamma2', 'gamma3', 'gamma4']) + #pinn.plot_pts() + pinn.train(10000, 100) + with open('poisson_history_{}_{}.txt'.format(args.id_run, args.features), 'w') as file_: + for i, losses in enumerate(pinn.history): + file_.write('{} {}\n'.format(i, sum(losses).item())) + pinn.save_state('pina.poisson') + + else: + pinn.load_state('pina.poisson') + pinn.plot(40) + + diff --git a/pina/__init__.py b/pina/__init__.py new file mode 100644 index 0000000..c0ec9ec --- /dev/null +++ b/pina/__init__.py @@ -0,0 +1,3 @@ +from .pinn import PINN +from .deep_feed_forward import DeepFeedForward +from .problem1d import Problem1D diff --git a/pina/adaptive_functions/__init__.py b/pina/adaptive_functions/__init__.py new file mode 100644 index 0000000..adc0f17 --- /dev/null +++ b/pina/adaptive_functions/__init__.py @@ -0,0 +1,6 @@ + +from .adaptive_tanh import AdaptiveTanh +from .adaptive_sin import AdaptiveSin +from .adaptive_cos import AdaptiveCos +from .adaptive_linear import AdaptiveLinear +from .adaptive_square import AdaptiveSquare diff --git a/pina/adaptive_functions/adaptive_cos.py b/pina/adaptive_functions/adaptive_cos.py new file mode 100644 index 0000000..0bae214 --- /dev/null +++ b/pina/adaptive_functions/adaptive_cos.py @@ -0,0 +1,50 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveCos(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self, alpha = None): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveCos, self).__init__() + #self.in_features = in_features + + # initialize alpha + if alpha == None: + self.alpha = Parameter(torch.tensor(1.0)) # create a tensor out of alpha + else: + self.alpha = Parameter(torch.tensor(alpha)) # create a tensor out of alpha + self.alpha.requiresGrad = True # set requiresGrad to true! + + self.scale = Parameter(torch.tensor(1.0)) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.tensor(0.0)) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + return self.scale * (torch.cos(self.alpha * x + self.translate)) diff --git a/pina/adaptive_functions/adaptive_exp.py b/pina/adaptive_functions/adaptive_exp.py new file mode 100644 index 0000000..d9b1732 --- /dev/null +++ b/pina/adaptive_functions/adaptive_exp.py @@ -0,0 +1,45 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveExp(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveExp, self).__init__() + + self.scale = Parameter(torch.normal(torch.tensor(1.0), torch.tensor(0.1))) # create a tensor out of alpha + self.scale.requiresGrad = True # set requiresGrad to true! + + self.alpha = Parameter(torch.normal(torch.tensor(1.0), torch.tensor(0.1))) # create a tensor out of alpha + self.alpha.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.normal(torch.tensor(0.0), torch.tensor(0.1))) # create a tensor out of alpha + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + return self.scale * (x + self.translate) diff --git a/pina/adaptive_functions/adaptive_linear.py b/pina/adaptive_functions/adaptive_linear.py new file mode 100644 index 0000000..ffe24f6 --- /dev/null +++ b/pina/adaptive_functions/adaptive_linear.py @@ -0,0 +1,42 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveLinear(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveLinear, self).__init__() + + self.scale = Parameter(torch.tensor(1.0)) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.tensor(0.0)) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + return self.scale * (x + self.translate) diff --git a/pina/adaptive_functions/adaptive_relu.py b/pina/adaptive_functions/adaptive_relu.py new file mode 100644 index 0000000..0e99ba9 --- /dev/null +++ b/pina/adaptive_functions/adaptive_relu.py @@ -0,0 +1,43 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveReLU(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveReLU, self).__init__() + + self.scale = Parameter(torch.rand(1)) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.rand(1)) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + #x += self.translate + return torch.relu(x+self.translate)*self.scale diff --git a/pina/adaptive_functions/adaptive_sin.py b/pina/adaptive_functions/adaptive_sin.py new file mode 100644 index 0000000..80793c2 --- /dev/null +++ b/pina/adaptive_functions/adaptive_sin.py @@ -0,0 +1,46 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveSin(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self, alpha = None): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveSin, self).__init__() + + # initialize alpha + self.alpha = Parameter(torch.normal(torch.tensor(1.0), torch.tensor(0.1))) # create a tensor out of alpha + self.alpha.requiresGrad = True # set requiresGrad to true! + + self.scale = Parameter(torch.normal(torch.tensor(1.0), torch.tensor(0.1))) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.normal(torch.tensor(0.0), torch.tensor(0.1))) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + return self.scale * (torch.sin(self.alpha * x + self.translate)) diff --git a/pina/adaptive_functions/adaptive_softplus.py b/pina/adaptive_functions/adaptive_softplus.py new file mode 100644 index 0000000..c2069ff --- /dev/null +++ b/pina/adaptive_functions/adaptive_softplus.py @@ -0,0 +1,42 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveSoftplus(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super().__init__() + + self.soft = torch.nn.Softplus() + + self.scale = Parameter(torch.rand(1)) + self.scale.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + #x += self.translate + return self.soft(x)*self.scale diff --git a/pina/adaptive_functions/adaptive_square.py b/pina/adaptive_functions/adaptive_square.py new file mode 100644 index 0000000..911673a --- /dev/null +++ b/pina/adaptive_functions/adaptive_square.py @@ -0,0 +1,42 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveSquare(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self, alpha = None): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveSquare, self).__init__() + + self.scale = Parameter(torch.tensor(1.0)) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.tensor(0.0)) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + return self.scale * (x + self.translate)**2 diff --git a/pina/adaptive_functions/adaptive_tanh.py b/pina/adaptive_functions/adaptive_tanh.py new file mode 100644 index 0000000..5b6bede --- /dev/null +++ b/pina/adaptive_functions/adaptive_tanh.py @@ -0,0 +1,52 @@ +import torch +from torch.nn.parameter import Parameter + +class AdaptiveTanh(torch.nn.Module): + ''' + Implementation of soft exponential activation. + Shape: + - Input: (N, *) where * means, any number of additional + dimensions + - Output: (N, *), same shape as the input + Parameters: + - alpha - trainable parameter + References: + - See related paper: + https://arxiv.org/pdf/1602.01321.pdf + Examples: + >>> a1 = soft_exponential(256) + >>> x = torch.randn(256) + >>> x = a1(x) + ''' + def __init__(self, alpha = None): + ''' + Initialization. + INPUT: + - in_features: shape of the input + - aplha: trainable parameter + aplha is initialized with zero value by default + ''' + super(AdaptiveTanh, self).__init__() + #self.in_features = in_features + + # initialize alpha + if alpha == None: + self.alpha = Parameter(torch.tensor(1.0)) # create a tensor out of alpha + else: + self.alpha = Parameter(torch.tensor(alpha)) # create a tensor out of alpha + + self.alpha.requiresGrad = True # set requiresGrad to true! + + self.scale = Parameter(torch.tensor(1.0)) + self.scale.requiresGrad = True # set requiresGrad to true! + + self.translate = Parameter(torch.tensor(0.0)) + self.translate.requiresGrad = True # set requiresGrad to true! + + def forward(self, x): + ''' + Forward pass of the function. + Applies the function to the input elementwise. + ''' + x += self.translate + return self.scale * (torch.exp(self.alpha * x) - torch.exp(-self.alpha * x))/(torch.exp(self.alpha * x) + torch.exp(-self.alpha * x)) diff --git a/pina/chebyshev.py b/pina/chebyshev.py new file mode 100644 index 0000000..6a9fa15 --- /dev/null +++ b/pina/chebyshev.py @@ -0,0 +1,6 @@ +from mpmath import chebyt, chop, taylor +import numpy as np +def chebyshev_roots(n): + """ Return the roots of *n* Chebyshev polynomials (between [-1, 1]) """ + return np.roots(chop(taylor(lambda x: chebyt(n, x), 0, n))[::-1]) + diff --git a/pina/cube.py b/pina/cube.py new file mode 100644 index 0000000..9c702a9 --- /dev/null +++ b/pina/cube.py @@ -0,0 +1,29 @@ +import numpy as np +from .chebyshev import chebyshev_roots + +class Cube(): + def __init__(self, bound): + self.bound = np.asarray(bound) + + def discretize(self, n, mode='random'): + + if mode == 'random': + pts = np.random.uniform(size=(n, self.bound.shape[0])) + elif mode == 'chebyshev': + pts = np.array([chebyshev_roots(n) *.5 + .5 for _ in range(self.bound.shape[0])]) + grids = np.meshgrid(*pts) + pts = np.hstack([grid.reshape(-1, 1) for grid in grids]) + elif mode == 'grid': + pts = np.array([np.linspace(0, 1, n) for _ in range(self.bound.shape[0])]) + grids = np.meshgrid(*pts) + pts = np.hstack([grid.reshape(-1, 1) for grid in grids]) + elif mode == 'lh' or mode == 'latin': + from scipy.stats import qmc + sampler = qmc.LatinHypercube(d=self.bound.shape[0]) + pts = sampler.random(n) + + # Scale pts + pts *= self.bound[:, 1] - self.bound[:, 0] + pts += self.bound[:, 0] + + return pts diff --git a/pina/deep_feed_forward.py b/pina/deep_feed_forward.py new file mode 100644 index 0000000..6723160 --- /dev/null +++ b/pina/deep_feed_forward.py @@ -0,0 +1,83 @@ + +from .problem import Problem +import torch +import torch.nn as nn +import numpy as np +from .cube import Cube +torch.pi = torch.acos(torch.zeros(1)).item() * 2 # which is 3.1415927410125732 +from torch.nn import Tanh, ReLU +#import torch.nn.utils.prune as prune +from pina.adaptive_functions import AdaptiveLinear +from pina.label_tensor import LabelTensor + +class DeepFeedForward(torch.nn.Module): + + def __init__(self, + inner_size=20, + n_layers=2, + func=nn.Tanh, + input_variables=None, + output_variables=None, + layers=None, + extra_features=None): + ''' + ''' + super(DeepFeedForward, self).__init__() + + if extra_features is None: + extra_features = [] + self.extra_features = nn.Sequential(*extra_features) + + if input_variables is None: input_variables = ['x'] + if output_variables is None: input_variables = ['y'] + + self.input_variables = input_variables + self.input_dimension = len(input_variables) + + self.output_variables = output_variables + self.output_dimension = len(output_variables) + + n_features = len(extra_features) + + if layers is None: layers = [inner_size] * n_layers + + tmp_layers = layers.copy() + tmp_layers.insert(0, self.input_dimension+n_features)#-1) + tmp_layers.append(self.output_dimension) + + self.layers = [] + for i in range(len(tmp_layers)-1): + self.layers.append(nn.Linear(tmp_layers[i], tmp_layers[i+1])) + + + + if isinstance(func, list): + self.functions = func + else: + self.functions = [func for _ in range(len(self.layers)-1)] + + + unique_list = [] + for layer, func in zip(self.layers[:-1], self.functions): + unique_list.append(layer) + if func is not None: unique_list.append(func()) + unique_list.append(self.layers[-1]) + + self.model = nn.Sequential(*unique_list) + + + def forward(self, x): + + nf = len(self.extra_features) + if nf == 0: + return LabelTensor(self.model(x), self.output_variables) + + # if self.extra_features + #input_ = torch.zeros(x.shape[0], nf+self.input_dimension, dtype=x.dtype, device=x.device) + input_ = torch.zeros(x.shape[0], nf+x.shape[1], dtype=x.dtype, device=x.device) + input_[:, :x.shape[1]] = x + for i, feature in enumerate(self.extra_features, start=self.input_dimension): + input_[:, i] = feature(x) + return LabelTensor(self.model(input_), self.output_variables) + + diff --git a/pina/label_tensor.py b/pina/label_tensor.py new file mode 100644 index 0000000..00a6448 --- /dev/null +++ b/pina/label_tensor.py @@ -0,0 +1,49 @@ +import torch + +class LabelTensor(): + + def __init__(self, x, labels): + + + if len(labels) != x.shape[1]: + print(len(labels), x.shape[1]) + raise ValueError + self.__labels = labels + self.tensor = x + + def __getitem__(self, key): + if key in self.labels: + return self.tensor[:, self.labels.index(key)] + else: + return self.tensor.__getitem__(key) + + def __repr__(self): + return self.tensor + + def __str__(self): + return self.tensor, self.labels + + @property + def labels(self): + return self.__labels + + @staticmethod + def hstack(labeltensor_list): + concatenated_tensor = torch.cat([lt.tensor for lt in labeltensor_list], axis=1) + concatenated_label = sum([lt.labels for lt in labeltensor_list], []) + return LabelTensor(concatenated_tensor, concatenated_label) + + + +if __name__ == "__main__": + import numpy as np + a = np.random.uniform(size=(20, 3)) + a = np.random.uniform(size=(20, 3)) + p = torch.from_numpy(a) + t = LabelTensor(p, labels=['u', 'p', 't']) + print(t) + print(t['u']) + t *= 2 + print(t['u']) + print(t[:, 0]) + diff --git a/pina/multi_deep_feed_forward.py b/pina/multi_deep_feed_forward.py new file mode 100644 index 0000000..3bde1e5 --- /dev/null +++ b/pina/multi_deep_feed_forward.py @@ -0,0 +1,27 @@ + +from .problem import Problem +import torch +import torch.nn as nn +import numpy as np +from .cube import Cube +torch.pi = torch.acos(torch.zeros(1)).item() * 2 # which is 3.1415927410125732 +from torch.nn import Tanh, ReLU +import torch.nn.utils.prune as prune +from pina.adaptive_functions import AdaptiveLinear +from pina.deep_feed_forward import DeepFeedForward + +class MultiDeepFeedForward(torch.nn.Module): + + def __init__(self, dff_dict): + ''' + ''' + super().__init__() + + if not isinstance(dff_dict, dict): + raise TypeError + + for name, constructor_args in dff_dict.items(): + setattr(self, name, DeepFeedForward(**constructor_args)) + + + diff --git a/pina/parametricproblem2d.py b/pina/parametricproblem2d.py new file mode 100644 index 0000000..beefd35 --- /dev/null +++ b/pina/parametricproblem2d.py @@ -0,0 +1,9 @@ +from .problem2d import Problem2D +import numpy as np + +class ParametricProblem2D(Problem2D): + + def __init__(self, variables=None, bc=None, params_bound=None, domain_bound=None): + + Problem2D.__init__(self, variables=variables, bc=bc, domain_bound=domain_bound) + self.params_domain = params_bound diff --git a/pina/pinn.py b/pina/pinn.py new file mode 100644 index 0000000..ce3bee9 --- /dev/null +++ b/pina/pinn.py @@ -0,0 +1,488 @@ +from mpmath import chebyt, chop, taylor + +from .problem import Problem +import torch +import torch.nn as nn +import numpy as np +from .cube import Cube +from .segment import Segment +from .deep_feed_forward import DeepFeedForward +from pina.label_tensor import LabelTensor +torch.pi = torch.acos(torch.zeros(1)).item() * 2 # which is 3.1415927410125732 + +class PINN(object): + + def __init__(self, + problem, + model, + optimizer=torch.optim.Adam, + lr=0.001, + regularizer=0.00001, + data_weight=1., + dtype=torch.float64, + device='cpu', + lr_accelerate=None, + error_norm='mse'): + ''' + :param Problem problem: the formualation of the problem. + :param dict architecture: a dictionary containing the information to + build the model. Valid options are: + - inner_size [int] the number of neurons in the hidden layers; by + default is 20. + - n_layers [int] the number of hidden layers; by default is 4. + - func [nn.Module or str] the activation function; passing a `str` + is possible to chose adaptive function (between 'adapt_tanh'); by + default is non-adaptive iperbolic tangent. + :param float lr: the learning rate; default is 0.001 + :param float regularizer: the coefficient for L2 regularizer term + :param type dtype: the data type to use for the model. Valid option are + `torch.float32` and `torch.float64` (`torch.float16` only on GPU); + default is `torch.float64`. + :param float lr_accelete: the coefficient that controls the learning + rate increase, such that, for all the epoches in which the loss is + decreasing, the learning_rate is update using + $learning_rate = learning_rate * lr_accelerate$. + When the loss stops to decrease, the learning rate is set to the + initial value [TODO test parameters] + + ''' + + self.problem = problem + + + # self._architecture = architecture if architecture else dict() + # self._architecture['input_dimension'] = self.problem.domain_bound.shape[0] + # self._architecture['output_dimension'] = len(self.problem.variables) + # if hasattr(self.problem, 'params_domain'): + # self._architecture['input_dimension'] += self.problem.params_domain.shape[0] + + self.accelerate = lr_accelerate + + self.error_norm = error_norm + + if device == 'cuda' and not torch.cuda.is_available(): + raise RuntimeError + self.device = torch.device(device) + + self.dtype = dtype + self.history = [] + + self.model = model + self.model.to(dtype=self.dtype, device=self.device) + + self.input_pts = {} + self.truth_values = {} + + + self.trained_epoch = 0 + self.optimizer = optimizer( + self.model.parameters(), lr=lr, weight_decay=regularizer) + + self.data_weight = data_weight + + @property + def problem(self): + return self._problem + + @problem.setter + def problem(self, problem): + if not isinstance(problem, Problem): + raise TypeError + self._problem = problem + + def get_data_residuals(self): + + data_residuals = [] + + for output in self.data_pts: + data_values_pred = self.model(self.data_pts[output]) + data_residuals.append(data_values_pred - self.data_values[output]) + + return torch.cat(data_residuals) + + def get_phys_residuals(self): + """ + """ + + residuals = [] + for equation in self.problem.equation: + residuals.append(equation(self.phys_pts, self.model(self.phys_pts))) + return residuals + + + def _compute_norm(self, vec): + """ + Compute the norm of the `vec` one-dimensional tensor based on the + `self.error_norm` attribute. + + .. todo: complete + + :param vec torch.tensor: the tensor + """ + if isinstance(self.error_norm, int): + return torch.sum(torch.abs(vec**self.error_norm))**(1./self.error_norm) + elif self.error_norm == 'mse': + return torch.mean(vec**2) + elif self.error_norm == 'me': + return torch.mean(torch.abs(vec)) + else: + raise RuntimeError + + def save_state(self, filename): + + checkpoint = { + 'epoch': self.trained_epoch, + 'model_state': self.model.state_dict(), + 'optimizer_state' : self.optimizer.state_dict(), + 'optimizer_class' : self.optimizer.__class__, + 'history' : self.history, + } + + # TODO save also architecture param? + #if isinstance(self.model, DeepFeedForward): + # checkpoint['model_class'] = self.model.__class__ + # checkpoint['model_structure'] = { + # } + torch.save(checkpoint, filename) + + def load_state(self, filename): + + checkpoint = torch.load(filename) + self.model.load_state_dict(checkpoint['model_state']) + + + self.optimizer = checkpoint['optimizer_class'](self.model.parameters()) + self.optimizer.load_state_dict(checkpoint['optimizer_state']) + + self.trained_epoch = checkpoint['epoch'] + self.history = checkpoint['history'] + + print(self.history) + return self + + + def span_pts(self, n, mode='grid', locations='all'): + ''' + + ''' + + if locations == 'all': + locations = [condition for condition in self.problem.conditions] + + + for location in locations: + manifold, func = self.problem.conditions[location].values() + if torch.is_tensor(manifold): + pts = manifold + else: + pts = manifold.discretize(n, mode) + + pts = torch.from_numpy(pts) + + self.input_pts[location] = LabelTensor(pts, self.problem.input_variables) + self.input_pts[location].tensor.to(dtype=self.dtype, device=self.device) + self.input_pts[location].tensor.requires_grad_(True) + self.input_pts[location].tensor.retain_grad() + + + def plot_pts(self, locations='all'): + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + if locations == 'all': + locations = [condition for condition in self.problem.conditions] + + for location in locations: + x, y = self.input_pts[location].tensor.T + #plt.plot(x.detach(), y.detach(), 'o', label=location) + np.savetxt('burgers_{}_pts.txt'.format(location), self.input_pts[location].tensor.detach(), header='x y', delimiter=' ') + + gggg + + plt.legend() + plt.show() + + + + def train(self, stop=100, frequency_print=2, trial=None): + + epoch = 0 + + while True: + + losses = [] + for condition_name in self.problem.conditions: + pts = self.input_pts[condition_name] + predicted = self.model(pts.tensor) + if isinstance(self.problem.conditions[condition_name]['func'], list): + for func in self.problem.conditions[condition_name]['func']: + residuals = func(pts, predicted) + losses.append(self._compute_norm(residuals)) + else: + residuals = self.problem.conditions[condition_name]['func'](pts, predicted) + losses.append(self._compute_norm(residuals)) + #print(condition_name, losses[-1]) + + self.optimizer.zero_grad() + sum(losses).backward() + self.optimizer.step() + + self.trained_epoch += 1 + if epoch % 50 == 0: + self.history.append([loss.detach().item() for loss in losses]) + epoch += 1 + + if trial: + import optuna + trial.report(loss[0].item()+loss[1].item(), epoch) + if trial.should_prune(): + raise optuna.exceptions.TrialPruned() + + if isinstance(stop, int): + if epoch == stop: + break + elif isinstance(stop, float): + if sum(losses) < stop: + break + + if epoch % frequency_print == 0: + print('[epoch {:05d}] {:.6e} '.format(self.trained_epoch, sum(losses).item()), end='') + for loss in losses: + print('{:.6e} '.format(loss), end='') + print() + + return sum(losses).item() + + + def error(self, dtype='l2', res=100): + + import numpy as np + if hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + + elif hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + try: + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.dtype, device=self.device) + Z_pred = self.model(unrolled_pts) + Z_pred = Z_pred.detach().numpy().reshape(grids_container[0].shape) + + if dtype == 'l2': + return np.linalg.norm(Z_pred - Z_true)/np.linalg.norm(Z_true) + else: + # TODO H1 + pass + except: + print("") + print("Something went wrong...") + print("Not able to compute the error. Please pass a data solution or a true solution") + + def plot(self, res, filename=None, variable=None): + ''' + ''' + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + + pts_container = [] + #for mn, mx in [[-1, 1], [-1, 1]]: + for mn, mx in [[0, 1], [0, 1]]: + #for mn, mx in [[-1, 1], [0, 1]]: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + unrolled_pts.to(dtype=self.dtype) + Z_pred = self.model(unrolled_pts) + + ####################################################### + # poisson + # Z_truth = self.problem.truth_solution(unrolled_pts[:, 0], unrolled_pts[:, 1]) + # Z_pred = Z_pred.tensor.detach().reshape(grids_container[0].shape) + # Z_truth = Z_truth.detach().reshape(grids_container[0].shape) + # err = np.abs(Z_pred-Z_truth) + + + # with open('poisson2_nofeat_plot.txt', 'w') as f_: + # f_.write('x y truth pred e\n') + # for (x, y), tru, pre, e in zip(unrolled_pts, Z_truth.reshape(-1, 1), Z_pred.reshape(-1, 1), err.reshape(-1, 1)): + # f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + # n = Z_pred.shape[1] + # plt.figure(figsize=(16, 6)) + # plt.subplot(1, 3, 1) + # plt.contourf(*grids_container, Z_truth) + # plt.colorbar() + # plt.subplot(1, 3, 2) + # plt.contourf(*grids_container, Z_pred) + # plt.colorbar() + # plt.subplot(1, 3, 3) + # plt.contourf(*grids_container, err) + # plt.colorbar() + # plt.show() + + ####################################################### + # burgers + import scipy + data = scipy.io.loadmat('Data/burgers_shock.mat') + data_solution = {'grid': np.meshgrid(data['x'], data['t']), 'grid_solution': data['usol'].T} + + grids_container = data_solution['grid'] + print(data_solution['grid_solution'].shape) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + unrolled_pts.to(dtype=self.dtype) + Z_pred = self.model(unrolled_pts) + Z_truth = data_solution['grid_solution'] + + Z_pred = Z_pred.tensor.detach().reshape(grids_container[0].shape) + print(Z_pred, Z_truth) + err = np.abs(Z_pred.numpy()-Z_truth) + + + with open('burgers_nofeat_plot.txt', 'w') as f_: + f_.write('x y truth pred e\n') + for (x, y), tru, pre, e in zip(unrolled_pts, Z_truth.reshape(-1, 1), Z_pred.reshape(-1, 1), err.reshape(-1, 1)): + f_.write('{} {} {} {} {}\n'.format(x.item(), y.item(), tru.item(), pre.item(), e.item())) + n = Z_pred.shape[1] + plt.figure(figsize=(16, 6)) + plt.subplot(1, 3, 1) + plt.contourf(*grids_container, Z_truth,vmin=-1, vmax=1) + plt.colorbar() + plt.subplot(1, 3, 2) + plt.contourf(*grids_container, Z_pred, vmin=-1, vmax=1) + plt.colorbar() + plt.subplot(1, 3, 3) + plt.contourf(*grids_container, err) + plt.colorbar() + plt.show() + + + # for i, output in enumerate(Z_pred.tensor.T, start=1): + # output = output.detach().numpy().reshape(grids_container[0].shape) + # plt.subplot(1, n, i) + # plt.contourf(*grids_container, output) + # plt.colorbar() + + if filename is None: + plt.show() + else: + plt.savefig(filename) + + def plot_params(self, res, param, filename=None, variable=None): + ''' + ''' + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + if hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + n_plot = 2 + elif hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + n_plot = 2 + else: + n_plot = 1 + + fig, axs = plt.subplots(nrows=1, ncols=n_plot, figsize=(n_plot*6,4)) + if not isinstance(axs, np.ndarray): axs = [axs] + + if hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + elif hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.type) + #print(unrolled_pts) + #print(param) + param_unrolled_pts = torch.cat((unrolled_pts, param.repeat(unrolled_pts.shape[0], 1)), 1) + if variable==None: + variable = self.problem.variables[0] + Z_pred = self.evaluate(param_unrolled_pts)[variable] + variable = "Solution" + else: + Z_pred = self.evaluate(param_unrolled_pts)[variable] + + Z_pred= Z_pred.detach().numpy().reshape(grids_container[0].shape) + set_pred = axs[0].contourf(*grids_container, Z_pred) + axs[0].set_title('PINN [trained epoch = {}]'.format(self.trained_epoch) + " " + variable) #TODO add info about parameter in the title + fig.colorbar(set_pred, ax=axs[0]) + + if n_plot == 2: + + set_true = axs[1].contourf(*grids_container, Z_true) + + axs[1].set_title('Truth solution') + fig.colorbar(set_true, ax=axs[1]) + + if filename is None: + plt.show() + else: + fig.savefig(filename + " " + variable) + + def plot_error(self, res, filename=None): + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + + fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(6,4)) + if not isinstance(axs, np.ndarray): axs = [axs] + + if hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + elif hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + try: + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.type) + + Z_pred = self.model(unrolled_pts) + Z_pred = Z_pred.detach().numpy().reshape(grids_container[0].shape) + set_pred = axs[0].contourf(*grids_container, abs(Z_pred - Z_true)) + axs[0].set_title('PINN [trained epoch = {}]'.format(self.trained_epoch) + "Pointwise Error") + fig.colorbar(set_pred, ax=axs[0]) + + if filename is None: + plt.show() + else: + fig.savefig(filename) + except: + print("") + print("Something went wrong...") + print("Not able to plot the error. Please pass a data solution or a true solution") + +''' +print(self.pred_loss.item(),loss.item(), self.old_loss.item()) +if self.accelerate is not None: + if self.pred_loss > loss and loss >= self.old_loss: + self.current_lr = self.original_lr + #print('restart') + elif (loss-self.pred_loss).item() < 0.1: + self.current_lr += .5*self.current_lr + #print('powa') + else: + self.current_lr -= .5*self.current_lr + #print(self.current_lr) + #self.current_lr = min(loss.item()*3, 0.02) + + for g in self.optimizer.param_groups: + g['lr'] = self.current_lr +''' diff --git a/pina/ppinn.py b/pina/ppinn.py new file mode 100644 index 0000000..f2eafe8 --- /dev/null +++ b/pina/ppinn.py @@ -0,0 +1,465 @@ +from mpmath import chebyt, chop, taylor + +from .problem import Problem +import torch +import torch.nn as nn +import numpy as np +from .cube import Cube +from .deep_feed_forward import DeepFeedForward +from pina.label_tensor import LabelTensor +from pina.pinn import PINN +torch.pi = torch.acos(torch.zeros(1)).item() * 2 # which is 3.1415927410125732 + +class ParametricPINN(PINN): + + def __init__(self, + problem, + model, + optimizer=torch.optim.Adam, + lr=0.001, + regularizer=0.00001, + data_weight=1., + dtype=torch.float64, + device='cpu', + lr_accelerate=None, + error_norm='mse'): + ''' + :param Problem problem: the formualation of the problem. + :param dict architecture: a dictionary containing the information to + build the model. Valid options are: + - inner_size [int] the number of neurons in the hidden layers; by + default is 20. + - n_layers [int] the number of hidden layers; by default is 4. + - func [nn.Module or str] the activation function; passing a `str` + is possible to chose adaptive function (between 'adapt_tanh'); by + default is non-adaptive iperbolic tangent. + :param float lr: the learning rate; default is 0.001 + :param float regularizer: the coefficient for L2 regularizer term + :param type dtype: the data type to use for the model. Valid option are + `torch.float32` and `torch.float64` (`torch.float16` only on GPU); + default is `torch.float64`. + :param float lr_accelete: the coefficient that controls the learning + rate increase, such that, for all the epoches in which the loss is + decreasing, the learning_rate is update using + $learning_rate = learning_rate * lr_accelerate$. + When the loss stops to decrease, the learning rate is set to the + initial value [TODO test parameters] + + ''' + + self.problem = problem + + + # self._architecture = architecture if architecture else dict() + # self._architecture['input_dimension'] = self.problem.domain_bound.shape[0] + # self._architecture['output_dimension'] = len(self.problem.variables) + # if hasattr(self.problem, 'params_domain'): + # self._architecture['input_dimension'] += self.problem.params_domain.shape[0] + + self.accelerate = lr_accelerate + + self.error_norm = error_norm + + if device == 'cuda' and not torch.cuda.is_available(): + raise RuntimeError + self.device = torch.device(device) + + self.dtype = dtype + self.history = [] + + self.model = model + self.model.to(dtype=self.dtype, device=self.device) + + self.input_pts = {} + self.truth_values = {} + + + self.trained_epoch = 0 + self.optimizer = optimizer( + self.model.parameters(), lr=lr, weight_decay=regularizer) + + self.data_weight = data_weight + + @property + def problem(self): + return self._problem + + @problem.setter + def problem(self, problem): + if not isinstance(problem, Problem): + raise TypeError + self._problem = problem + + def get_data_residuals(self): + + data_residuals = [] + + for output in self.data_pts: + data_values_pred = self.model(self.data_pts[output]) + data_residuals.append(data_values_pred - self.data_values[output]) + + return torch.cat(data_residuals) + + def get_phys_residuals(self): + """ + """ + + residuals = [] + for equation in self.problem.equation: + residuals.append(equation(self.phys_pts, self.model(self.phys_pts))) + return residuals + + + def _compute_norm(self, vec): + """ + Compute the norm of the `vec` one-dimensional tensor based on the + `self.error_norm` attribute. + + .. todo: complete + + :param vec torch.tensor: the tensor + """ + if isinstance(self.error_norm, int): + return torch.sum(torch.abs(vec**self.error_norm))**(1./self.error_norm) + elif self.error_norm == 'mse': + return torch.mean(vec**2) + elif self.error_norm == 'me': + return torch.mean(torch.abs(vec)) + else: + raise RuntimeError + + def save_state(self, filename): + + checkpoint = { + 'epoch': self.trained_epoch, + 'model_state': self.model.state_dict(), + 'optimizer_state' : self.optimizer.state_dict(), + 'optimizer_class' : self.optimizer.__class__, + } + + # TODO save also architecture param? + #if isinstance(self.model, DeepFeedForward): + # checkpoint['model_class'] = self.model.__class__ + # checkpoint['model_structure'] = { + # } + torch.save(checkpoint, filename) + + def load_state(self, filename): + + checkpoint = torch.load(filename) + self.model.load_state_dict(checkpoint['model_state']) + + self.optimizer = checkpoint['optimizer_class'](self.model.parameters()) + self.optimizer.load_state_dict(checkpoint['optimizer_state']) + + self.trained_epoch = checkpoint['epoch'] + return self + + + def span_pts(self, n, mode='grid', locations='all'): + ''' + + ''' + + if locations == 'all': + locations = [condition for condition in self.problem.conditions] + + + for location in locations: + manifold, func = self.problem.conditions[location].values() + if torch.is_tensor(manifold): + pts = manifold + else: + pts = manifold.discretize(n, mode) + + pts = torch.from_numpy(pts) + + self.input_pts[location] = LabelTensor(pts, self.problem.input_variables) + self.input_pts[location].tensor.to(dtype=self.dtype, device=self.device) + self.input_pts[location].tensor.requires_grad_(True) + self.input_pts[location].tensor.retain_grad() + + + def train(self, stop=100, frequency_print=2, trial=None): + + epoch = 0 + + ## TODO for elliptic + # parameters = torch.cat(torch.linspace( + # self.problem.parameter_domain[0, 0], + # self.problem.parameter_domain[0, 1], + # 5) + + ## for param laplacian + #parameters = torch.rand(50, 2)*2-1 + parameters = torch.from_numpy(Cube([[-1, 1], [-1, 1]]).discretize(5, 'grid')) + # alpha_p = torch.logspace(start=-2, end=0, steps=10) + # mu_p = torch.linspace(0.5, 3, 5) + # g1_, g2_ = torch.meshgrid(alpha_p, mu_p) + # parameters = torch.cat([g2_.reshape(-1, 1), g1_.reshape(-1, 1)], axis=1) + print(parameters) + + while True: + + losses = [] + + for condition_name in self.problem.conditions: + + pts = self.input_pts[condition_name] + pts = torch.cat([ + pts.tensor.repeat_interleave(parameters.shape[0], dim=0), + torch.tile(parameters, (pts.tensor.shape[0], 1)) + ], axis=1) + pts = LabelTensor(pts, self.problem.input_variables + self.problem.parameters) + predicted = self.model(pts.tensor) + #predicted = self.model(pts) + + if isinstance(self.problem.conditions[condition_name]['func'], list): + for func in self.problem.conditions[condition_name]['func']: + residuals = func(pts, None, predicted) + losses.append(self._compute_norm(residuals)) + else: + residuals = self.problem.conditions[condition_name]['func'](pts, None, predicted) + losses.append(self._compute_norm(residuals)) + self.optimizer.zero_grad() + sum(losses).backward() + self.optimizer.step() + + #for p in parameters: + # pts = self.input_pts[condition_name] + # #pts = torch.cat([pts.tensor, p.double().repeat(pts.tensor.shape[0]).reshape(-1, 2)], axis=1) + # #pts = torch.cat([pts.tensor, p.double().repeat(pts.tensor.shape[0]).reshape(-1, 1)], axis=1) + # #print(self.problem.input_variables) + # # print(self.problem.parameters) + # # print(pts.shape) + # print(pts.tensor.repeat_interleave(parameters.shape[0])) + # # print(pts) + # # gg + # a = torch.cat([ + # pts.tensor.repeat_interleave(parameters.shape[0], dim=0), + # torch.tile(parameters, (pts.tensor.shape[0], 1)) + # ], axis=1) + # for i in a: + # print(i.detach()) + # ttt + # pts = LabelTensor(pts, self.problem.input_variables + self.problem.parameters) + + + # ffff + # print(pts.labels) + + # predicted = self.model(pts.tensor) + # #predicted = self.model(pts) + # if isinstance(self.problem.conditions[condition_name]['func'], list): + # for func in self.problem.conditions[condition_name]['func']: + # residuals = func(pts, LabelTensor(p.reshape(1, -1), ['mu', 'alpha']), predicted) + # tmp_losses.append(self._compute_norm(residuals)) + # else: + # residuals = self.problem.conditions[condition_name]['func'](pts, p, predicted) + # tmp_losses.append(self._compute_norm(residuals)) + #losses.append(sum(tmp_losses)) + + + self.trained_epoch += 1 + #if epoch % 10 == 0: + # self.history.append(losses) + + epoch += 1 + + if trial: + import optuna + rial.report(loss[0].item()+loss[1].item(), epoch) + if trial.should_prune(): + raise optuna.exceptions.TrialPruned() + + if isinstance(stop, int): + if epoch == stop: + break + elif isinstance(stop, float): + if loss[0].item() + loss[1].item() < stop: + break + + if epoch % frequency_print == 0: + print('[epoch {:05d}] {:.6e} '.format(self.trained_epoch, sum(losses).item()), end='') + for loss in losses: + print('{:.6e} '.format(loss), end='') + print() + + return sum(losses).item() + + + def error(self, dtype='l2', res=100): + + import numpy as np + if hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + + elif hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + try: + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.dtype, device=self.device) + Z_pred = self.model(unrolled_pts) + Z_pred = Z_pred.detach().numpy().reshape(grids_container[0].shape) + + if dtype == 'l2': + return np.linalg.norm(Z_pred - Z_true)/np.linalg.norm(Z_true) + else: + # TODO H1 + pass + except: + print("") + print("Something went wrong...") + print("Not able to compute the error. Please pass a data solution or a true solution") + + def plot(self, res, param, filename=None, variable=None): + ''' + ''' + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + pts_container = [] + for mn, mx in [[-1, 1], [-1, 1]]: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T + unrolled_pts = torch.cat([unrolled_pts, param.double().repeat(unrolled_pts.shape[0]).reshape(-1, 2)], axis=1) + + unrolled_pts.to(dtype=self.dtype) + unrolled_pts = LabelTensor(unrolled_pts, ['x1', 'x2', 'mu1', 'mu2']) + + Z_pred = self.model(unrolled_pts.tensor) + n = Z_pred.tensor.shape[1] + plt.figure(figsize=(6*n, 6)) + + for i, output in enumerate(Z_pred.tensor.T, start=1): + output = output.detach().numpy().reshape(grids_container[0].shape) + plt.subplot(1, n, i) + plt.contourf(*grids_container, output) + plt.colorbar() + + if filename is None: + plt.show() + else: + plt.savefig(filename) + + def plot_params(self, res, param, filename=None, variable=None): + ''' + ''' + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + if hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + n_plot = 2 + elif hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + n_plot = 2 + else: + n_plot = 1 + + fig, axs = plt.subplots(nrows=1, ncols=n_plot, figsize=(n_plot*6,4)) + if not isinstance(axs, np.ndarray): axs = [axs] + + if hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + elif hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + grids_container = np.meshgrid(*pts_container) + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.type) + #print(unrolled_pts) + #print(param) + param_unrolled_pts = torch.cat((unrolled_pts, param.repeat(unrolled_pts.shape[0], 1)), 1) + if variable==None: + variable = self.problem.variables[0] + Z_pred = self.evaluate(param_unrolled_pts)[variable] + variable = "Solution" + else: + Z_pred = self.evaluate(param_unrolled_pts)[variable] + + Z_pred= Z_pred.detach().numpy().reshape(grids_container[0].shape) + set_pred = axs[0].contourf(*grids_container, Z_pred) + axs[0].set_title('PINN [trained epoch = {}]'.format(self.trained_epoch) + " " + variable) #TODO add info about parameter in the title + fig.colorbar(set_pred, ax=axs[0]) + + if n_plot == 2: + + set_true = axs[1].contourf(*grids_container, Z_true) + + axs[1].set_title('Truth solution') + fig.colorbar(set_true, ax=axs[1]) + + if filename is None: + plt.show() + else: + fig.savefig(filename + " " + variable) + + def plot_error(self, res, filename=None): + import matplotlib + matplotlib.use('GTK3Agg') + import matplotlib.pyplot as plt + + + fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(6,4)) + if not isinstance(axs, np.ndarray): axs = [axs] + + if hasattr(self.problem, 'data_solution') and self.problem.data_solution is not None: + grids_container = self.problem.data_solution['grid'] + Z_true = self.problem.data_solution['grid_solution'] + elif hasattr(self.problem, 'truth_solution') and self.problem.truth_solution is not None: + pts_container = [] + for mn, mx in self.problem.domain_bound: + pts_container.append(np.linspace(mn, mx, res)) + + grids_container = np.meshgrid(*pts_container) + Z_true = self.problem.truth_solution(*grids_container) + try: + unrolled_pts = torch.tensor([t.flatten() for t in grids_container]).T.to(dtype=self.type) + + Z_pred = self.model(unrolled_pts) + Z_pred = Z_pred.detach().numpy().reshape(grids_container[0].shape) + set_pred = axs[0].contourf(*grids_container, abs(Z_pred - Z_true)) + axs[0].set_title('PINN [trained epoch = {}]'.format(self.trained_epoch) + "Pointwise Error") + fig.colorbar(set_pred, ax=axs[0]) + + if filename is None: + plt.show() + else: + fig.savefig(filename) + except: + print("") + print("Something went wrong...") + print("Not able to plot the error. Please pass a data solution or a true solution") + +''' +print(self.pred_loss.item(),loss.item(), self.old_loss.item()) +if self.accelerate is not None: + if self.pred_loss > loss and loss >= self.old_loss: + self.current_lr = self.original_lr + #print('restart') + elif (loss-self.pred_loss).item() < 0.1: + self.current_lr += .5*self.current_lr + #print('powa') + else: + self.current_lr -= .5*self.current_lr + #print(self.current_lr) + #self.current_lr = min(loss.item()*3, 0.02) + + for g in self.optimizer.param_groups: + g['lr'] = self.current_lr +''' diff --git a/pina/problem.py b/pina/problem.py new file mode 100644 index 0000000..05bb66e --- /dev/null +++ b/pina/problem.py @@ -0,0 +1,49 @@ +import torch + +class Problem(object): + + def __init__(self, *args, **kwargs): + raise NotImplemented + + @property + def variables(self): + return self._variables + + @variables.setter + def variables(self, variables): + if variables is None: + variables = ['var'] + self._variables = variables + + @property + def boundary_conditions(self): + return self._bc + + @boundary_conditions.setter + def boundary_conditions(self, bc): + if isinstance(bc, (list, tuple)): + bc = {'var': bc} + self._bc = bc + + @property + def spatial_dimensions(self): + return self._spatial_dimensions + + @staticmethod + def grad(output_, input_): + gradients = torch.autograd.grad( + output_, + input_.tensor, + grad_outputs=torch.ones(output_.size()).to( + dtype=input_.tensor.dtype, + device=input_.tensor.device), + create_graph=True, retain_graph=True, allow_unused=True)[0] + from pina.label_tensor import LabelTensor + return LabelTensor(gradients, input_.labels) + + + + def __str__(self): + s = '' + #s = 'Variables: {}\n'.format(self.variables) + return s diff --git a/pina/problem1d.py b/pina/problem1d.py new file mode 100644 index 0000000..8df84cd --- /dev/null +++ b/pina/problem1d.py @@ -0,0 +1,11 @@ + +from .problem import Problem +import numpy as np + +class Problem1D(Problem): + + def __init__(self, variables=None, bc=None): + self._spatial_dimensions = 1 + self.variables = variables + print(bc) + self.bc = bc diff --git a/pina/problem2d.py b/pina/problem2d.py new file mode 100644 index 0000000..f54e0b0 --- /dev/null +++ b/pina/problem2d.py @@ -0,0 +1,16 @@ +from .problem import Problem + + +class Problem2D(Problem): + + spatial_dimensions = 2 + + @property + def boundary_condition(self): + return self._boundary_condition + + @boundary_condition.setter + def boundary_condition(self, bc): + self._boundary_condition = bc + + diff --git a/pina/segment.py b/pina/segment.py new file mode 100644 index 0000000..5414ba9 --- /dev/null +++ b/pina/segment.py @@ -0,0 +1,27 @@ +import torch +import numpy as np +from .chebyshev import chebyshev_roots + +class Segment(): + def __init__(self, p1, p2): + self.p1 = p1 + self.p2 = p2 + + def discretize(self, n, mode='random'): + pts = [] + + if mode == 'random': + iterator = np.random.uniform(0, 1, n) + elif mode == 'grid': + iterator = np.linspace(0, 1, n) + elif mode == 'chebyshev': + iterator = chebyshev_roots(n) * .5 + .5 + + for k in iterator: + x = self.p1[0] + k*(self.p2[0]-self.p1[0]) + y = self.p1[1] + k*(self.p2[1]-self.p1[1]) + pts.append((x, y)) + + pts = np.array(pts) + return pts + diff --git a/pina/tdproblem1d.py b/pina/tdproblem1d.py new file mode 100644 index 0000000..4df4431 --- /dev/null +++ b/pina/tdproblem1d.py @@ -0,0 +1,27 @@ +import numpy as np +from .problem1d import Problem1D +from .segment import Segment + +class TimeDepProblem1D(Problem1D): + + def __init__(self, variables=None, bc=None, initial=None, tend=None, domain_bound=None): + self.variables = variables + self._spatial_dimensions = 2 + self.tend = tend + self.tstart = 0 + if domain_bound is None: + bound_pts = [bc[0] for bc in self.boundary_conditions] + domain_bound = np.array([ + [min(bound_pts), max(bound_pts)], + [self.tstart, self.tend ] + ]) + + self.domain_bound = np.array([[-1, 1],[0, 1]])#domain_bound + print(domain_bound) + self.boundary_conditions = ( + (Segment((bc[0][0], self.tstart), (bc[1][0], self.tstart)), initial), + (Segment((bc[0][0], self.tstart), (bc[0][0], self.tend)), bc[0][1]), + (Segment((bc[1][0], self.tstart), (bc[1][0], self.tend)), bc[1][1]) + ) + + diff --git a/problems/burgers.py b/problems/burgers.py new file mode 100644 index 0000000..3d90ce5 --- /dev/null +++ b/problems/burgers.py @@ -0,0 +1,78 @@ +import numpy as np +import scipy.io +import torch + +from pina.problem import Problem +from pina.segment import Segment +from pina.cube import Cube +from pina.tdproblem1d import TimeDepProblem1D + +def tmp_grad(output_, input_): + return torch.autograd.grad( + output_, + input_.tensor, + grad_outputs=torch.ones(output_.size()).to( + dtype=input_.tensor.dtype, + device=input_.tensor.device), + create_graph=True, retain_graph=True, allow_unused=True)[0] + +class Burgers1D(TimeDepProblem1D): + + def __init__(self): + + + def burger_equation(input_, output_): + + grad_u = self.grad(output_['u'], input_) + grad_x, grad_t = tmp_grad(output_['u'], input_).T + gradgrad_u_x = self.grad(grad_u['x'], input_) + grad_xx = tmp_grad(grad_x, input_)[:, 0] + #print(grad_t, grad_u['t']) + + #rrrr + return grad_u['t'] + output_['u']*grad_u['x'] - (0.01/torch.pi)*gradgrad_u_x['x'] + + + def nil_dirichlet(input_, output_): + u_expected = 0.0 + return output_['u'] - u_expected + + def initial_condition(input_, output_): + u_expected = -torch.sin(torch.pi*input_['x']) + return output_['u'] - u_expected + + + + self.conditions = { + 'gamma1': {'location': Segment((-1, 0), (-1, 1)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment(( 1, 0), ( 1, 1)), 'func': nil_dirichlet}, + 'initia': {'location': Segment((-1, 0), ( 1, 0)), 'func': initial_condition}, + 'D': {'location': Cube([[-1, 1],[0,1]]), 'func': burger_equation} + } + + self.input_variables = ['x', 't'] + self.output_variables = ['u'] + self.spatial_domain = Cube([[0, 1]]) + self.temporal_domain = Cube([[0, 1]]) + +bc = ( + (-1, lambda x: torch.zeros(x.shape[0], 1)), + ( 1, lambda x: torch.zeros(x.shape[0], 1)) +) + +initial = lambda x: -np.sin(np.pi*x[:,0]).reshape(-1, 1) + +def equation(x, fx): + grad_x, grad_t = Problem.grad(fx, x).T + grad_xx = Problem.grad(grad_x, x)[:, 0] + a = grad_t + fx.flatten()*grad_x - (0.01/torch.pi)*grad_xx + return a + + +burgers = TimeDepProblem1D(bc=bc, initial=initial, tend=1, domain_bound=[-1, 1]) +burgers.equation = equation + +# read data for errors and plots +data = scipy.io.loadmat('Data/burgers_shock.mat') +data_solution = {'grid': np.meshgrid(data['x'], data['t']), 'grid_solution': data['usol'].T} +burgers.data_solution = data_solution diff --git a/problems/elliptic_optimal_control.py b/problems/elliptic_optimal_control.py new file mode 100644 index 0000000..adef451 --- /dev/null +++ b/problems/elliptic_optimal_control.py @@ -0,0 +1,49 @@ +import numpy as np +import torch +from pina.problem import Problem +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D + +xmin, xmax, ymin, ymax = -1, 1, -1, 1 + +class EllipticOptimalControl(Problem2D): + + def __init__(self, alpha=1): + + def term1(input_, output_): + grad_p = self.grad(output_['p'], input_) + gradgrad_p_x1 = self.grad(grad_p['x1'], input_) + gradgrad_p_x2 = self.grad(grad_p['x2'], input_) + yd = 2.0 + return output_['y'] - yd - (gradgrad_p_x1['x1'] + gradgrad_p_x2['x2']) + + def term2(input_, output_): + grad_y = self.grad(output_['y'], input_) + gradgrad_y_x1 = self.grad(grad_y['x1'], input_) + gradgrad_y_x2 = self.grad(grad_y['x2'], input_) + return - (gradgrad_y_x1['x1'] + gradgrad_y_x2['x2']) - output_['u'] + + def term3(input_, output_): + return output_['p'] - output_['u']*alpha + + + def nil_dirichlet(input_, output_): + y_value = 0.0 + p_value = 0.0 + return torch.abs(output_['y'] - y_value) + torch.abs(output_['p'] - p_value) + + self.conditions = { + 'gamma1': {'location': Segment((xmin, ymin), (xmax, ymin)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment((xmax, ymin), (xmax, ymax)), 'func': nil_dirichlet}, + 'gamma3': {'location': Segment((xmax, ymax), (xmin, ymax)), 'func': nil_dirichlet}, + 'gamma4': {'location': Segment((xmin, ymax), (xmin, ymin)), 'func': nil_dirichlet}, + 'D1': {'location': Cube([[xmin, xmax], [ymin, ymax]]), 'func': [term1, term2, term3]}, + #'D2': {'location': Cube([[0, 1], [0, 1]]), 'func': term2}, + #'D3': {'location': Cube([[0, 1], [0, 1]]), 'func': term3} + } + + self.input_variables = ['x1', 'x2'] + self.output_variables = ['u', 'p', 'y'] + self.spatial_domain = Cube([[xmin, xmax], [xmin, xmax]]) + diff --git a/problems/laplacian_optimal_control_parametric.py b/problems/laplacian_optimal_control_parametric.py new file mode 100644 index 0000000..b6e01ea --- /dev/null +++ b/problems/laplacian_optimal_control_parametric.py @@ -0,0 +1,55 @@ +import numpy as np +import torch +from pina.problem import Problem +from pina.segment import Segment +from pina.parametricproblem2d import ParametricProblem2D + +bc = { + 'y': ( + (Segment((0, 0), (4, 0)), lambda x: torch.ones(x.shape[0], 1)), + (Segment((4, 0), (4, 1)), lambda x: torch.ones(x.shape[0], 1)), + (Segment((4, 1), (0, 1)), lambda x: torch.ones(x.shape[0], 1)), + (Segment((0, 1), (0, 0)), lambda x: torch.ones(x.shape[0], 1)), + ), + 'p': ( + (Segment((0, 0), (4, 0)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment((4, 0), (4, 1)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment((4, 1), (0, 1)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment((0, 1), (0, 0)), lambda x: torch.zeros(x.shape[0], 1)), + ) +} + +# optimal control parameters and data +alpha = 1e-5 +# yd = 10*x[:, 0]*(1-x[:, 0])*x[:, 1]*(1-x[:, 1]) +# three variables +# state y = f[0] +# control u = f[1] +# adjoint p = f[2] + +# the three variables + +def adjoint_eq(x, f): + grad_x, grad_y = Problem.grad(f['p'], x)[:, :2].T + grad_xx = Problem.grad(grad_x, x)[:, 0] + grad_yy = Problem.grad(grad_y, x)[:, 1] + return - grad_xx - grad_yy - f['y'] + 1*(x[:, 0] <= 1) + x[:, 2]*(x[:, 0] > 1) + +def control_eq(x, f): + return alpha*f['u'] - f['p'] + +def state_eq(x, f): + grad_x, grad_y = Problem.grad(f['y'], x)[:, :2].T + grad_xx = Problem.grad(grad_x, x)[:, 0] + grad_yy = Problem.grad(grad_y, x)[:, 1] + return - grad_xx - grad_yy - f['u'] + +def equation(x, f): + return state_eq(x, f) + control_eq(x, f) + adjoint_eq(x, f) + +laplace = ParametricProblem2D( + variables=['y', 'u', 'p'], + bc=bc, + domain_bound=np.array([[0, 4],[0, 1]]), + params_bound=np.array([[0.5, 2.5]])) +laplace.equation = equation diff --git a/problems/laplacian_parametric.py b/problems/laplacian_parametric.py new file mode 100644 index 0000000..0cf64e6 --- /dev/null +++ b/problems/laplacian_parametric.py @@ -0,0 +1,29 @@ +import numpy as np +import torch +from pina.problem import Problem +from pina.segment import Segment +from pina.parametricproblem2d import ParametricProblem2D + +bc = ( + (Segment((-1, -1), ( 1, -1)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment(( 1, -1), ( 1, 1)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment(( 1, 1), (-1, 1)), lambda x: torch.zeros(x.shape[0], 1)), + (Segment((-1, 1), (-1, -1)), lambda x: torch.zeros(x.shape[0], 1)), +) + +params_domain = np.array([ + [-1.0, 1.0], + [-1.0, 1.0]]) + +def equation(x, fx): + grad_x, grad_y = Problem.grad(fx, x)[:, :2].T + grad_xx = Problem.grad(grad_x, x)[:, 0] + grad_yy = Problem.grad(grad_y, x)[:, 1] + a = grad_xx + grad_yy - torch.exp(- 2*(x[:, 0] - x[:, 2])**2 - 2*(x[:, 1] - x[:, 3])**2) + return a + + +laplace = ParametricProblem2D(bc=bc, domain_bound=params_domain, params_bound=params_domain) + +laplace.equation = equation + diff --git a/problems/parametric_elliptic_optimal_control.py b/problems/parametric_elliptic_optimal_control.py new file mode 100644 index 0000000..22c9940 --- /dev/null +++ b/problems/parametric_elliptic_optimal_control.py @@ -0,0 +1,53 @@ +import numpy as np +import torch +from pina.problem import Problem +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D + +xmin, xmax, ymin, ymax = -1, 1, -1, 1 + +class ParametricEllipticOptimalControl(Problem2D): + + def __init__(self, alpha=1): + + def term1(input_, param_, output_): + grad_p = self.grad(output_['p'], input_) + gradgrad_p_x1 = self.grad(grad_p['x1'], input_) + gradgrad_p_x2 = self.grad(grad_p['x2'], input_) + return output_['y'] - param_ - (gradgrad_p_x1['x1'] + gradgrad_p_x2['x2']) + + def term2(input_, param_, output_): + grad_y = self.grad(output_['y'], input_) + gradgrad_y_x1 = self.grad(grad_y['x1'], input_) + gradgrad_y_x2 = self.grad(grad_y['x2'], input_) + return - (gradgrad_y_x1['x1'] + gradgrad_y_x2['x2']) - output_['u_param'] + + def term3(input_, param_, output_): + return output_['p'] - output_['u_param']*alpha + + + def term(input_, param_, output_): + return term1( input_, param_, output_) +term2( input_, param_, output_) + term3( input_, param_, output_) + + def nil_dirichlet(input_, param_, output_): + y_value = 0.0 + p_value = 0.0 + return torch.abs(output_['y'] - y_value) + torch.abs(output_['p'] - p_value) + + self.conditions = { + 'gamma1': {'location': Segment((xmin, ymin), (xmax, ymin)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment((xmax, ymin), (xmax, ymax)), 'func': nil_dirichlet}, + 'gamma3': {'location': Segment((xmax, ymax), (xmin, ymax)), 'func': nil_dirichlet}, + 'gamma4': {'location': Segment((xmin, ymax), (xmin, ymin)), 'func': nil_dirichlet}, + 'D1': {'location': Cube([[xmin, xmax], [ymin, ymax]]), 'func': term}, + #'D2': {'location': Cube([[0, 1], [0, 1]]), 'func': term2}, + #'D3': {'location': Cube([[0, 1], [0, 1]]), 'func': term3} + } + + self.input_variables = ['x1', 'x2'] + self.output_variables = ['u', 'p', 'y'] + self.parameters = ['mu'] + self.spatial_domain = Cube([[xmin, xmax], [xmin, xmax]]) + self.parameter_domain = np.array([[0.5, 3]]) + diff --git a/problems/parametric_elliptic_optimal_control_alpha_variable.py b/problems/parametric_elliptic_optimal_control_alpha_variable.py new file mode 100644 index 0000000..dd12ebe --- /dev/null +++ b/problems/parametric_elliptic_optimal_control_alpha_variable.py @@ -0,0 +1,52 @@ +import numpy as np +import torch +from pina.problem import Problem +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D + +xmin, xmax, ymin, ymax = -1, 1, -1, 1 + +class ParametricEllipticOptimalControl(Problem2D): + + def __init__(self, alpha=1): + + def term1(input_, param_, output_): + grad_p = self.grad(output_['p'], input_) + gradgrad_p_x1 = self.grad(grad_p['x1'], input_) + gradgrad_p_x2 = self.grad(grad_p['x2'], input_) + #print('mu', input_['mu']) + return output_['y'] - input_['mu'] - (gradgrad_p_x1['x1'] + gradgrad_p_x2['x2']) + + def term2(input_, param_, output_): + grad_y = self.grad(output_['y'], input_) + gradgrad_y_x1 = self.grad(grad_y['x1'], input_) + gradgrad_y_x2 = self.grad(grad_y['x2'], input_) + return - (gradgrad_y_x1['x1'] + gradgrad_y_x2['x2']) - output_['u_param'] + + def term3(input_, param_, output_): + #print('a', input_['alpha'], output_['p'], output_['u_param']) + return output_['p'] - output_['u_param']*input_['alpha'] + + + def nil_dirichlet(input_, param_, output_): + y_value = 0.0 + p_value = 0.0 + return torch.abs(output_['y'] - y_value) + torch.abs(output_['p'] - p_value) + + self.conditions = { + 'gamma1': {'location': Segment((xmin, ymin), (xmax, ymin)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment((xmax, ymin), (xmax, ymax)), 'func': nil_dirichlet}, + 'gamma3': {'location': Segment((xmax, ymax), (xmin, ymax)), 'func': nil_dirichlet}, + 'gamma4': {'location': Segment((xmin, ymax), (xmin, ymin)), 'func': nil_dirichlet}, + 'D1': {'location': Cube([[xmin, xmax], [ymin, ymax]]), 'func': [term1, term2]}, + #'D2': {'location': Cube([[0, 1], [0, 1]]), 'func': term2}, + #'D3': {'location': Cube([[0, 1], [0, 1]]), 'func': term3} + } + + self.input_variables = ['x1', 'x2'] + self.output_variables = ['u', 'p', 'y'] + self.parameters = ['mu', 'alpha'] + self.spatial_domain = Cube([[xmin, xmax], [xmin, xmax]]) + self.parameter_domain = np.array([[0.5, 3], [0.0001, 1]]) + diff --git a/problems/parametric_poisson.py b/problems/parametric_poisson.py new file mode 100644 index 0000000..c42d339 --- /dev/null +++ b/problems/parametric_poisson.py @@ -0,0 +1,43 @@ +import numpy as np +import torch +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D +from pina.problem import Problem + + +class ParametricPoisson2DProblem(Problem2D): + + def __init__(self): + + def laplace_equation(input_, param_, output_): + grad_u = self.grad(output_['u'], input_) + gradgrad_u_x = self.grad(grad_u['x'], input_) + gradgrad_u_y = self.grad(grad_u['y'], input_) + force_term = torch.exp( + - 2*(input_['x'] - input_['mu1'])**2 - + 2*(input_['y'] - input_['mu2'])**2 + ) + return gradgrad_u_x['x'] + gradgrad_u_y['y'] - force_term + + def nil_dirichlet(input_, param_, output_): + value = 0.0 + return output_['u'] - value + + self.conditions = { + 'gamma1': {'location': Segment((-1, -1), ( 1, -1)),'func': nil_dirichlet}, + 'gamma2': {'location': Segment(( 1, -1), ( 1, 1)),'func': nil_dirichlet}, + 'gamma3': {'location': Segment(( 1, 1), (-1, 1)),'func': nil_dirichlet}, + 'gamma4': {'location': Segment((-1, 1), (-1, -1)),'func': nil_dirichlet}, + 'D': {'location': Cube([[-1, 1], [-1, 1]]), 'func': laplace_equation} + } + + self.input_variables = ['x', 'y'] + self.output_variables = ['u'] + self.parameters = ['mu1', 'mu2'] + #self.truth_solution = poisson_sol + self.spatial_domain = Cube([[0, 1], [0, 1]]) + self.parameter_domain = np.array([[-1, 1], [-1, 1]]) + + + #self.check() # Check the problem is correctly set diff --git a/problems/poisson2D.py b/problems/poisson2D.py new file mode 100644 index 0000000..ea299e4 --- /dev/null +++ b/problems/poisson2D.py @@ -0,0 +1,41 @@ +import numpy as np +import torch +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D +from pina.problem import Problem + + +class Poisson2DProblem(Problem2D): + + def __init__(self): + + def laplace_equation(input_, output_): + grad_u = self.grad(output_['u'], input_) + gradgrad_u_x = self.grad(grad_u['x'], input_) + gradgrad_u_y = self.grad(grad_u['y'], input_) + force_term = (torch.sin(input_['x']*torch.pi) + * torch.sin(input_['y']*torch.pi)) + return gradgrad_u_x['x'] + gradgrad_u_y['y'] - force_term + + def nil_dirichlet(input_, output_): + value = 0.0 + return output_['u'] - value + + self.conditions = { + 'gamma1': {'location': Segment((0, 0), (1, 0)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment((1, 0), (1, 1)), 'func': nil_dirichlet}, + 'gamma3': {'location': Segment((1, 1), (0, 1)), 'func': nil_dirichlet}, + 'gamma4': {'location': Segment((0, 1), (0, 0)), 'func': nil_dirichlet}, + 'D': {'location': Cube([[0, 1], [0, 1]]), 'func': laplace_equation} + } + + def poisson_sol(x, y): + return -(np.sin(x*np.pi)*np.sin(y*np.pi))/(2*np.pi**2) + + self.input_variables = ['x', 'y'] + self.output_variables = ['u'] + self.truth_solution = poisson_sol + self.spatial_domain = Cube([[0, 1], [0, 1]]) + + #self.check() # Check the problem is correctly set diff --git a/problems/poisson_2.py b/problems/poisson_2.py new file mode 100644 index 0000000..cbafb7c --- /dev/null +++ b/problems/poisson_2.py @@ -0,0 +1,43 @@ +import numpy as np +import torch +from pina.segment import Segment +from pina.cube import Cube +from pina.problem2d import Problem2D +from pina.problem import Problem + + +class Poisson2DProblem(Problem2D): + + def __init__(self): + + def laplace_equation(input_, output_): + grad_u = self.grad(output_['u'], input_) + gradgrad_u_x = self.grad(grad_u['x'], input_) + gradgrad_u_y = self.grad(grad_u['y'], input_) + #force_term = (torch.sin(input_['x']*torch.pi) + # * torch.sin(input_['y']*torch.pi)) + force_term = -2*(input_['y']*(1-input_['y']) + + input_['x']*(1-input_['x'])) + return gradgrad_u_x['x'] + gradgrad_u_y['y'] - force_term + + def nil_dirichlet(input_, output_): + value = 0.0 + return output_['u'] - value + + self.conditions = { + 'gamma1': {'location': Segment((0, 0), (1, 0)), 'func': nil_dirichlet}, + 'gamma2': {'location': Segment((1, 0), (1, 1)), 'func': nil_dirichlet}, + 'gamma3': {'location': Segment((1, 1), (0, 1)), 'func': nil_dirichlet}, + 'gamma4': {'location': Segment((0, 1), (0, 0)), 'func': nil_dirichlet}, + 'D': {'location': Cube([[0, 1], [0, 1]]), 'func': laplace_equation} + } + + def poisson_sol(x, y): + return x*(1-x)*y*(1-y) + + self.input_variables = ['x', 'y'] + self.output_variables = ['u'] + self.truth_solution = poisson_sol + self.spatial_domain = Cube([[0, 1], [0, 1]]) + + #self.check() # Check the problem is correctly set