Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18769) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18770) @@ -54,4 +54,7 @@ case 4: GetAlpha2Temp(palpha2,gauss); + break; + case 5: + GetAlpha2WaterLayer(palpha2,gauss); break; default: @@ -250,20 +253,17 @@ *palpha2=alpha2; }/*}}}*/ -void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ - - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. - * FrictionGetAlphaComplement is used in control methods on drag, and it computes: - * alpha_complement= Neff ^r * vel ^s*/ - - if(this->law!=1)_error_("not supported"); +void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /*This routine calculates the basal friction coefficient + alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ /*diverse: */ IssmDouble r,s; - IssmDouble vx,vy,vz,vmag; - IssmDouble drag_p,drag_q; + IssmDouble drag_p, drag_q; IssmDouble Neff; + IssmDouble thickness,bed; + IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient; - IssmDouble bed,thickness; - IssmDouble alpha_complement; + IssmDouble alpha2; /*Recover parameters: */ @@ -285,4 +285,67 @@ 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"); + } + + /*Check to prevent dividing by zero if vmag==0*/ + if(vmag==0. && (s-1.)<0.) alpha2=0.; + else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.)); + _assert_(!xIsNan(alpha2)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ +void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ + + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. + * FrictionGetAlphaComplement is used in control methods on drag, and it computes: + * alpha_complement= Neff ^r * vel ^s*/ + + if(this->law!=1)_error_("not supported"); + + /*diverse: */ + IssmDouble r,s; + IssmDouble vx,vy,vz,vmag; + IssmDouble drag_p,drag_q; + IssmDouble Neff; + IssmDouble drag_coefficient; + IssmDouble bed,thickness; + IssmDouble alpha_complement; + + /*Recover parameters: */ + element->GetInputValue(&drag_p,FrictionPEnum); + element->GetInputValue(&drag_q,FrictionQEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + + //compute r and q coefficients: */ + r=drag_q/drag_p; + s=1./drag_p; + + //From bed and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*bed); + if(Neff<0)Neff=0; + //We need the velocity magnitude to evaluate the basal stress: switch(dim){ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18769) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18770) @@ -34,4 +34,5 @@ void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss); void GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss); }; Index: /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m (revision 18770) +++ /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m (revision 18770) @@ -0,0 +1,54 @@ +%FRICTIONWATERLAYER class definition +% +% Usage: +% frictionwaterlayer=frictionwaterlayer(); + +classdef frictiontemp + properties (SetAccess=public) + coefficient = NaN; + p = NaN; + q = NaN; + water_layer = NaN; + end + methods + function obj = frictionwaterlayer(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + obj=structtoobj(frictiontemp(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + %Early return + if ~ismember(StressbalanceAnalysisEnum(),analyses) & ~ismember(ThermalAnalysisEnum(),analyses), return; end + + md = checkfield(md,'fieldname','friction.coefficient','forcing',1,'NaN',1); + md = checkfield(md,'fieldname','friction.q','NaN',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','friction.p','NaN',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','thermal.spctemperature','forcing',1,'>=',0.); + + end % }}} + function disp(obj) % {{{ + disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(bed+water_layer), r=q/p and s=1/p)')); + fielddisplay(obj,'coefficient','frictiontemp coefficient [SI]'); + fielddisplay(obj,'p','p exponent'); + fielddisplay(obj,'q','q exponent'); + fielddisplay(obj,'water_layer','water thickness at the base of the ice (m)'); + end % }}} + function marshall(obj,md,fid) % {{{ + + WriteData(fid,'enum',FrictionLawEnum,'data',5,'format','Integer'); + WriteData(fid,'class','friction','object',obj,'fieldname','coefficient','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'class','friction','object',obj,'fieldname','p','format','DoubleMat','mattype',2); + WriteData(fid,'class','friction','object',obj,'fieldname','q','format','DoubleMat','mattype',2); + WriteData(fid,'class','friction','object',obj,'fieldname','water_layer','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + end % }}} + end +end