Friction Class Reference

#include <Friction.h>

Public Member Functions
	Friction ()
	Friction (Element *element_in, int dim_in)
	~Friction ()
void	Echo (void)
void	GetAlphaComplement (IssmDouble *alpha_complement, Gauss *gauss)
void	GetAlphaHydroComplement (IssmDouble *alpha_complement, Gauss *gauss)
void	GetAlphaTempComplement (IssmDouble *alpha_complement, Gauss *gauss)
void	GetAlphaViscousComplement (IssmDouble *alpha_complement, Gauss *gauss)
void	GetAlpha2 (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Coulomb (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Hydro (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Josh (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Shakti (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Temp (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Viscous (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2WaterLayer (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Weertman (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2WeertmanTemp (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2PISM (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Schoof (IssmDouble *palpha2, Gauss *gauss)
void	GetAlpha2Tsai (IssmDouble *palpha2, Gauss *gauss)
IssmDouble	EffectivePressure (Gauss *gauss)
IssmDouble	VelMag (Gauss *gauss)

Data Fields
Element *	element
int	dim
int	law

Detailed Description

Definition at line 13 of file Friction.h.

Constructor & Destructor Documentation

Friction() [1/2]

Friction::Friction

(

)

Definition at line 18 of file Friction.cpp.

                  {/*{{{*/
   this->element=NULL;
   this->dim=0;
   this->law=0;
 
}

Friction() [2/2]

Friction::Friction	(	Element *	element_in,
		int	dim_in
	)

Definition at line 25 of file Friction.cpp.

                                                {/*{{{*/
 
   this->element=element_in;
   this->dim=dim_in;
   element_in->FindParam(&this->law,FrictionLawEnum);
}

~Friction()

Friction::~Friction

(

)

Definition at line 32 of file Friction.cpp.

                   {/*{{{*/
}

Member Function Documentation

Echo()

void Friction::Echo

(

void

)

Definition at line 37 of file Friction.cpp.

                       {/*{{{*/
   _printf_("Friction:\n");
   _printf_("   dim: " << this->dim<< "\n");
}

GetAlphaComplement()

void Friction::GetAlphaComplement	(	IssmDouble *	alpha_complement,
		Gauss *	gauss
	)

Definition at line 42 of file Friction.cpp.

                                                                            {/*{{{*/
 
   switch(this->law){
      case 1:
         GetAlphaViscousComplement(palpha_complement,gauss);
         break;
      case 3:
         GetAlphaHydroComplement(palpha_complement,gauss);
         break;
      case 4:
         GetAlphaTempComplement(palpha_complement,gauss);
         break;
     default:
         _error_("not supported");
   }
 
 
   /*Checks*/
   _assert_(!xIsNan<IssmDouble>(*palpha_complement));
   _assert_(!xIsInf<IssmDouble>(*palpha_complement));
 
}/*}}}*/

GetAlphaHydroComplement()

void Friction::GetAlphaHydroComplement	(	IssmDouble *	alpha_complement,
		Gauss *	gauss
	)

Definition at line 64 of file Friction.cpp.

                                                                                 {/*{{{*/
 
   /*diverse: */
   IssmDouble  q_exp;
   IssmDouble  C_param;
   IssmDouble  As;
   IssmDouble  n;
   IssmDouble  alpha;
   IssmDouble  Chi,Gamma;
   IssmDouble  alpha_complement;
 
   /*Recover parameters: */
   element->GetInputValue(&q_exp,gauss,FrictionQEnum);
   element->GetInputValue(&C_param,gauss,FrictionCEnum);
   element->GetInputValue(&As,gauss,FrictionAsEnum);
   element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
 
   /*Get effective pressure and velocity magnitude*/
   IssmDouble Neff = EffectivePressure(gauss);
   IssmDouble vmag = VelMag(gauss);
 
   if (q_exp==1){
      alpha=1;
   }
   else{
      alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
   }
   Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
   Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
   /*Check to prevent dividing by zero if vmag==0*/
   if(vmag==0.) alpha_complement=0.;
   else  if(Neff==0.) alpha_complement=0.;
   else  alpha_complement=-(C_param*Neff/(n*vmag)) *
               pow(Gamma,((1.-n)/n)) *
               (Gamma/As - (alpha*q_exp*pow(Chi,q_exp-1.)* Gamma * Gamma/As));
 
   /*Assign output pointers:*/
   *palpha_complement=alpha_complement;
}

GetAlphaTempComplement()

void Friction::GetAlphaTempComplement	(	IssmDouble *	alpha_complement,
		Gauss *	gauss
	)

Definition at line 104 of file Friction.cpp.

                                                                                {/*{{{*/
   /*Here, we want to parameterize the friction as a function of temperature
    *
    * alpha2 = alpha2_viscous * 1/f(T)
    *
    * where f(T) = exp((T-Tpmp)/gamma)
    */
 
   /*Intermediaries: */
   IssmDouble  f,T,pressure,Tpmp,gamma;
   IssmDouble  alpha_complement;
 
   /*Get viscous part*/
   this->GetAlphaViscousComplement(&alpha_complement,gauss);
 
   /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
   element->GetInputValue(&T,gauss,TemperatureEnum);
   element->GetInputValue(&pressure,gauss,PressureEnum);
   Tpmp = element->TMeltingPoint(pressure);
 
   /*Compute scaling parameter*/
   element->parameters->FindParam(&gamma,FrictionGammaEnum);
   alpha_complement = alpha_complement/ (exp((T-Tpmp)/gamma)+1e-3);
 
   /*Assign output pointers:*/
   *palpha_complement=alpha_complement;
}/*}}}*/

GetAlphaViscousComplement()

void Friction::GetAlphaViscousComplement	(	IssmDouble *	alpha_complement,
		Gauss *	gauss
	)

Definition at line 131 of file Friction.cpp.

                                                                                   {/*{{{*/
 
   /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
    * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    * alpha_complement= Neff ^r * vel ^s*/
 
   /*diverse: */
   IssmDouble  r,s;
   IssmDouble  drag_p,drag_q;
   IssmDouble  drag_coefficient;
   IssmDouble  alpha_complement;
 
   /*Recover parameters: */
   element->GetInputValue(&drag_p,gauss,FrictionPEnum);
   element->GetInputValue(&drag_q,gauss,FrictionQEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 
   /*compute r and q coefficients: */
   r=drag_q/drag_p;
   s=1./drag_p;
 
   /*Get effective pressure*/
   IssmDouble Neff = EffectivePressure(gauss);
   IssmDouble vmag = VelMag(gauss);
 
   /*Check to prevent dividing by zero if vmag==0*/
   if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
   else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
 
   /*Assign output pointers:*/
   *palpha_complement=alpha_complement;
}

GetAlpha2()

void Friction::GetAlpha2	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 164 of file Friction.cpp.

                                                         {/*{{{*/
 
   switch(this->law){
      case 1:
         GetAlpha2Viscous(palpha2,gauss);
         break;
      case 2:
         GetAlpha2Weertman(palpha2,gauss);
         break;
      case 3:
         GetAlpha2Hydro(palpha2,gauss);
         break;
      case 4:
         GetAlpha2Temp(palpha2,gauss);
         break;
      case 5:
         GetAlpha2WaterLayer(palpha2,gauss);
         break;
      case 6:
         GetAlpha2WeertmanTemp(palpha2,gauss);
         break;
      case 7:
         GetAlpha2Coulomb(palpha2,gauss);
         break;
      case 8:
         GetAlpha2Shakti(palpha2,gauss);
         break;
      case 9:
         GetAlpha2Josh(palpha2,gauss);
         break;
      case 10:
         GetAlpha2PISM(palpha2,gauss);
         break;
      case 11:
         GetAlpha2Schoof(palpha2,gauss);
         break;
      case 12:
         GetAlpha2Tsai(palpha2,gauss);
         break;
     default:
         _error_("Friction law "<< this->law <<" not supported");
   }
 
   /*Checks*/
   _assert_(!xIsNan<IssmDouble>(*palpha2));
   _assert_(!xIsInf<IssmDouble>(*palpha2));
   _assert_(*palpha2>=0);
 
}/*}}}*/

GetAlpha2Coulomb()

void Friction::GetAlpha2Coulomb	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 213 of file Friction.cpp.

                                                                {/*{{{*/
 
   /*This routine calculates the basal friction coefficient
     alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
 
   /*diverse: */
   IssmDouble  r,s;
   IssmDouble  drag_p, drag_q;
   IssmDouble  thickness,base,floatation_thickness,sealevel;
   IssmDouble  drag_coefficient,drag_coefficient_coulomb;
   IssmDouble  alpha2,alpha2_coulomb;
 
   /*Recover parameters: */
   element->GetInputValue(&drag_p,gauss,FrictionPEnum);
   element->GetInputValue(&drag_q,gauss,FrictionQEnum);
   element->GetInputValue(&thickness, gauss,ThicknessEnum);
   element->GetInputValue(&base, gauss,BaseEnum);
   element->GetInputValue(&sealevel, gauss,SealevelEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
   element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble gravity   = element->FindParam(ConstantsGEnum);
 
   //compute r and q coefficients: */
   r=drag_q/drag_p;
   s=1./drag_p;
 
   /*Get effective pressure*/
   IssmDouble Neff = EffectivePressure(gauss);
   IssmDouble vmag = VelMag(gauss);
 
   /*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.));
   }
 
   floatation_thickness=0;
   if(base<0) floatation_thickness=-(rho_water/rho_ice)*base;
   if(vmag==0.){
      alpha2_coulomb=0.;
   }
   else{
      alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
   }
 
   if(alpha2_coulomb<alpha2) alpha2=alpha2_coulomb;
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Hydro()

void Friction::GetAlpha2Hydro	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 267 of file Friction.cpp.

                                                              {/*{{{*/
 
   /*This routine calculates the basal friction coefficient
      Based on Gagliardini 2007, needs a good effective pressure computation
      Not tested so far so use at your own risks
     alpha2= NeffC[Chi/(1+alpha*Chi^q)]^(1/n)*|vel|^(-1)  with
       Chi=|vel|/(C^n*Neff^n*As)
       alpha=(q-1)^(q-1)/q^q */
 
   /*diverse: */
   IssmDouble  q_exp;
   IssmDouble  C_param;
   IssmDouble  As;
   IssmDouble  n;
   IssmDouble  alpha;
   IssmDouble  Chi,Gamma;
   IssmDouble  alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&q_exp,gauss,FrictionQEnum);
   element->GetInputValue(&C_param,gauss,FrictionCEnum);
   element->GetInputValue(&As,gauss,FrictionAsEnum);
   element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
 
   /*Get effective pressure*/
   IssmDouble Neff = EffectivePressure(gauss);
   IssmDouble vmag = VelMag(gauss);
 
   //compute alpha and Chi coefficients: */
   if (q_exp==1){
      alpha=1;
   }
   else{
      alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
   }
   Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
   Gamma=(Chi/(1. + alpha * pow(Chi,q_exp)));
   /*Check to prevent dividing by zero if vmag==0*/
   if(vmag==0.) alpha2=0.;
   else  if (Neff==0) alpha2=0.0;
   else  alpha2=Neff * C_param * pow(Gamma,1./n) * pow(vmag,-1);
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Josh()

void Friction::GetAlpha2Josh	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 372 of file Friction.cpp.

                                                             {/*{{{*/
   /*Here, we want to parameterize the friction as a function of temperature
    *
    * alpha2 = alpha2_viscous * 1/f(T)
    *
    * where f(T) = exp((T-Tpmp)/gamma)
    */
 
   /*Intermediaries: */
   IssmDouble  T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient;
   IssmDouble  alpha2,time,gamma,ref,alp_new,alphascaled;
   const IssmDouble yts = 365*24*3600.;
 
   /*Get viscous part*/
   this->GetAlpha2Viscous(&alpha2,gauss);
 
   /*Get delta Refs*/
   element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
   /*New*/
   /*element->GetInputValue(&deltaTrefsfc,gauss,FrictionSurfaceTemperatureEnum);
    *    element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum);
    *       */
 
   /*Compute delta T*/
   element->GetInputValue(&T,gauss,TemperatureEnum);
   element->GetInputValue(&pressure,gauss,PressureEnum);
   Tpmp = element->TMeltingPoint(pressure);
   deltaT = T-Tpmp;
 
 
   /*Compute gamma*/
   element->parameters->FindParam(&time,TimeEnum);
   element->parameters->FindParam(&gamma,FrictionGammaEnum);
 
   ref = exp(deltaTref/gamma);
   alp_new = ref/exp(deltaT/gamma);
 
   alphascaled = sqrt(alp_new)*drag_coefficient;
   if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient);
 
   alp_new=alp_new*alpha2;
 
   /*Assign output pointers:*/
   *palpha2=alp_new;
}/*}}}*/

GetAlpha2Shakti()

void Friction::GetAlpha2Shakti	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 312 of file Friction.cpp.

                                                               {/*{{{*/
 
   /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/
 
   /*diverse: */
   IssmDouble  pressure_ice,pressure_water;
   IssmDouble  Neff;
   IssmDouble  drag_coefficient;
   IssmDouble  base,thickness,head,sealevel;
   IssmDouble  alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&thickness, gauss,ThicknessEnum);
   element->GetInputValue(&base, gauss,BaseEnum);
   element->GetInputValue(&head, gauss,HydrologyHeadEnum);
   element->GetInputValue(&sealevel, gauss,SealevelEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
   IssmDouble rho_water   = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble gravity     = element->FindParam(ConstantsGEnum);
 
   //From base and thickness, compute effective pressure when drag is viscous:
   pressure_ice   = rho_ice*gravity*thickness;
   pressure_water = rho_water*gravity*(head-base+sealevel);
   Neff=pressure_ice-pressure_water;
   if(Neff<0.) Neff=0.;
 
   alpha2=drag_coefficient*drag_coefficient*Neff;
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}

GetAlpha2Temp()

void Friction::GetAlpha2Temp	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 345 of file Friction.cpp.

                                                             {/*{{{*/
   /*Here, we want to parameterize the friction as a function of temperature
    *
    * alpha2 = alpha2_viscous * 1/f(T)
    *
    * where f(T) = exp((T-Tpmp)/gamma)
    */
 
   /*Intermediaries: */
   IssmDouble  f,T,pressure,Tpmp,gamma;
   IssmDouble  alpha2;
 
   /*Get viscous part*/
   this->GetAlpha2Viscous(&alpha2,gauss);
 
   /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
   element->GetInputValue(&T,gauss,TemperatureEnum);
   element->GetInputValue(&pressure,gauss,PressureEnum);
   Tpmp = element->TMeltingPoint(pressure);
 
   /*Compute scaling parameter*/
   element->parameters->FindParam(&gamma,FrictionGammaEnum);
   alpha2 = alpha2 / (exp((T-Tpmp)/gamma)+1e-3);
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Viscous()

void Friction::GetAlpha2Viscous	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 418 of file Friction.cpp.

                                                                {/*{{{*/
 
   /*This routine calculates the basal friction coefficient
     alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
 
   /*diverse: */
   IssmDouble  r,s;
   IssmDouble  drag_p, drag_q;
   IssmDouble  drag_coefficient;
   IssmDouble  alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&drag_p,gauss,FrictionPEnum);
   element->GetInputValue(&drag_q,gauss,FrictionQEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 
   /*compute r and q coefficients: */
   r=drag_q/drag_p;
   s=1./drag_p;
 
   /*Get effective pressure and basal velocity*/
   IssmDouble Neff = EffectivePressure(gauss);
   IssmDouble vmag = VelMag(gauss);
 
   /*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.));
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2WaterLayer()

void Friction::GetAlpha2WaterLayer	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 449 of file Friction.cpp.

                                                                   {/*{{{*/
 
   /*This routine calculates the basal friction coefficient
     alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
 
   /*diverse: */
   IssmDouble  r,s;
   IssmDouble  drag_p, drag_q;
   IssmDouble  Neff,F;
   IssmDouble  thickness,base,sealevel;
   IssmDouble  drag_coefficient,water_layer;
   IssmDouble  alpha2;
 
   /*Recover parameters: */
   element->parameters->FindParam(&F,FrictionFEnum);
   element->GetInputValue(&drag_p,gauss,FrictionPEnum);
   element->GetInputValue(&drag_q,gauss,FrictionQEnum);
   element->GetInputValue(&thickness, gauss,ThicknessEnum);
   element->GetInputValue(&base, gauss,BaseEnum);
   element->GetInputValue(&sealevel, gauss,SealevelEnum);
   element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
   element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
   IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
   IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble gravity     = element->FindParam(ConstantsGEnum);
 
   //compute r and q coefficients: */
   r=drag_q/drag_p;
   s=1./drag_p;
 
   //From base and thickness, compute effective pressure when drag is viscous:
   if(base>0) base=0;
   if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(base-sealevel);
   else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
   else _error_("negative water layer thickness");
   if(Neff<0) Neff=0;
 
   IssmDouble vmag = VelMag(gauss);
 
   /*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.));
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Weertman()

void Friction::GetAlpha2Weertman	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 495 of file Friction.cpp.

                                                                 {/*{{{*/
 
   /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
 
   /*diverse: */
   IssmDouble  C,m;
   IssmDouble  alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&C,gauss,FrictionCEnum);
   element->GetInputValue(&m,gauss,FrictionMEnum);
 
   /*Get velocity magnitude*/
   IssmDouble vmag = VelMag(gauss);
 
   /*Check to prevent dividing by zero if vmag==0*/
   if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
   else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2WeertmanTemp()

void Friction::GetAlpha2WeertmanTemp	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 517 of file Friction.cpp.

                                                                     {/*{{{*/
   /*Here, we want to parameterize the friction as a function of temperature
    *
    * alpha2 = alpha2_weertman * 1/f(T)
    *
    * where f(T) = exp((T-Tpmp)/gamma)
    */
 
   /*Intermediaries: */
   IssmDouble  f,T,pressure,Tpmp,gamma;
   IssmDouble  alpha2;
 
   /*Get viscous part*/
   this->GetAlpha2Weertman(&alpha2,gauss);
 
   /*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
   element->GetInputValue(&T,gauss,TemperatureEnum);
   element->GetInputValue(&pressure,gauss,PressureEnum);
   Tpmp = element->TMeltingPoint(pressure);
 
   /*Compute scaling parameter*/
   element->parameters->FindParam(&gamma,FrictionGammaEnum);
   alpha2 = alpha2 / exp((T-Tpmp)/gamma);
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2PISM()

void Friction::GetAlpha2PISM	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 544 of file Friction.cpp.

                                                             {/*{{{*/
   /*Here, we want to parameterize the friction using a pseudoplastic friction law,
    * computing the basal shear stress as
    *
    * alpha2 = tau_c (u_b/(abs(u_b)^(1-q)*u_0^q))
    *
    * The yield stress tau_c is a function of the effective pressure N
    * using a Mohr-Coloumb criterion, so that
    * tau_c = tan(phi)*N,
    * where phi is the till friction angle, representing sediment strength
    *
    * The effective pressure is given by Eq. (5) in Aschwanden et al. 2016:
    *
    * N = delta * P0 * 10^((e_0/Cc)(1-(W/Wmax)))
    *
    * W is calculated by a non-conserving hydrology model in HydrologyPismAnalysis.cpp
    *
    * see Aschwanden et al. 2016 and Bueler and Brown, 2009 for more details
    */
 
   /*compute ice overburden pressure P0*/
   IssmDouble thickness,base,P0;
   element->GetInputValue(&thickness, gauss,ThicknessEnum);
   element->GetInputValue(&base, gauss,BaseEnum);
   //element->GetInputValue(&sealevel, gauss,SealevelEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble gravity   = element->FindParam(ConstantsGEnum);
   P0 = gravity*rho_ice*thickness;
 
   /*Compute effective pressure*/
   IssmDouble  N,delta,W,Wmax,e0,Cc;
   element->parameters->FindParam(&delta,FrictionDeltaEnum);
   element->parameters->FindParam(&e0,FrictionVoidRatioEnum);
   element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum);
   element->GetInputValue(&W,gauss,WatercolumnEnum);
   element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum);
 
   /*Check that water column height is within 0 and upper bound, correct if needed*/
   if(W>Wmax) W=Wmax;
   if(W<0)    W=0.;
 
   N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
 
   /*Get till friction angles, defined by user [deg]*/
   IssmDouble phi;
   element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
 
   /*Convert till friction angle from user-defined deg to rad, which Matlab uses*/
   phi = phi*PI/180.;
 
   /*Compute yield stress following a Mohr-Colomb criterion*/
   IssmDouble tau_c = N*tan(phi);
 
   /*Compute basal speed*/
   IssmDouble ub;
   element->GetInputValue(&ub,gauss,VelEnum);
 
   /*now compute alpha^2*/
   IssmDouble u0,q;
   element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum);
   element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum);
   IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q));
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Schoof()

void Friction::GetAlpha2Schoof	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 610 of file Friction.cpp.

                                                               {/*{{{*/
 
   /*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,alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
   element->GetInputValue(&C,gauss,FrictionCEnum);
   element->GetInputValue(&m,gauss,FrictionMEnum);
 
   /*Get effective pressure and velocity magnitude*/
   IssmDouble N  = EffectivePressure(gauss);
   IssmDouble ub = VelMag(gauss);
 
   /*Compute alpha^2*/
   if(ub<1e-10){
      alpha2 = 0.;
   }
   else{
      alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
   }
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

GetAlpha2Tsai()

void Friction::GetAlpha2Tsai	(	IssmDouble *	palpha2,
		Gauss *	gauss
	)

Definition at line 643 of file Friction.cpp.

                                                             {/*{{{*/
 
   /*This routine calculates the basal friction coefficient
    *
    * alpha2= min(C |ub|^m , f N ) / |ub|
    *
    * */
 
   /*diverse: */
   IssmDouble  C,f,m,alpha2;
 
   /*Recover parameters: */
   element->GetInputValue(&f,gauss,FrictionfEnum);
   element->GetInputValue(&C,gauss,FrictionCEnum);
   element->GetInputValue(&m,gauss,FrictionMEnum);
 
   /*Get effective pressure and velocity magnitude*/
   IssmDouble N  = EffectivePressure(gauss);
   IssmDouble ub = VelMag(gauss);
 
   /*Compute alpha^2*/
   if(ub<1e-10){
      alpha2 = 0.;
   }
   else{
      alpha2= C*pow(ub,m);
 
      if(alpha2>f*N) alpha2 = f*N;
 
      alpha2 = alpha2/ub;
   }
 
   /*Assign output pointers:*/
   *palpha2=alpha2;
}/*}}}*/

EffectivePressure()

IssmDouble Friction::EffectivePressure

(

Gauss *

gauss

)

Definition at line 679 of file Friction.cpp.

                                                  {/*{{{*/
   /*Get effective pressure as a function of  flag */
 
   /*diverse: */
   int         coupled_flag;
   IssmDouble  thickness,base,sealevel;
   IssmDouble  p_ice,p_water;
   IssmDouble  Neff,Neff_limit;
 
   /*Recover parameters: */
   element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
   element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
 
   /*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
   switch(coupled_flag){
      case 0:{
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
         element->GetInputValue(&base, gauss,BaseEnum);
         element->GetInputValue(&sealevel, gauss,SealevelEnum);
         IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
         p_ice   = gravity*rho_ice*thickness;
         p_water = rho_water*gravity*(sealevel-base);
         Neff = p_ice - p_water;
      }
         break;
      case 1:{
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
         p_ice   = gravity*rho_ice*thickness;
         p_water = 0.;
         Neff = p_ice - p_water;
      }
         break;
      case 2:{
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
         element->GetInputValue(&base, gauss,BaseEnum);
         element->GetInputValue(&sealevel, gauss,SealevelEnum);
         IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
         p_ice   = gravity*rho_ice*thickness;
         p_water = max(0.,rho_water*gravity*(sealevel-base));
         Neff = p_ice - p_water;
      }
         break;
      case 3:{
         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
         p_ice   = gravity*rho_ice*thickness;
      }
         break;
      case 4:{
         element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble gravity   = element->FindParam(ConstantsGEnum);
         p_ice   = gravity*rho_ice*thickness;
      }
         break;
      default:
         _error_("not supported");
   }
 
   /*Make sure Neff is positive*/
   if(Neff<Neff_limit*p_ice) Neff=Neff_limit*p_ice;
 
   /*Return effective pressure*/
   return Neff;
 
}/*}}}*/

VelMag()

IssmDouble Friction::VelMag

(

Gauss *

gauss

)

Definition at line 754 of file Friction.cpp.

                                       {/*{{{*/
   /*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;
 
}/*}}}*/

Field Documentation

element

Element* Friction::element

Definition at line 16 of file Friction.h.

dim

int Friction::dim

Definition at line 17 of file Friction.h.

law

int Friction::law

Definition at line 18 of file Friction.h.

---

The documentation for this class was generated from the following files:

• src/c/classes/Loads/Friction.h
• src/c/classes/Loads/Friction.cpp