diff --git a/pinnicle/physics/constants.py b/pinnicle/physics/constants.py index 823ad17..e9bca03 100644 --- a/pinnicle/physics/constants.py +++ b/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 diff --git a/pinnicle/physics/stressbalance.py b/pinnicle/physics/stressbalance.py index 395f87a..4949267 100644 --- a/pinnicle/physics/stressbalance.py +++ b/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] #}}} - +# }}} diff --git a/pyproject.toml b/pyproject.toml index 1016e2e..5956a25 100644 --- a/pyproject.toml +++ b/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 diff --git a/tests/test_physics.py b/tests/test_physics.py index 57ecaa6..beb0de6 100644 --- a/tests/test_physics.py +++ b/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 = {} diff --git a/tests/test_pinn.py b/tests/test_pinn.py index 892a657..1a27cd2 100644 --- a/tests/test_pinn.py +++ b/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) diff --git a/tests/test_plotting.py b/tests/test_plotting.py index b8b0c9b..598a221 100644 --- a/tests/test_plotting.py +++ b/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")