Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictioncoulomb.m (revision 28029) +++ /issm/trunk-jpl/src/m/classes/frictioncoulomb.m (revision 28030) @@ -6,10 +6,10 @@ classdef frictioncoulomb properties (SetAccess=public) - coefficient = NaN; - p = NaN; - q = NaN; - coefficientcoulomb = NaN; - coupling = 0; - effective_pressure = NaN; + coefficient = NaN; + p = NaN; + q = NaN; + coefficientcoulomb = NaN; + coupling = 0; + effective_pressure = NaN; effective_pressure_limit = 0; end @@ -27,5 +27,5 @@ error('not implemented yet'); otherwise - error('not supported yet'); + error('not supported yet'); end end % }}} @@ -52,5 +52,6 @@ md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]); - md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0); + md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0); + switch self.coupling case 0 @@ -60,5 +61,5 @@ error('not implemented yet'); otherwise - error('not supported yet'); + error('not supported yet'); end end % }}} @@ -71,10 +72,11 @@ fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]'); fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(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)'); + 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 % }}} function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; WriteData(fid,prefix,'name','md.friction.law','data',7,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); WriteData(fid,prefix,'object',self,'fieldname','coefficientcoulomb','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2); @@ -85,9 +87,9 @@ case 0 case 1 - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); case 2 error('not implemented yet'); otherwise - error('not supported yet'); + error('not supported yet'); end Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictioncoulomb.py (revision 28029) +++ /issm/trunk-jpl/src/m/classes/frictioncoulomb.py (revision 28030) @@ -7,6 +7,6 @@ -class frictioncoulomb(object): - """FRICTIONCOULOMB class definition +class frictioncoulomb(): + """frictioncoulomb class definition Usage: @@ -14,5 +14,5 @@ """ - def __init__(self): # {{{ + def __init__(self, *args): # {{{ self.coefficient = np.nan self.coefficientcoulomb = np.nan @@ -23,5 +23,10 @@ self.effective_pressure_limit = 0 - self.setdefaultparameters() + nargs = len(args) + + if nargs == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') # }}} @@ -77,6 +82,8 @@ def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 7, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficientcoulomb', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2) @@ -85,5 +92,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double') if self.coupling == 1: - WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) elif self.coupling == 2: raise ValueError('not implemented yet') Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb2.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictioncoulomb2.py (revision 28030) +++ /issm/trunk-jpl/src/m/classes/frictioncoulomb2.py (revision 28030) @@ -0,0 +1,92 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class frictioncoulomb2(object): + """frictioncoulomb2 class definition + + Usage: + frictioncoulomb2 = frictioncoulomb2() + """ + + def __init__(self): # {{{ + self.coefficient = np.nan + self.coefficientcoulomb = np.nan + self.p = np.nan + self.q = np.nan + self.coupling = 0 + self.effective_pressure = np.nan + self.effective_pressure_limit = 0 + + self.setdefaultparameters() + # }}} + + def __repr__(self): # {{{ + s = 'Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n' + s += 'coefficientcoulomb^2 * Neff), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n' + s += '{}\n'.format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]")) + s += '{}\n'.format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]")) + s += '{}\n'.format(fielddisplay(self, "p", "p exponent")) + s += '{}\n'.format(fielddisplay(self, "q", "q exponent")) + s += '{}\n'.format(fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(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', 'Neff 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.coefficient = project3d(md, 'vector', self.coefficient, 'type', 'node', 'layer', 1) + self.coefficientcoulomb = project3d(md, 'vector', self.coefficientcoulomb, '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 == 1: + self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1) + elif self.coupling == 2: + raise ValueError('not implemented yet') + elif self.coupling > 2: + raise ValueError('not supported yet') + 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.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'friction.coefficientcoulomb', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'friction.q', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'friction.p', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2]) + md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0) + if self.coupling == 1: + md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + elif self.coupling == 2: + raise ValueError('not implemented yet') + elif self.coupling > 2: + raise ValueError('not supported yet') + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 7, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficientcoulomb', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2) + WriteData(fid, prefix, 'object', self, 'fieldname', 'q', '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 == 1: + 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 == 2: + raise ValueError('not implemented yet') + elif self.coupling > 2: + raise ValueError('not supported yet') + # }}} Index: /issm/trunk-jpl/src/m/classes/frictionregcoulomb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionregcoulomb.m (revision 28029) +++ /issm/trunk-jpl/src/m/classes/frictionregcoulomb.m (revision 28030) @@ -39,14 +39,14 @@ end % }}} function disp(self) % {{{ - %See Joughin et al. 2019 (equivalent form by Matt Trevers, poster at AGU 2022) https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL082526 - disp('Regularized Couloub friction law (Joughin et al., 2019) parameters:'); - disp(' Regularized Couloub friction law reads:'); + %See Joughin et al. 2019 (equivalent form by Matt Trevers, poster at AGU 2022) https://agupubs.onlinelåibrary.wiley.com/doi/full/10.1029/2019GL082526 + disp('Regularized Coulomb friction law (Joughin et al., 2019) parameters:'); + disp(' Regularized Coulomb friction law reads:'); disp(' C^2 |u|^(1/m) '); disp(' tau_b = - ____________________________'); - disp(' (|u|/u0 + 1)^(1/m) '); + disp(' (|u|/u0 + 1)^(1/m) '); disp(' '); fielddisplay(self,'C','friction coefficient [SI]'); fielddisplay(self,'m','m exponent (set to m=3 in original paper)'); - fielddisplay(self,'u0','velocity controling plastic limit'); + fielddisplay(self,'u0','velocity controlling plastic limit'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -54,5 +54,5 @@ WriteData(fid,prefix,'name','md.friction.law','data',14,'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','C','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); WriteData(fid,prefix,'class','friction','object',self,'fieldname','u0','format','Double','scale',1/yts); WriteData(fid,prefix,'class','friction','object',self,'fieldname','m','format','DoubleMat','mattype',2); Index: /issm/trunk-jpl/src/m/classes/frictionregcoulomb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionregcoulomb.py (revision 28030) +++ /issm/trunk-jpl/src/m/classes/frictionregcoulomb.py (revision 28030) @@ -0,0 +1,73 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class frictionregcoulomb(*args): + """frictionregcoulomb class definition + + Usage: + frictionregcoulomb = frictionregcoulomb() + """ + + def __init__(self, *args): # {{{ + self.C = np.nan + self.u0 = 0 + self.m = np.nan + + nargs = len(args) + + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + # See ./frictionregcoulomb.m + raise Exception('not implemented') + else: + raise Exception('constructor not supported') + # }}} + + def __repr__(self): # {{{ + # See Joughin et al. 2019 (equivalent form by Matt Trevers, poster at AGU 2022) https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL082526 + s = 'Regularized Coulomb friction law (Joughin et al., 2019) parameters:\n' + s += ' Regularized Coulomb friction law reads:\n' + s += ' C^2 |u|^(1/m) \n' + s += ' tau_b = - ____________________________\n' + s += ' (|u|/u0 + 1)^(1/m) \n' + s += '\n' + s += '{}\n'.format(fielddisplay(self, 'C', 'friction coefficient [SI]')) + s += '{}\n'.format(fielddisplay(self, 'm', 'm exponent (set to m = 3 in original paper)')) + s += '{}\n'.format(fielddisplay(self, 'u0', 'velocity controlling plastic limit')) + return s + # }}} + + def setdefaultparameters(self): # {{{ + self.u0 = 1000 + return self + # }}} + + def extrude(self, md): # {{{ + self.C = project3d(md, 'vector', self.C, '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.u0', 'NaN', 1, 'Inf', 1, '>', 0., 'numel', 1) + md = checkfield(md, 'fieldname', 'friction.m', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofelements, 1]); + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + + WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 14, 'format', 'Integer') + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'C', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'u0', 'format', 'Double', 'scale', 1 / yts) + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'm', 'format', 'DoubleMat', 'mattype', 2) + # }}} Index: /issm/trunk-jpl/src/m/classes/frictionregcoulomb2.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionregcoulomb2.m (revision 28029) +++ /issm/trunk-jpl/src/m/classes/frictionregcoulomb2.m (revision 28030) @@ -2,11 +2,11 @@ % % Usage: -% frictionregcoulomb=frictionregcoulomb(); +% frictionregcoulomb2=frictionregcoulomb2(); classdef frictionregcoulomb2 properties (SetAccess=public) - C = NaN; - K = NaN; - m = NaN; + C = NaN; + K = NaN; + m = NaN; effective_pressure_limit = 0; end @@ -43,13 +43,13 @@ function disp(self) % {{{ %See Zoet and Iverson 2020 or Choi et al., 2022 - disp('Regularized Couloub friction law 2 parameters:'); - disp(' Regularized Couloub friction law reads:'); + disp('Regularized Coulomb friction law 2 parameters:'); + disp(' Regularized Coulomb friction law reads:'); disp(' C N |u|^(1/m) '); disp(' tau_b = - ____________________________'); - disp(' (|u| + (K*N)^m)^(1/m) '); + disp(' (|u| + (K*N)^m)^(1/m) '); disp(' '); fielddisplay(self,'C','friction coefficient [SI]'); fielddisplay(self,'m','m exponent'); - fielddisplay(self,'K','(K*N)^m to be velocity controling plastic limit'); + fielddisplay(self,'K','(K*N)^m to be velocity controlling plastic limit'); 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 % }}} @@ -58,6 +58,6 @@ WriteData(fid,prefix,'name','md.friction.law','data',15,'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','K','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','C','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','K','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',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/frictionregcoulomb2.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionregcoulomb2.py (revision 28030) +++ /issm/trunk-jpl/src/m/classes/frictionregcoulomb2.py (revision 28030) @@ -0,0 +1,77 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class frictionregcoulomb2(*args): + """frictionregcoulomb2 class definition + + Usage: + frictionregcoulomb2 = frictionregcoulomb2() + """ + + def __init__(self, *args): # {{{ + self.C = np.nan + self.K = np.nan + self.m = np.nan + self.effective_pressure_limit = 0 + + nargs = len(args) + + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + # See ./frictionregcoulomb2.m + raise Exception('not implemented') + else: + raise Exception('constructor not supported') + # }}} + + def __repr__(self): # {{{ + # See Zoet and Iverson 2020 or Choi et al., 2022 + s = 'Regularized Coulomb friction law 2 parameters:\n' + s += ' Regularized Coulomb friction law reads:\n' + s += ' C N |u|^(1/m) \n' + s += ' tau_b = - ____________________________\n' + s += ' (|u| + (K*N)^m)^(1/m) \n' + s += '\n' + s += '{}\n'.format(fielddisplay(self, 'C', 'friction coefficient [SI]')) + s += '{}\n'.format(fielddisplay(self, 'm', 'm exponent')) + s += '{}\n'.format(fielddisplay(self, 'K', '(K * N) ^ m to be velocity controlling plastic limit')) + s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff 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.m = project3d(md, 'vector', self.m, 'type', 'element') + self.K = project3d(md, 'vector', self.K, 'type', 'node') + 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.K', 'NaN', 1, 'Inf', 1, '>', 0.) + md = checkfield(md, 'fieldname', 'friction.m', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofelements, 1]); + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + + WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 15, 'format', 'Integer') + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'C', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'K', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', 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') + # }}}