Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26186) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26187) @@ -921,7 +921,13 @@ break; case 11: + iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum); iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum); iomodel->FetchDataToInput(inputs,elements,"md.friction.Cmax",FrictionCmaxEnum); + if(FrictionCoupling==3){ + iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);} + else if(FrictionCoupling==4){ + iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum); + } break; case 12: @@ -938,5 +944,4 @@ else if(FrictionCoupling==4){ iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum); - } break; @@ -1047,5 +1052,5 @@ break; case 11: - parameters->AddObject(new IntParam(FrictionCouplingEnum,2)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum)); break; Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26186) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 26187) @@ -616,5 +616,5 @@ if(fread(record_name,record_name_size*sizeof(char),1,fid)!=0){}; - _error_("error while looking in binary file. String " << record_name << " a string of size "<=',0); - switch self.coupling - case 0 - case 1 - case 2 - case 3 - md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); - case 4 - - otherwise - error('not supported yet'); + if self.coupling==3 + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); end end % }}} @@ -92,14 +76,6 @@ WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer'); WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double'); - switch self.coupling - case 0 - case 1 - case 2 - case 3 - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - case 4 - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - otherwise - error('not supported yet'); + if self.coupling==3 || self.coupling==4 + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end end % }}} Index: /issm/trunk-jpl/src/m/classes/friction.py =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.py (revision 26186) +++ /issm/trunk-jpl/src/m/classes/friction.py (revision 26187) @@ -23,5 +23,4 @@ self.setdefaultparameters() #}}} - def __repr__(self): # {{{ s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n' @@ -35,27 +34,20 @@ return s #}}} - def setdefaultparameters(self): # {{{ self.effective_pressure_limit = 0 return self #}}} - def extrude(self, md): # {{{ self.coefficient = project3d(md, 'vector', self.coefficient, 'type', 'node', 'layer', 1) self.p = project3d(md, 'vector', self.p, 'type', 'element') self.q = project3d(md, 'vector', self.q, 'type', 'element') - # If self.coupling == 0: #doesnt work with empty loop, so just skip it? if self.coupling in[3, 4]: self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1) - elif self.coupling > 4: - raise ValueError('md.friction.coupling larger than 4, not supported yet') return self #}}} - def defaultoutputs(self, md): # {{{ list = [] return list #}}} - def checkconsistency(self, md, solution, analyses): # {{{ # Early return @@ -71,9 +63,6 @@ if self.coupling == 3: md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - elif self.coupling > 4: - raise ValueError('not supported yet') return md # }}} - def marshall(self, prefix, md, fid): # {{{ WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer') @@ -89,9 +78,5 @@ WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coupling', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double') - if self.coupling == 3: + if self.coupling == 3 or self.coupling == 4: WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) - elif self.coupling == 4: - WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) - elif self.coupling > 4: - raise ValueError('not supported yet') # }}} Index: /issm/trunk-jpl/src/m/classes/frictionschoof.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 26186) +++ /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 26187) @@ -6,7 +6,9 @@ classdef frictionschoof properties (SetAccess=public) - C = NaN; - Cmax = NaN; - m = NaN; + C = NaN; + Cmax = NaN; + m = NaN; + coupling = 0; + effective_pressure = NaN; effective_pressure_limit = 0; end @@ -26,7 +28,11 @@ self.Cmax = project3d(md,'vector',self.Cmax,'type','node'); self.m = project3d(md,'vector',self.m,'type','element'); + if self.coupling==3 || self.coupling==4 + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1); + end end % }}} function self = setdefaultparameters(self) % {{{ + self.coupling = 0; self.effective_pressure_limit = 0; @@ -36,8 +42,12 @@ %Early return if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end - md = checkfield(md,'fieldname','friction.C','timeseries',1,'NaN',1,'Inf',1,'>',0.); + 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); + md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0:4]); + if self.coupling==3 + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1); + end end % }}} function disp(self) % {{{ @@ -52,4 +62,6 @@ fielddisplay(self,'Cmax','Iken''s bound (typically between 0.17 and 0.84) [SI]'); fielddisplay(self,'m','m exponent (generally taken as m = 1/n = 1/3)'); + fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]'); + fielddisplay(self,'coupling','Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: use coupled model (not implemented yet)'); fielddisplay(self,'effective_pressure_limit','Neff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'); end % }}} @@ -62,4 +74,8 @@ 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'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer'); + if self.coupling==3 || self.coupling==4 + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + end end % }}} end Index: /issm/trunk-jpl/src/m/classes/frictionschoof.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionschoof.py (revision 26186) +++ /issm/trunk-jpl/src/m/classes/frictionschoof.py (revision 26187) @@ -16,8 +16,10 @@ def __init__(self, *args): # {{{ - self.C = np.nan - self.Cmax = np.nan - self.m = np.nan - self.effective_pressure_limit = 0 + self.C = np.nan + self.Cmax = np.nan + self.m = np.nan + self.coupling = 0 + self.effective_pressure = np.nan + self.effective_pressure_limit = 0 nargs = len(args) @@ -29,5 +31,4 @@ raise Exception('constructor not supported') #}}} - def __repr__(self): # {{{ # See Brondex et al. 2019 https://www.the-cryosphere.net/13/177/2019/ @@ -41,20 +42,21 @@ 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, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)')) + s += '{}\n'.format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]')) 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') + if self.coupling in[3, 4]: + self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1) return self #}}} - def checkconsistency(self, md, solution, analyses): # {{{ # Early return @@ -65,7 +67,9 @@ 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) + md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2, 3, 4]) + if self.coupling == 3: + md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1) return md # }}} - def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts @@ -75,5 +79,7 @@ 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, 'class', 'friction', 'object', self, 'fieldname', 'coupling', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double') + if self.coupling == 3 or self.coupling == 4: + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) # }}} -