Index: /issm/trunk-jpl/src/m/classes/calving.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calving.py (revision 25805) +++ /issm/trunk-jpl/src/m/classes/calving.py (revision 25806) @@ -1,21 +1,19 @@ +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() @@ -23,8 +21,7 @@ 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 #}}} @@ -40,5 +37,5 @@ 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 Index: /issm/trunk-jpl/src/m/classes/frictionschoof.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionschoof.py (revision 25806) +++ /issm/trunk-jpl/src/m/classes/frictionschoof.py (revision 25806) @@ -0,0 +1,79 @@ +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') + # }}} + Index: /issm/trunk-jpl/src/m/classes/frictionshakti.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionshakti.py (revision 25805) +++ /issm/trunk-jpl/src/m/classes/frictionshakti.py (revision 25806) @@ -28,5 +28,5 @@ 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 #}}} Index: /issm/trunk-jpl/src/m/classes/inversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.py (revision 25805) +++ /issm/trunk-jpl/src/m/classes/inversion.py (revision 25806) @@ -54,7 +54,7 @@ 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]')) @@ -160,5 +160,5 @@ 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) Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25805) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25806) @@ -22,5 +22,5 @@ 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 @@ -66,5 +66,5 @@ self.ocean_area_scaling = 0 self.isgrd = 1 - self.computesealevelchange = 0 + self.computesealevelchange = 1 # Numerical discretization accuracy @@ -74,9 +74,9 @@ self.runfrequency = 1 + # Fractional contribution + self.glfraction = 1 + # Horizontal displacement? (not on by default) self.horiz = 0 - - # Fractional contribution - self.glfraction = 1 #}}} Index: /issm/trunk-jpl/src/m/inversions/supportedcontrols.py =================================================================== --- /issm/trunk-jpl/src/m/inversions/supportedcontrols.py (revision 25805) +++ /issm/trunk-jpl/src/m/inversions/supportedcontrols.py (revision 25806) @@ -1,2 +1,15 @@ def supportedcontrols(): - return ['BalancethicknessThickeningRate', 'FrictionCoefficient', 'FrictionAs', 'MaterialsRheologyBbar', 'DamageDbar', 'Vx', 'Vy'] + return [ + 'BalancethicknessThickeningRate', + 'FrictionCoefficient', + 'FrictionC', + 'FrictionAs', + 'MaterialsRheologyBbar', + 'DamageDbar', + 'Vx', + 'Vy', + 'Thickness', + 'BalancethicknessOmega', + 'BalancethicknessApparentMassbalance', + 'MaterialsRheologyB' + ] Index: /issm/trunk-jpl/src/m/inversions/supportedcostfunctions.py =================================================================== --- /issm/trunk-jpl/src/m/inversions/supportedcostfunctions.py (revision 25805) +++ /issm/trunk-jpl/src/m/inversions/supportedcostfunctions.py (revision 25806) @@ -1,2 +1,2 @@ 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)) Index: /issm/trunk-jpl/src/m/qmu/postqmu.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/postqmu.m (revision 25805) +++ /issm/trunk-jpl/src/m/qmu/postqmu.m (revision 25806) @@ -16,5 +16,5 @@ 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); Index: /issm/trunk-jpl/test/NightlyRun/runme.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.m (revision 25805) +++ /issm/trunk-jpl/test/NightlyRun/runme.m (revision 25806) @@ -107,5 +107,5 @@ %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]); Index: /issm/trunk-jpl/test/NightlyRun/test350.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test350.py (revision 25805) +++ /issm/trunk-jpl/test/NightlyRun/test350.py (revision 25806) @@ -1,14 +1,15 @@ #Test Name: SquareSheetHydrologyShakti +import numpy as np + +from frictionshakti import frictionshakti +from MatlabFuncs import oshostname +from model import * from operator import itemgetter -import numpy as np -from model import * -from socket import gethostname -from triangle import triangle -from setmask import setmask 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.) Index: /issm/trunk-jpl/test/NightlyRun/test356.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test356.py (revision 25806) +++ /issm/trunk-jpl/test/NightlyRun/test356.py (revision 25806) @@ -0,0 +1,56 @@ +#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 +] Index: /issm/trunk-jpl/test/NightlyRun/test444.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test444.py (revision 25805) +++ /issm/trunk-jpl/test/NightlyRun/test444.py (revision 25806) @@ -7,14 +7,16 @@ import numpy as np + +from ContourToMesh import * +from dmeth_params_set import * +from MatlabFuncs import oshostname from model import * +from parameterize import * +from partitioner import * +from regionaloutput import * +from setflowequation import * +from setmask import * +from solve import * from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * -from partitioner import * -from dmeth_params_set import * -from ContourToMesh import * -from regionaloutput import * #model not consistent: equality thickness = surface-base violated Index: /issm/trunk-jpl/test/NightlyRun/test481.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test481.py (revision 25806) +++ /issm/trunk-jpl/test/NightlyRun/test481.py (revision 25806) @@ -0,0 +1,60 @@ +#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 +]