Merge pull request #39 from ISSMteam/add_variable_rheology_B

add variable rheology B, functions `pde` in physics can not be covered by codecov due to tensorflow function

pinnicle/physics/constants.py

@@ -9,8 +9,8 @@ def __init__(self, **kwargs):
         self.yts    = 3600.0*24*365     # year to second (s)
 
         # Typical range of variables
-        self.variable_lb = {'u': -1.0e4/self.yts, 'v':-1.0e4/self.yts, 's':-1.0e3, 'H':10.0,   'C':0.01,  'a': -5.0/self.yts}
-        self.variable_ub = {'u': 1.0e4/self.yts,  'v':1.0e4/self.yts,  's':3.6e3,  'H':3500.0, 'C':1.0e4, 'a':  5.0/self.yts}
+        self.variable_lb = {'u': -1.0e4/self.yts, 'v':-1.0e4/self.yts, 's':-1.0e3, 'H':10.0,   'C':0.01,  'a': -5.0/self.yts, 'B': 7.0e7}
+        self.variable_ub = {'u': 1.0e4/self.yts,  'v':1.0e4/self.yts,  's':3.6e3,  'H':3500.0, 'C':1.0e4, 'a':  5.0/self.yts, 'B': 7.0e8}
         # add more if needed
         self.variable_lb['u_base'] = -1.0e4/self.yts
         self.variable_ub['u_base'] =  1.0e4/self.yts

pinnicle/physics/stressbalance.py

@@ -2,6 +2,7 @@
 from . import EquationBase, Constants
 from ..parameter import EquationParameter
 
+# SSA constant B {{{
 class SSAEquationParameter(EquationParameter, Constants):
     """ default parameters for SSA
     """
@@ -25,7 +26,6 @@ def set_default(self):
                 'n': 3.0,               # exponent of Glen's flow law
                 'B':1.26802073401e+08   # -8 degree C, cuffey
                 }
-
 class SSA(EquationBase): #{{{
     """ SSA on 2D problem with uniform B
     """
@@ -83,7 +83,92 @@ def pde(self, nn_input_var, nn_output_var):
         f2 = sigma21 + sigma22 - alpha*v/(u_norm+1e-30) - self.rhoi*self.g*H*s_y
 
         return [f1, f2] #}}}
+#}}}
+# SSA variable B {{{
+class SSAVariableBEquationParameter(EquationParameter, Constants):
+    """ default parameters for SSA, with spatially varying rheology B
+    """
+    _EQUATION_TYPE = 'SSA_VB' 
+    def __init__(self, param_dict={}):
+        # load necessary constants
+        Constants.__init__(self)
+        super().__init__(param_dict)
+
+    def set_default(self):
+        self.input = ['x', 'y']
+        self.output = ['u', 'v', 's', 'H', 'C', 'B']
+        self.output_lb = [self.variable_lb[k] for k in self.output]
+        self.output_ub = [self.variable_ub[k] for k in self.output]
+        self.data_weights = [1.0e-8*self.yts**2.0, 1.0e-8*self.yts**2.0, 1.0e-6, 1.0e-6, 1.0e-8, 1e-16]
+        self.residuals = ["f"+self._EQUATION_TYPE+"1", "f"+self._EQUATION_TYPE+"2"]
+        self.pde_weights = [1.0e-10, 1.0e-10]
+
+        # scalar variables: name:value
+        self.scalar_variables = {
+                'n': 3.0,               # exponent of Glen's flow law
+                }
+class SSAVariableB(EquationBase): 
+    """ SSA on 2D problem with spatially varying B
+    """
+    _EQUATION_TYPE = 'SSA_VB' 
+    def __init__(self, parameters=SSAVariableBEquationParameter()):
+        super().__init__(parameters)
+
+    def pde(self, nn_input_var, nn_output_var):
+        """ residual of SSA 2D PDEs
+
+        Args:
+            nn_input_var: global input to the nn
+            nn_output_var: global output from the nn
+        """
+        # no cover: start
+        # get the ids
+        xid = self.local_input_var["x"]
+        yid = self.local_input_var["y"]
+
+        uid = self.local_output_var["u"]
+        vid = self.local_output_var["v"]
+        sid = self.local_output_var["s"]
+        Hid = self.local_output_var["H"]
+        Cid = self.local_output_var["C"]
+        Bid = self.local_output_var["B"]
 
+        # unpacking normalized output
+        u, v, H, C, B = nn_output_var[:, uid:uid+1], nn_output_var[:, vid:vid+1], nn_output_var[:, Hid:Hid+1], nn_output_var[:, Cid:Cid+1], nn_output_var[:, Bid:Bid+1]
+
+        # spatial derivatives
+        u_x = dde.grad.jacobian(nn_output_var, nn_input_var, i=uid, j=xid)
+        v_x = dde.grad.jacobian(nn_output_var, nn_input_var, i=vid, j=xid)
+        s_x = dde.grad.jacobian(nn_output_var, nn_input_var, i=sid, j=xid)
+        u_y = dde.grad.jacobian(nn_output_var, nn_input_var, i=uid, j=yid)
+        v_y = dde.grad.jacobian(nn_output_var, nn_input_var, i=vid, j=yid)
+        s_y = dde.grad.jacobian(nn_output_var, nn_input_var, i=sid, j=yid)
+
+        eta = 0.5*B *(u_x**2.0 + v_y**2.0 + 0.25*(u_y+v_x)**2.0 + u_x*v_y+1.0e-15)**(0.5*(1.0-self.n)/self.n)
+        # stress tensor
+        etaH = eta * H
+        B11 = etaH*(4*u_x + 2*v_y)
+        B22 = etaH*(4*v_y + 2*u_x)
+        B12 = etaH*(  u_y +   v_x)
+
+        # Getting the other derivatives
+        sigma11 = dde.grad.jacobian(B11, nn_input_var, i=0, j=xid)
+        sigma12 = dde.grad.jacobian(B12, nn_input_var, i=0, j=yid)
+
+        sigma21 = dde.grad.jacobian(B12, nn_input_var, i=0, j=xid)
+        sigma22 = dde.grad.jacobian(B22, nn_input_var, i=0, j=yid)
+
+        # compute the basal stress
+        u_norm = (u**2+v**2)**0.5
+        alpha = C**2 * (u_norm)**(1.0/self.n)
+
+        f1 = sigma11 + sigma12 - alpha*u/(u_norm+1e-30) - self.rhoi*self.g*H*s_x
+        f2 = sigma21 + sigma22 - alpha*v/(u_norm+1e-30) - self.rhoi*self.g*H*s_y
+
+        return [f1, f2] # }}} 
+        # no cover: stop
+#}}}
+# MOLHO constant B{{{
 class MOLHOEquationParameter(EquationParameter, Constants):
     """ default parameters for MOLHO
     """
@@ -208,4 +293,4 @@ def pde(self, nn_input_var, nn_output_var):
         f4 = sigma41 + sigma42 + mu4*vshear - self.rhoi*self.g*H*s_y*(self.n+1.0)/(self.n+2.0)
 
         return [f1, f2, f3, f4] #}}}
-
+# }}}

pyproject.toml

@@ -52,19 +52,21 @@ exclude = ["DATA", "Models", "docs*", "examples*"]
 [tool.coverage.report]
 # Regexes for lines to exclude from consideration
 exclude_also = [
-    # Don't complain about missing debug-only code:
-    "def __repr__",
-    "if self\\.debug",
+  # 1. Comments to turn coverage on and off:
+  "no cover: start(?s:.)*?no cover: stop",
+  # Don't complain about missing debug-only code:
+  "def __repr__",
+  "if self\\.debug",
 
-    # Don't complain if tests don't hit defensive assertion code:
-    "raise AssertionError",
-    "raise NotImplementedError",
+  # Don't complain if tests don't hit defensive assertion code:
+  "raise AssertionError",
+  "raise NotImplementedError",
 
-    # Don't complain about abstract methods, they aren't run:
-    "@(abc\\.)?abstractmethod",
+  # Don't complain about abstract methods, they aren't run:
+  "@(abc\\.)?abstractmethod",
 
-	 # Don't complain pytest.mark.skip:
-	 "@pytest.mark.skip",
-    ]
+  # Don't complain pytest.mark.skip:
+  "@pytest.mark.skip",
+]
 
 ignore_errors = true

tests/test_physics.py

@@ -52,6 +52,21 @@ def test_Physics_SSA():
     assert len(phy.data_weights) == 5
     assert len(phy.pde_weights) == 2
 
+def test_Physics_SSAVB():
+    SSA = {}
+    SSA["scalar_variables"] = {"n":3}
+    hp = {}
+    hp["equations"] = {"SSA_VB":SSA}
+    phy = Physics(PhysicsParameter(hp))
+
+    assert phy.input_var == ['x', 'y']
+    assert phy.output_var == ['u', 'v', 's', 'H', 'C','B']
+    assert phy.residuals == ['fSSA_VB1', 'fSSA_VB2']
+    assert len(phy.output_lb) == 6
+    assert len(phy.output_ub) == 6
+    assert len(phy.data_weights) == 6
+    assert len(phy.pde_weights) == 2
+
 def test_Physics_MC():
     MC = {}
     MC["scalar_variables"] = {"B":1.26802073401e+08}
@@ -152,13 +167,17 @@ def test_operator():
     SSA["scalar_variables"] = {"B":1.26802073401e+08}
 
     hp = {}
-    hp["equations"] = {"MC":MC, 'SSA':SSA}
+    hp["equations"] = {"MC":MC, "SSA":SSA, "SSA_VB":{}, "MOLHO":{}}
     phy = Physics(PhysicsParameter(hp))
 
     assert phy.operator('mc')
     assert phy.operator('Mc')
     assert phy.operator('SSA')
     assert phy.operator('ssa')
+    assert phy.operator('SSA_VB')
+    assert phy.operator('ssa_vb')
+    assert phy.operator('molho')
+    assert phy.operator('MOLHO')
 
 def test_Physics_dummy():
     dummy = {}

tests/test_pinn.py

@@ -2,7 +2,7 @@
 import pinnicle as pinn
 import numpy as np
 import deepxde as dde
-from pinnicle.utils import data_misfit, plot_nn, diffplot
+from pinnicle.utils import data_misfit, plot_nn
 import pytest
 
 dde.config.set_default_float('float64')
@@ -199,7 +199,7 @@ def test_only_callbacks(tmp_path):
 def test_plot(tmp_path):
     hp["save_path"] = str(tmp_path)
     hp["is_save"] = True
-    issm["data_size"] = {"u":4000, "v":4000, "s":4000, "H":4000, "C":None}
+    issm["data_size"] = {"u":10, "v":10, "s":10, "H":10, "C":None}
     hp["data"] = {"ISSM": issm}
     experiment = pinn.PINN(params=hp)
     experiment.compile()
@@ -217,18 +217,57 @@ def test_plot(tmp_path):
     assert len(im_data) == 5
     assert im_data['u'].shape == (10,10) 
 
-def test_diffplot(tmp_path):
-    hp["save_path"] = str(tmp_path)
-    hp["is_save"] = True
-    issm["data_size"] = {"u":4000, "v":4000, "s":4000, "H":4000, "C":None}
+def test_SSA_pde_function():
+    SSA = {}
+    SSA["n"] = {"n":3}
+    hp["equations"] = {"SSA":SSA}
+    hp["num_collocation_points"] = 10
+    issm["data_size"] = {"u":10, "v":10, "s":10, "H":10, "C":None, "vel":10}
     hp["data"] = {"ISSM": issm}
     experiment = pinn.PINN(params=hp)
     experiment.compile()
+    y = experiment.model.predict(experiment.model_data.X['u'], operator=experiment.physics.operator("SSA"))
+    assert len(y) == 2
+    assert y[0].shape == (10,1)
+    assert y[1].shape == (10,1)
 
-    fig, axs = diffplot(experiment, 'H')
-    assert fig is not None
-    assert axs.shape == (3,)
-    fig, axs = diffplot(experiment, ['u', 'v'], feat_title='vel')
-    assert fig is not None
-    assert axs.shape == (3,)
+def test_SSA_VB_pde_function():
+    SSA = {}
+    SSA["n"] = {"n":3}
+    hp["equations"] = {"SSA_VB":SSA}
+    hp["num_collocation_points"] = 10
+    issm["data_size"] = {"u":10, "v":10, "s":10, "H":10, "C":None, "vel":10}
+    hp["data"] = {"ISSM": issm}
+    experiment = pinn.PINN(params=hp)
+    experiment.compile()
+    y = experiment.model.predict(experiment.model_data.X['u'], operator=experiment.physics.operator("SSA_VB"))
+    assert len(y) == 2
+    assert y[0].shape == (10,1)
+    assert y[1].shape == (10,1)
 
+def test_MOLHO_pde_function():
+    MOLHO = {}
+    MOLHO["n"] = {"n":3}
+    hp["equations"] = {"MOLHO":MOLHO}
+    hp["num_collocation_points"] = 10
+    issm["data_size"] = {"u":10, "v":10, "s":10, "H":10, "C":None, "vel":10}
+    hp["data"] = {"ISSM": issm}
+    experiment = pinn.PINN(params=hp)
+    experiment.compile()
+    y = experiment.model.predict(experiment.model_data.X['u'], operator=experiment.physics.operator("MOLHO"))
+    assert len(y) == 4
+    assert y[0].shape == (10,1)
+    assert y[3].shape == (10,1)
+
+def test_MC_pde_function():
+    MC = {}
+    MC["n"] = {"n":3}
+    hp["equations"] = {"MC":MC}
+    hp["num_collocation_points"] = 10
+    issm["data_size"] = {"u":10, "v":10, "s":10, "H":10, "C":None, "vel":10}
+    hp["data"] = {"ISSM": issm}
+    experiment = pinn.PINN(params=hp)
+    experiment.compile()
+    y = experiment.model.predict(experiment.model_data.X['u'], operator=experiment.physics.operator("MC"))
+    assert len(y) == 1
+    assert y[0].shape == (10,1)

tests/test_plotting.py

@@ -2,7 +2,7 @@
 import pinnicle as pinn
 import numpy as np
 import deepxde as dde
-from pinnicle.utils import plot_similarity, plot_residuals, tripcolor_similarity, tripcolor_residuals
+from pinnicle.utils import plot_similarity, plot_residuals, tripcolor_similarity, tripcolor_residuals, diffplot
 import matplotlib.pyplot as plt
 import pytest
 
@@ -155,6 +155,20 @@ def test_triresiduals(tmp_path):
     assert (fig is not None) and (np.size(axs)==2)
     fig, axs = tripcolor_residuals(experiment, cbar_limits=[-7e3, 7e3])
     assert (fig is not None) and (np.size(axs)==2)
-
     plt.close("all") 
 
+def test_diffplot(tmp_path):
+    hp["save_path"] = str(tmp_path)
+    hp["is_save"] = True
+    issm["data_size"] = {"u":100, "v":100, "s":100, "H":100, "C":None}
+    hp["data"] = {"ISSM": issm}
+    experiment = pinn.PINN(params=hp)
+    experiment.compile()
+
+    fig, axs = diffplot(experiment, 'H')
+    assert fig is not None
+    assert axs.shape == (3,)
+    fig, axs = diffplot(experiment, ['u', 'v'], feat_title='vel')
+    assert fig is not None
+    assert axs.shape == (3,)
+    plt.close("all")