Refactoring code

examples/problems/burgers.py

@@ -0,0 +1,49 @@
+import numpy as np
+import scipy.io
+import torch
+
+from pina.segment import Segment
+from pina.cube import Cube
+from pina.problem import TimeDependentProblem, Problem1D
+from pina.operators import grad
+
+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(TimeDependentProblem, Problem1D):
+
+    input_variables = ['x', 't']
+    output_variables = ['u']
+    spatial_domain = Cube([[-1, 1]])
+    temporal_domain = Cube([[0, 1]])
+
+    def burger_equation(input_, output_):
+        grad_u = grad(output_['u'], input_)
+        grad_x, grad_t = tmp_grad(output_['u'], input_).T
+        gradgrad_u_x = grad(grad_u['x'], input_)
+        grad_xx = tmp_grad(grad_x, input_)[:, 0]
+        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
+
+
+
+    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}
+    }

examples/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]])
+

examples/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]])
+

examples/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]])
+

examples/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

examples/problems/poisson2D.py

@@ -0,0 +1,35 @@
+import numpy as np
+import torch
+from pina.segment import Segment
+from pina.cube import Cube
+from pina.problem import Problem2D
+from pina.operators import grad, div, nabla
+
+
+class Poisson2D(Problem2D):
+
+    input_variables = ['x', 'y']
+    output_variables = ['u']
+    spatial_domain = Cube([[0, 1], [0, 1]])
+
+    def laplace_equation(input_, output_):
+        force_term = (torch.sin(input_['x']*torch.pi)
+                        * torch.sin(input_['y']*torch.pi))
+        return nabla(output_['u'], input_).flatten() - force_term
+
+    def nil_dirichlet(input_, output_):
+        value = 0.0
+        return output_['u'] - value
+
+    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(self, x, y):
+        return -(np.sin(x*np.pi)*np.sin(y*np.pi))/(2*np.pi**2)
+
+    truth_solution = poisson_sol

examples/run_burgers.py

@@ -9,34 +9,20 @@ 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
+from pina import Plotter
 
 
 class myFeature(torch.nn.Module):
     """
-    Feature: sin(x)
+    Feature: sin(pi*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")
@@ -51,13 +37,13 @@ if __name__ == "__main__":
 
     burgers_problem = Burgers1D()
     model = DeepFeedForward(
-        layers=[20, 10, 5],
-        #layers=[8, 4, 2],
+        layers=[30, 20, 10, 5],
+        #layers=[8, 8, 8],
         #layers=[16, 8, 4, 4],
         #layers=[20, 4, 4, 4],
         output_variables=burgers_problem.output_variables,
         input_variables=burgers_problem.input_variables,
-        func=Tanh,
+        func=Softplus,
         extra_features=feat
     )
 
@@ -70,46 +56,11 @@ if __name__ == "__main__":
         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.span_pts(2000, 'latin', ['D'])
+        pinn.span_pts(150, 'random', ['gamma1', 'gamma2', 'initia'])
+        pinn.train(5000, 100)
         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()
+        plotter = Plotter()
+        plotter.plot(pinn)

examples/run_poisson.py

@@ -2,15 +2,16 @@ import sys
 import numpy as np
 import torch
 import argparse
-from pina.pinn import PINN
+from pina 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 problems.poisson2D import Poisson2D
 from pina.deep_feed_forward import DeepFeedForward
 
 from pina.adaptive_functions import AdaptiveSin, AdaptiveCos, AdaptiveTanh
 
+from pina import Plotter
 
 class myFeature(torch.nn.Module):
     """
@@ -54,17 +55,18 @@ if __name__ == "__main__":
 
     if args.s:
 
-        pinn.span_pts(10, 'grid', ['D'])
-        pinn.span_pts(10, 'grid', ['gamma1', 'gamma2', 'gamma3', 'gamma4'])
+        pinn.span_pts(20, 'grid', ['D'])
+        pinn.span_pts(20, 'grid', ['gamma1', 'gamma2', 'gamma3', 'gamma4'])
         #pinn.plot_pts()
-        pinn.train(10000, 100)
+        pinn.train(1000, 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()))
+                file_.write('{} {}\n'.format(i, sum(losses)))
         pinn.save_state('pina.poisson')
 
     else:
         pinn.load_state('pina.poisson')
-        pinn.plot(40)
+        plotter = Plotter()
+        plotter.plot(pinn)