Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 21739) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 21740) @@ -174,7 +174,11 @@ switch(frictionlaw){ case 1: + iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); + if (FrictionCoupling==1){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + } break; case 2: @@ -187,5 +191,5 @@ iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if (FrictionCoupling==0){ + if (FrictionCoupling==1){ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); } @@ -245,5 +249,5 @@ iomodel->FindConstant(&frictionlaw,"md.friction.law"); if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); - if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21739) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21740) @@ -781,7 +781,11 @@ switch(frictionlaw){ case 1: + iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); + if(FrictionCoupling==1){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + } break; case 2: @@ -794,5 +798,5 @@ iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if(FrictionCoupling==0){ + if(FrictionCoupling==1){ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); } @@ -893,5 +897,5 @@ iomodel->FindConstant(&frictionlaw,"md.friction.law"); if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); - if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + if(frictionlaw==3 || frictionlaw==1) 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)); Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 21739) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 21740) @@ -157,7 +157,11 @@ switch(frictionlaw){ case 1: + iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling"); iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); + if (FrictionCoupling==1){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + } break; case 2: @@ -170,6 +174,6 @@ iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); - if (FrictionCoupling==0){ - iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); + if (FrictionCoupling==1){ + iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum); } break; @@ -220,5 +224,6 @@ iomodel->FindConstant(&frictionlaw,"md.friction.law"); if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); - if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum)); + }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21739) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21740) @@ -70,4 +70,5 @@ IssmDouble vx,vy,vz,vmag; IssmDouble alpha_complement; + IssmDouble base,sealevel,thickness; /*Recover parameters: */ @@ -77,11 +78,25 @@ element->GetInputValue(&As,gauss,FrictionAsEnum); element->GetInputValue(&n,gauss,MaterialsRheologyNEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&base, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); + + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); - if (CoupledFlag==1){ - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); - } - else{ - element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); + switch(CoupledFlag){ + case 0: + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); + break; + case 1: + element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); + break; + case 2: + element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + break; + default: + _error_("not supported"); } @@ -159,9 +174,10 @@ void Friction::GetAlphaViscousComplement(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. + /* 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: */ + int CoupledFlag; IssmDouble r,s; IssmDouble vx,vy,vz,vmag; @@ -169,103 +185,6 @@ IssmDouble Neff; IssmDouble drag_coefficient; - IssmDouble bed,thickness,sealevel; + IssmDouble base,thickness,sealevel; 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(&sealevel, gauss,SealevelEnum); - 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-sealevel)); - if(Neff<0)Neff=0; - - //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"); - } - - /*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));_assert_(!xIsNan(alpha_complement)); - - /*Assign output pointers:*/ - *palpha_complement=alpha_complement; -} -/*}}}*/ -void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - - 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: - GetAlpha2Sommers(palpha2,gauss); - break; - case 9: - GetAlpha2Josh(palpha2,gauss); - break; - default: - _error_("Friction law "<< this->law <<" not supported"); - } - -}/*}}}*/ -void Friction::GetAlpha2Coulomb(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*base, r=q/p and s=1/p**/ - - /*diverse: */ - IssmDouble r,s; - IssmDouble drag_p, drag_q; - IssmDouble Neff; - IssmDouble thickness,base,bed,floatation_thickness,sealevel; - IssmDouble vx,vy,vz,vmag; - IssmDouble drag_coefficient,drag_coefficient_coulomb; - IssmDouble alpha2,alpha2_coulomb; /*Recover parameters: */ @@ -275,5 +194,113 @@ element->GetInputValue(&base, gauss,BaseEnum); element->GetInputValue(&sealevel, gauss,SealevelEnum); - element->GetInputValue(&bed, gauss,BedEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); + //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, or get Neff from forcing or coupled to hydrologymodel: + switch(CoupledFlag){ + case 0: + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); + break; + case 1: + element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); + break; + case 2: + element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + break; + default: + _error_("not supported"); + } + if(Neff<0)Neff=0; + + //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"); + } + + /*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));_assert_(!xIsNan(alpha_complement)); + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; +} +/*}}}*/ +void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + 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: + GetAlpha2Sommers(palpha2,gauss); + break; + case 9: + GetAlpha2Josh(palpha2,gauss); + break; + default: + _error_("Friction law "<< this->law <<" not supported"); + } + +}/*}}}*/ +void Friction::GetAlpha2Coulomb(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*base, r=q/p and s=1/p**/ + + /*diverse: */ + IssmDouble r,s; + IssmDouble drag_p, drag_q; + IssmDouble Neff; + IssmDouble thickness,base,floatation_thickness,sealevel; + IssmDouble vx,vy,vz,vmag; + IssmDouble drag_coefficient,drag_coefficient_coulomb; + IssmDouble alpha2,alpha2_coulomb; + + /*Recover parameters: */ + element->GetInputValue(&drag_p,FrictionPEnum); + element->GetInputValue(&drag_q,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); @@ -315,5 +342,5 @@ floatation_thickness=0; - if(bed<0) floatation_thickness=-rho_water/rho_ice*bed; + 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_water*gravity*(thickness-floatation_thickness)/vmag; @@ -349,4 +376,5 @@ IssmDouble vx,vy,vz,vmag; IssmDouble alpha2; + IssmDouble base,thickness,sealevel; /*Recover parameters: */ @@ -355,13 +383,26 @@ element->GetInputValue(&As,gauss,FrictionAsEnum); element->GetInputValue(&n,gauss,MaterialsRheologyNEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&base, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); - element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); - if (CoupledFlag==1){ - element->GetInputValue(&Neff,gauss,EffectivePressureEnum); - } - else{ - element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); - } - + switch(CoupledFlag){ + case 0: + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); + break; + case 1: + element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); + break; + case 2: + element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + break; + default: + _error_("not supported"); + } + if(Neff<0)Neff=0; @@ -407,5 +448,5 @@ void Friction::GetAlpha2Sommers(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-bed)*/ + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/ /*diverse: */ @@ -413,10 +454,10 @@ IssmDouble Neff; IssmDouble drag_coefficient; - IssmDouble bed,thickness,head,sealevel; + IssmDouble base,thickness,head,sealevel; IssmDouble alpha2; /*Recover parameters: */ element->GetInputValue(&thickness, gauss,ThicknessEnum); - element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&base, gauss,BaseEnum); element->GetInputValue(&head, gauss,HydrologyHeadEnum); element->GetInputValue(&sealevel, gauss,SealevelEnum); @@ -426,7 +467,7 @@ IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); - //From bed and thickness, compute effective pressure when drag is viscous: + //From base and thickness, compute effective pressure when drag is viscous: pressure_ice = rho_ice*gravity*thickness; - pressure_water = rho_water*gravity*(head-bed+sealevel); + pressure_water = rho_water*gravity*(head-base+sealevel); Neff=pressure_ice-pressure_water; if(Neff<0.) Neff=0.; @@ -511,7 +552,8 @@ /*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**/ + 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: */ + int CoupledFlag; IssmDouble r,s; IssmDouble drag_p, drag_q; @@ -532,11 +574,23 @@ IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); - + element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum); //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: - Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); + //From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing: + switch(CoupledFlag){ + case 0: + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); + break; + case 1: + element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum); + break; + case 2: + element->GetInputValue(&Neff,gauss,EffectivePressureEnum); + break; + default: + _error_("not supported"); + } if(Neff<0)Neff=0; @@ -572,5 +626,5 @@ /*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**/ + 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: */ @@ -578,5 +632,5 @@ IssmDouble drag_p, drag_q; IssmDouble Neff,F; - IssmDouble thickness,bed,sealevel; + IssmDouble thickness,base,sealevel; IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient,water_layer; @@ -588,5 +642,5 @@ element->GetInputValue(&drag_q,FrictionQEnum); element->GetInputValue(&thickness, gauss,ThicknessEnum); - element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&base, gauss,BaseEnum); element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); @@ -600,7 +654,7 @@ s=1./drag_p; - //From bed and thickness, compute effective pressure when drag is viscous: - if(bed>0) bed=0; - if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(bed-sealevel); + //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");