source: issm/oecreview/Archive/24684-25833/ISSM-25805-25806.diff

../trunk-jpl/src/m/classes/frictionschoof.py

@@ +1 @@
+import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
+from project3d import project3d
+from structtoobj import structtoobj
+from WriteData import WriteData
+
+
+class frictionschoof(object):
+    """FRICTIONSCHOOF class definition
+
+    Usage:
+        friction = frictionschoof()
+    """
+
+    def __init__(self, *args):  # {{{
+        self.C                          = np.nan
+        self.Cmax                       = np.nan
+        self.m                          = np.nan
+        self.effective_pressure_limit   = 0
+        
+        nargs = len(args)
+        if nargs == 0:
+            self.setdefaultparameters()
+        elif nargs == 1:
+            self = structtoobj(self, args[0])
+        else:
+            raise Exception('constructor not supported')
+    #}}}
+
+    def __repr__(self):  # {{{
+        # See Brondex et al. 2019 https://www.the-cryosphere.net/13/177/2019/
+        s = 'Schoof sliding law parameters:\n'
+        s += '   Schoof\'s sliding law reads:\n'
+        s += '                         C^2 |u_b|^(m-1)                \n'
+        s += '      tau_b = - _____________________________   u_b   \n'
+        s += '               (1+(C^2/(Cmax N))^1/m |u_b| )^m          \n'
+        s += '\n'
+        s += "{}\n".format(fielddisplay(self, 'C', 'friction coefficient [SI]'))
+        s += "{}\n".format(fielddisplay(self, 'Cmax', 'Iken\'s bound (typically between 0.17 and 0.84) [SI]'))
+        s += "{}\n".format(fielddisplay(self, 'm', 'm exponent (generally taken as m = 1/n = 1/3)'))
+        s += "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'fNeff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'))
+        return s
+    #}}}
+
+    def setdefaultparameters(self):  # {{{
+        self.effective_pressure_limit = 0
+        return self
+    #}}}
+
+    def extrude(self, md):  # {{{
+        self.C = project3d(md, 'vector', self.C, 'type', 'node')
+        self.Cmax = project3d(md, 'vector', self.Cmax, 'type', 'node')
+        self.m = project3d(md, 'vector', self.m, 'type', 'element')
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses):  # {{{
+        # Early return
+        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
+            return md
+        md = checkfield(md, 'fieldname', 'friction.C', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>',0.)
+        md = checkfield(md, 'fieldname', 'friction.Cmax', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 0.)
+        md = checkfield(md, 'fieldname', 'friction.m', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofelements, 1])
+        md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
+        return md
+    # }}}
+
+    def marshall(self, prefix, md, fid):  # {{{
+        yts = md.constants.yts
+
+        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 11, 'format', 'Integer')
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'C', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'Cmax', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'm', 'format', 'DoubleMat', 'mattype', 2)
+        WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double')
+    # }}}
+

../trunk-jpl/src/m/classes/frictionshakti.py

@@ -27 +27 @@
     def __repr__(self):  # {{{
         s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n'
         s += '(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n'
-        s += "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
+        s += '{}\n'.format(fielddisplay(self, 'coefficient', 'friction coefficient [SI]'))
         return s
     #}}}
 

../trunk-jpl/src/m/classes/calving.py

@@ +1 @@
+import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
 from project3d import project3d
-from checkfield import checkfield
 from WriteData import WriteData
 
 
 class calving(object):
-    """
-    CALVING class definition
+    """CALVING class definition
 
-       Usage:
-          calving = calving()
+    Usage:
+        calving = calving()
     """
 
     def __init__(self):  # {{{
-
-        self.calvingrate = float('NaN')
-
-    #set defaults
+        self.calvingrate = np.nan
         self.setdefaultparameters()
 
     #}}}
 
     def __repr__(self):  # {{{
-        string = '   Calving parameters:'
-        string = "%s\n%s" % (string, fielddisplay(self, 'calvingrate', 'calving rate at given location [m / a]'))
-
-        return string
+        s = '   Calving parameters:'
+        s += '{}\n'.format(fielddisplay(self, 'calvingrate', 'calving rate at given location [m / a]'))
+        return s
     #}}}
 
     def extrude(self, md):  # {{{
@@ -39 +36 @@
 
     def checkconsistency(self, md, solution, analyses):  # {{{
         #Early return
-        if (solution != 'TransientSolution') or (not md.transient.ismovingfront):
+        if solution != 'TransientSolution' or not md.transient.ismovingfront:
             return md
 
         md = checkfield(md, 'fieldname', 'calving.calvingrate', '>=', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)

../trunk-jpl/src/m/classes/inversion.py

@@ -53 +53 @@
         s += '{}\n'.format(fielddisplay(self, 'step_threshold', 'decrease threshold for misfit, default is 30%'))
         s += '{}\n'.format(fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
         s += '{}\n'.format(fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
-        s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
-        s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
-        s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
+        s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m/yr]'))
+        s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m/yr]'))
+        s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m/yr]'))
         s += '{}\n'.format(fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
         s += '{}\n'.format(fielddisplay(self, 'surface_obs', 'observed surface elevation [m]'))
         s += '{}\n'.format('Available cost functions:')
@@ -159 +159 @@
         if not self.iscontrol:
             return
         WriteData(fid, prefix, 'object', self, 'fieldname', 'nsteps', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter_per_step', 'format', 'DoubleMat', 'mattype', 3)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter_per_step', 'format', 'IntMat', 'mattype', 3)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_functions_coefficients', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'gradient_scaling', 'format', 'DoubleMat', 'mattype', 3)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_function_threshold', 'format', 'Double')

../trunk-jpl/src/m/classes/solidearthsettings.py

@@ -21 +21 @@
         self.rotation               = 0
         self.ocean_area_scaling     = 0
         self.runfrequency           = 1 # How many time steps we skip before we run grd_core
-        self.computesealevelchange  = 0 # Will grd_core compute sea level?
+        self.computesealevelchange  = 1 # Will grd_core compute sea level?
         self.isgrd                  = 1 # Will GRD patterns be computed?
         self.degacc                 = 0 # Degree increment for resolution of Green tables
         self.horiz                  = 0 # Compute horizontal displacement?
@@ -65 +65 @@
         self.rotation = 1
         self.ocean_area_scaling = 0
         self.isgrd = 1
-        self.computesealevelchange = 0
+        self.computesealevelchange = 1
 
         # Numerical discretization accuracy
         self.degacc = .01
@@ -73 +73 @@
         # How many time steps we skip before we run solidearthsettings solver during transient
         self.runfrequency = 1
 
+        # Fractional contribution
+        self.glfraction = 1
+
         # Horizontal displacement? (not on by default)
         self.horiz = 0
-
-        # Fractional contribution
-        self.glfraction = 1
     #}}}
 
     def checkconsistency(self, md, solution, analyses): # {{{

../trunk-jpl/src/m/inversions/supportedcostfunctions.py

@@ -1 +1 @@
 def supportedcostfunctions():
-    return [101, 102, 103, 104, 105, 201, 501, 502, 503, 504, 505]
+    return list(range(101, 105 + 1)) + [201] + list(range(501, 508 + 1)) + [510] + list(range(601, 604 + 1))

../trunk-jpl/src/m/inversions/supportedcontrols.py

@@ -1 +1 @@
 def supportedcontrols():
-    return ['BalancethicknessThickeningRate', 'FrictionCoefficient', 'FrictionAs', 'MaterialsRheologyBbar', 'DamageDbar', 'Vx', 'Vy']
+    return [
+        'BalancethicknessThickeningRate',
+        'FrictionCoefficient',
+        'FrictionC',
+        'FrictionAs',
+        'MaterialsRheologyBbar',
+        'DamageDbar',
+        'Vx',
+        'Vy',
+        'Thickness',
+        'BalancethicknessOmega',
+        'BalancethicknessApparentMassbalance',
+        'MaterialsRheologyB'
+    ]

../trunk-jpl/src/m/qmu/postqmu.m

@@ -15 +15 @@
            disp(sprintf('%s',fline));
            fline=fgetl(fide);
        end
-       warning(['Dakota returned error in ''' qmuerrfile ' file. qmu directory retained.'])
+       warning(['Dakota returned error in ''' qmuerrfile ''' file. QMU directory retained.'])
    end
    status=fclose(fide);
 end

../trunk-jpl/test/NightlyRun/test356.py

@@ +1 @@
+#Test Name: TransientFrictionSchoof
+import numpy as np
+
+from frictionschoof import frictionschoof
+from MatlabFuncs import oshostname
+from model import *
+from parameterize import parameterize
+from setflowequation import setflowequation
+from setmask import setmask
+from solve import solve
+from transient import transient
+from triangle import triangle
+
+
+md = triangle(model(), '../Exp/Square.exp', 200000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareSheetConstrained.py')
+md = setflowequation(md, 'SSA', 'all')
+
+# Use Schoof's law
+Cmax = 0.8
+md.friction = frictionschoof()
+md.friction.m    = 1.0 / 3.0 * np.ones((md.mesh.numberofelements, 1))
+md.friction.Cmax = Cmax * np.ones((md.mesh.numberofvertices, 1))
+md.friction.C = 200 * np.ones((md.mesh.numberofvertices, 1))
+
+# Control parameters
+md.inversion.iscontrol = 1
+md.inversion.control_parameters = ['FrictionC']
+md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices, 1))
+md.inversion.max_parameters = 10000. * np.ones((md.mesh.numberofvertices, 1))
+md.inversion.nsteps = 2
+md.inversion.cost_functions = [102, 501]
+md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 2))
+md.inversion.cost_functions_coefficients[:, 1] = 2e-7
+md.inversion.gradient_scaling = 3. * np.ones((md.inversion.nsteps, 1))
+md.inversion.maxiter_per_step = 2 * np.ones((md.inversion.nsteps, 1))
+md.inversion.step_threshold = 0.3 * np.ones((md.inversion.nsteps, 1))
+md.inversion.vx_obs = md.initialization.vx
+md.inversion.vy_obs= md.initialization.vy
+
+md.cluster = generic('name', oshostname(), 'np', 3)
+md = solve(md, 'Stressbalance')
+
+#Fields and tolerances to track changes
+field_names = ['Gradient', 'Misfits', 'FrictionC', 'Pressure', 'Vel', 'Vx', 'Vy']
+field_tolerances = [1e-12, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
+field_values = [
+    md.results.StressbalanceSolution.Gradient1,
+    md.results.StressbalanceSolution.J,
+    md.results.StressbalanceSolution.FrictionC,
+    md.results.StressbalanceSolution.Pressure,
+    md.results.StressbalanceSolution.Vel,
+    md.results.StressbalanceSolution.Vx,
+    md.results.StressbalanceSolution.Vy
+]

../trunk-jpl/test/NightlyRun/test481.py

@@ +1 @@
+#Test Name: TransientFrictionSchoof
+import numpy as np
+
+from frictionschoof import frictionschoof
+from MatlabFuncs import oshostname
+from model import *
+from parameterize import parameterize
+from setflowequation import setflowequation
+from setmask import setmask
+from solve import solve
+from transient import transient
+from triangle import triangle
+
+md = triangle(model(), '../Exp/Square.exp', 150000.)
+md = setmask(md, '../Exp/SquareShelf.exp', '')
+md = parameterize(md, '../Par/SquareSheetShelf.py')
+md.extrude(4, 1)
+md = setflowequation(md, 'HO', 'all')
+md.transient.isthermal = 0
+md.friction = frictionschoof(md.friction)
+md.friction.C    = pow(20.e4, 0.5) * np.ones((md.mesh.numberofvertices, 1))
+md.friction.Cmax = 0.5 * np.ones((md.mesh.numberofvertices, 1))
+md.friction.m    = 1./3.* np.ones((md.mesh.numberofelements, 1))
+md.cluster = generic('name', oshostname(), 'np', 3)
+md = solve(md, 'Transient')
+
+#Fields and tolerances to track changes
+field_names = [
+    'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1',
+    'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2',
+    'Vx3', 'Vy3', 'Vel3', 'Pressure3', 'Bed3', 'Surface3', 'Thickness3'
+]
+field_tolerances = [
+    2e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09,
+    1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09,
+    2e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09
+]
+field_values = [
+    md.results.TransientSolution[0].Vx,
+    md.results.TransientSolution[0].Vy,
+    md.results.TransientSolution[0].Vel,
+    md.results.TransientSolution[0].Pressure,
+    md.results.TransientSolution[0].Base,
+    md.results.TransientSolution[0].Surface,
+    md.results.TransientSolution[0].Thickness,
+    md.results.TransientSolution[1].Vx,
+    md.results.TransientSolution[1].Vy,
+    md.results.TransientSolution[1].Vel,
+    md.results.TransientSolution[1].Pressure,
+    md.results.TransientSolution[1].Base,
+    md.results.TransientSolution[1].Surface,
+    md.results.TransientSolution[1].Thickness,
+    md.results.TransientSolution[2].Vx,
+    md.results.TransientSolution[2].Vy,
+    md.results.TransientSolution[2].Vel,
+    md.results.TransientSolution[2].Pressure,
+    md.results.TransientSolution[2].Base,
+    md.results.TransientSolution[2].Surface,
+    md.results.TransientSolution[2].Thickness
+]

../trunk-jpl/test/NightlyRun/test350.py

@@ -1 +1 @@
 #Test Name: SquareSheetHydrologyShakti
-from operator import itemgetter
 import numpy as np
+
+from frictionshakti import frictionshakti
+from MatlabFuncs import oshostname
 from model import *
-from socket import gethostname
-from triangle import triangle
-from setmask import setmask
+from operator import itemgetter
 from parameterize import parameterize
 from setflowequation import setflowequation
+from setmask import setmask
 from solve import solve
-from frictionshakti import frictionshakti
 from transient import transient
+from triangle import triangle
 
 md = triangle(model(), '../Exp/Square.exp', 50000.)
 md.mesh.x = md.mesh.x / 1000

../trunk-jpl/test/NightlyRun/runme.m

@@ -106 +106 @@
 % }}}
 %Process Ids according to benchmarks{{{
 if strcmpi(benchmark,'nightly'),
-        test_ids=intersect(test_ids,[1:999]);
+        test_ids=intersect(test_ids,union([1:999]));
 elseif strcmpi(benchmark,'validation'),
         test_ids=intersect(test_ids,[1001:1999]);
 elseif strcmpi(benchmark,'ismip'),

../trunk-jpl/test/NightlyRun/test444.py

@@ -6 +6 @@
 #
 
 import numpy as np
+
+from ContourToMesh import *
+from dmeth_params_set import *
+from MatlabFuncs import oshostname
 from model import *
-from triangle import *
-from setmask import *
 from parameterize import *
+from partitioner import *
+from regionaloutput import *
 from setflowequation import *
+from setmask import *
 from solve import *
-from partitioner import *
-from dmeth_params_set import *
-from ContourToMesh import *
-from regionaloutput import *
+from triangle import *
 
 #model not consistent:  equality thickness = surface-base violated