Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 18772) @@ -107,4 +107,10 @@ iomodel->FetchDataToInput(elements,PressureEnum); iomodel->FetchDataToInput(elements,TemperatureEnum); + break; + case 5: + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum); break; default: Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18772) @@ -274,4 +274,10 @@ iomodel->FetchDataToInput(elements,TemperatureEnum); break; + case 5: + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum); + break; default: _error_("not supported"); Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 18772) @@ -88,4 +88,10 @@ iomodel->FetchDataToInput(elements,FrictionCEnum); iomodel->FetchDataToInput(elements,FrictionMEnum); + break; + case 5: + iomodel->FetchDataToInput(elements,FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,FrictionPEnum); + iomodel->FetchDataToInput(elements,FrictionQEnum); + iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum); break; default: Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18772) @@ -264,6 +264,70 @@ IssmDouble thickness,bed; IssmDouble vx,vy,vz,vmag; + IssmDouble drag_coefficient,water_layer; + IssmDouble alpha2; + + /*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); + element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum); + 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+water_layer)); + 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 alpha2; + IssmDouble bed,thickness; + IssmDouble alpha_complement; /*Recover parameters: */ @@ -285,67 +349,4 @@ 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/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18771) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18772) @@ -96,4 +96,5 @@ FrictionLawEnum, FrictionGammaEnum, + FrictionWaterLayerEnum, GeometryHydrostaticRatioEnum, HydrologyModelEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18772) @@ -104,4 +104,5 @@ case FrictionLawEnum : return "FrictionLaw"; case FrictionGammaEnum : return "FrictionGamma"; + case FrictionWaterLayerEnum : return "FrictionWaterLayer"; case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio"; case HydrologyModelEnum : return "HydrologyModel"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18771) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18772) @@ -104,4 +104,5 @@ else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum; + else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; else if (strcmp(name,"HydrologyModel")==0) return HydrologyModelEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum; else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum; - else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum; + if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum; + else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum; else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"MeshY")==0) return MeshYEnum; else if (strcmp(name,"MeshZ")==0) return MeshZEnum; - else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"DomainType")==0) return DomainTypeEnum; + if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum; + else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; - else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; + if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; - else if (strcmp(name,"StringParam")==0) return StringParamEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Seg")==0) return SegEnum; + if (strcmp(name,"StringParam")==0) return StringParamEnum; + else if (strcmp(name,"Seg")==0) return SegEnum; else if (strcmp(name,"SegInput")==0) return SegInputEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"P2")==0) return P2Enum; else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; - else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"P2xP1")==0) return P2xP1Enum; + if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; + else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; else if (strcmp(name,"P1xP2")==0) return P1xP2Enum; else if (strcmp(name,"P1xP3")==0) return P1xP3Enum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; - else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; + if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum; + else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; Index: /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m (revision 18771) +++ /issm/trunk-jpl/src/m/classes/frictionwaterlayer.m (revision 18772) @@ -4,5 +4,5 @@ % frictionwaterlayer=frictionwaterlayer(); -classdef frictiontemp +classdef frictionwaterlayer properties (SetAccess=public) coefficient = NaN; @@ -17,5 +17,5 @@ obj=setdefaultparameters(obj); case 1 - obj=structtoobj(frictiontemp(),varargin{1}); + obj=structtoobj(frictionwaterlayer(),varargin{1}); otherwise error('constructor not supported'); Index: /issm/trunk-jpl/test/NightlyRun/test438.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test438.m (revision 18772) +++ /issm/trunk-jpl/test/NightlyRun/test438.m (revision 18772) @@ -0,0 +1,37 @@ +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'../Exp/SquareShelf.exp',''); +md=parameterize(md,'../Par/SquareSheetShelf.par'); +md=setflowequation(md,'SSA','all'); +md.friction=frictionwaterlayer(md.friction); +md.friction.water_layer=zeros(md.mesh.numberofvertices,2); +md.friction.water_layer(:,2)=1; +md.friction.water_layer(md.mesh.numberofvertices+1,:)=[1 2]; +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,TransientSolutionEnum()); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test439.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test439.m (revision 18772) +++ /issm/trunk-jpl/test/NightlyRun/test439.m (revision 18772) @@ -0,0 +1,38 @@ +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'../Exp/SquareShelf.exp',''); +md=parameterize(md,'../Par/SquareSheetShelf.par'); +md=extrude(md,4,1); +md=setflowequation(md,'HO','all'); +md.friction=frictionwaterlayer(md.friction); +md.friction.water_layer=zeros(md.mesh.numberofvertices,2); +md.friction.water_layer(:,2)=1; +md.friction.water_layer(md.mesh.numberofvertices+1,:)=[1 2]; +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,TransientSolutionEnum()); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'}; +field_tolerances={1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + };