Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23838) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23839) @@ -870,4 +870,7 @@ iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.); break; + case 11: + iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum); + iomodel->FetchDataToInput(elements,"md.friction.Cmax",FrictionCmaxEnum); default: _error_("friction law "<< frictionlaw <<" not supported"); @@ -931,13 +934,41 @@ int frictionlaw; iomodel->FindConstant(&frictionlaw,"md.friction.law"); - if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); - if(frictionlaw==3 || frictionlaw==4 || frictionlaw==1 || frictionlaw==7) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); - if(frictionlaw==5) parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum)); - if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); - if(frictionlaw==10){ - parameters->AddObject(iomodel->CopyConstantObject("md.friction.pseudoplasticity_exponent",FrictionPseudoplasticityExponentEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.friction.threshold_speed",FrictionThresholdSpeedEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.friction.delta",FrictionDeltaEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.friction.void_ratio",FrictionVoidRatioEnum)); + switch(frictionlaw){ + case 1: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + break; + case 2: + break; + case 3: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + break; + case 4: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + break; + case 5: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum)); + break; + case 6: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); + break; + case 7: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + break; + case 8: + break; + case 9: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); + break; + case 10: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.pseudoplasticity_exponent",FrictionPseudoplasticityExponentEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.threshold_speed",FrictionThresholdSpeedEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.delta",FrictionDeltaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.void_ratio",FrictionVoidRatioEnum)); + break; + case 11: + parameters->AddObject(iomodel->CopyConstantObject("md.friction.m",FrictionMEnum)); + break; + default: _error_("Friction law "<GetInputValue(&n,gauss,MaterialsRheologyNEnum); - /*Get effective pressure*/ + /*Get effective pressure and velocity magnitude*/ IssmDouble Neff = EffectivePressure(gauss); - - //We need the velocity magnitude to evaluate the basal stress: - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } - // vmag=100./(3600.*24.*365.); + IssmDouble vmag = VelMag(gauss); if (q_exp==1){ @@ -155,5 +133,4 @@ /*diverse: */ IssmDouble r,s; - IssmDouble vx,vy,vz,vmag; IssmDouble drag_p,drag_q; IssmDouble drag_coefficient; @@ -171,25 +148,5 @@ /*Get effective pressure*/ IssmDouble Neff = EffectivePressure(gauss); - - //We need the velocity magnitude to evaluate the basal stress: - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } + IssmDouble vmag = VelMag(gauss); /*Check to prevent dividing by zero if vmag==0*/ @@ -237,4 +194,7 @@ GetAlpha2PISM(palpha2,gauss); break; + case 11: + GetAlpha2Schoof(palpha2,gauss); + break; default: _error_("Friction law "<< this->law <<" not supported"); @@ -251,5 +211,4 @@ IssmDouble drag_p, drag_q; IssmDouble thickness,base,floatation_thickness,sealevel; - IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient,drag_coefficient_coulomb; IssmDouble alpha2,alpha2_coulomb; @@ -273,25 +232,5 @@ /*Get effective pressure*/ IssmDouble Neff = EffectivePressure(gauss); - - /*Compute velocity magnitude*/ - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } + IssmDouble vmag = VelMag(gauss); /*Check to prevent dividing by zero if vmag==0*/ @@ -325,9 +264,6 @@ IssmDouble As; IssmDouble n; - IssmDouble alpha; IssmDouble Chi,Gamma; - - IssmDouble vx,vy,vz,vmag; IssmDouble alpha2; @@ -340,27 +276,6 @@ /*Get effective pressure*/ IssmDouble Neff = EffectivePressure(gauss); - - /*Compute velocity magnitude*/ - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } - - // vmag=100./(3600.*24.*365.); + IssmDouble vmag = VelMag(gauss); + //compute alpha and Chi coefficients: */ if (q_exp==1){ @@ -496,5 +411,4 @@ IssmDouble r,s; IssmDouble drag_p, drag_q; - IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient; IssmDouble alpha2; @@ -509,27 +423,7 @@ s=1./drag_p; - /*Get effective pressure*/ + /*Get effective pressure and basal velocity*/ IssmDouble Neff = EffectivePressure(gauss); - - /*Compute velocity magnitude*/ - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } + IssmDouble vmag = VelMag(gauss); /*Check to prevent dividing by zero if vmag==0*/ @@ -551,5 +445,4 @@ IssmDouble Neff,F; IssmDouble thickness,base,sealevel; - IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient,water_layer; IssmDouble alpha2; @@ -579,23 +472,5 @@ if(Neff<0) Neff=0; - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } + IssmDouble vmag = VelMag(gauss); /*Check to prevent dividing by zero if vmag==0*/ @@ -613,5 +488,4 @@ /*diverse: */ IssmDouble C,m; - IssmDouble vx,vy,vz,vmag; IssmDouble alpha2; @@ -620,23 +494,6 @@ element->GetInputValue(&m,FrictionMEnum); - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } + /*Get velocity magnitude*/ + IssmDouble vmag = VelMag(gauss); /*Check to prevent dividing by zero if vmag==0*/ @@ -741,4 +598,36 @@ _assert_(!xIsNan(alpha2)); _assert_(!xIsInf(alpha2)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ +void Friction::GetAlpha2Schoof(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /*This routine calculates the basal friction coefficient + * + * C |u_b|^(m-1) + * alpha2= __________________________ + * (1+(C/(Cmax N))^1/m |u_b| )^m + * + * */ + + /*diverse: */ + IssmDouble C,Cmax,m; + + /*Recover parameters: */ + element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum); + element->GetInputValue(&C,gauss,FrictionCEnum); + element->GetInputValue(&m,FrictionMEnum); + + /*Get effective pressure and velocity magnitude*/ + IssmDouble N = EffectivePressure(gauss); + IssmDouble ub = VelMag(gauss); + + /*Compute alpha^2*/ + IssmDouble alpha2= (C*pow(ub,m-1.)) / pow(1.+ pow(C/(Cmax*N),1./m)*ub,m); + + /*Checks*/ + _assert_(!xIsNan(alpha2)); + _assert_(alpha2>=0); /*Assign output pointers:*/ @@ -810,2 +699,32 @@ }/*}}}*/ +IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/ + /*Get effective pressure as a function of flag */ + + + /*diverse*/ + IssmDouble vx,vy,vz,vmag; + + switch(dim){ + case 1: + element->GetInputValue(&vx,gauss,VxEnum); + vmag=sqrt(vx*vx); + break; + case 2: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + vmag=sqrt(vx*vx+vy*vy); + break; + case 3: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + element->GetInputValue(&vz,gauss,VzEnum); + vmag=sqrt(vx*vx+vy*vy+vz*vz); + break; + default: + _error_("not supported"); + } + + return vmag; + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 23838) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 23839) @@ -41,6 +41,8 @@ void GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2PISM(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Schoof(IssmDouble* palpha2,Gauss* gauss); IssmDouble EffectivePressure(Gauss* gauss); + IssmDouble VelMag(Gauss* gauss); }; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23838) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23839) @@ -502,4 +502,5 @@ FrictionAsEnum, FrictionCEnum, + FrictionCmaxEnum, FrictionCoefficientcoulombEnum, FrictionCoefficientEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23838) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23839) @@ -508,4 +508,5 @@ case FrictionAsEnum : return "FrictionAs"; case FrictionCEnum : return "FrictionC"; + case FrictionCmaxEnum : return "FrictionCmax"; case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb"; case FrictionCoefficientEnum : return "FrictionCoefficient"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23838) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23839) @@ -520,4 +520,5 @@ else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum; else if (strcmp(name,"FrictionC")==0) return FrictionCEnum; + else if (strcmp(name,"FrictionCmax")==0) return FrictionCmaxEnum; else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum; else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum; else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; - else if (strcmp(name,"SmbC")==0) return SmbCEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum; + if (strcmp(name,"SmbC")==0) return SmbCEnum; + else if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum; else if (strcmp(name,"SmbDailyrainfall")==0) return SmbDailyrainfallEnum; else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; - else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; + if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; + else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum; - else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; + else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; - else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"None")==0) return NoneEnum; + if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; + else if (strcmp(name,"None")==0) return NoneEnum; else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; else if (strcmp(name,"P2")==0) return P2Enum; - else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"P2xP4")==0) return P2xP4Enum; + if (strcmp(name,"P2xP1")==0) return P2xP1Enum; + else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"MeshZ")==0) return MeshZEnum; else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; - else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; + if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; + else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; Index: /issm/trunk-jpl/src/m/classes/frictionschoof.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 23839) +++ /issm/trunk-jpl/src/m/classes/frictionschoof.m (revision 23839) @@ -0,0 +1,62 @@ +%FRICTIONSCHOOF class definition +% +% Usage: +% frictionschoof=frictionschoof(); + +classdef frictionschoof + properties (SetAccess=public) + C = NaN; + Cmax = NaN; + m = 0.; + end + methods + function self = frictionschoof(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = extrude(self,md) % {{{ + md.friction.C = project3d(md,'vector',md.friction.C,'type','node','layer',1); + md.friction.Cmax = project3d(md,'vector',md.friction.Cmax,'type','node','layer',1); + end % }}} + function self = setdefaultparameters(self) % {{{ + + %default m is 1/n = 1/3 + self.m = 1./3.; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %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.Cmax','timeseries',1,'NaN',1,'Inf',1,'>',0.); + md = checkfield(md,'fieldname','friction.m','NaN',1,'Inf',1,'numel',1,'>',0.); + end % }}} + function disp(self) % {{{ + %See Brondex et al. 2019 https://www.the-cryosphere.net/13/177/2019/ + disp('Schoof sliding law parameters:'); + disp(' Schoof''s sliding law reads:'); + disp(' C |u_b|^(m-1) '); + disp(' tau_b = - _____________________________ u_b '); + disp(' (1+(C/(Cmax N))^1/m |u_b| )^m '); + disp(' '); + fielddisplay(self,'C','friction coefficient [SI]'); + 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)'); + end % }}} + function 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','Double'); + + + end % }}} + end +end