Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 18731)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 18732)
@@ -26,4 +26,9 @@
if(numoutputs)parameters->AddObject(new StringArrayParam(ThermalRequestedOutputsEnum,requestedoutputs,numoutputs));
iomodel->DeleteData(&requestedoutputs,numoutputs,ThermalRequestedOutputsEnum);
+
+ /*Deal with friction parameters*/
+ int frictionlaw;
+ iomodel->Constant(&frictionlaw,FrictionLawEnum);
+ if(frictionlaw==4) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
}/*}}}*/
void EnthalpyAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -95,4 +100,11 @@
iomodel->FetchDataToInput(elements,FrictionCEnum);
iomodel->FetchDataToInput(elements,FrictionMEnum);
+ break;
+ case 4:
+ iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
+ iomodel->FetchDataToInput(elements,FrictionPEnum);
+ iomodel->FetchDataToInput(elements,FrictionQEnum);
+ iomodel->FetchDataToInput(elements,PressureEnum);
+ iomodel->FetchDataToInput(elements,TemperatureEnum);
break;
default:
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18731)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18732)
@@ -130,4 +130,9 @@
iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
+ /*Deal with friction parameters*/
+ int frictionlaw;
+ iomodel->Constant(&frictionlaw,FrictionLawEnum);
+ if(frictionlaw==4) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
+
}/*}}}*/
void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -261,4 +266,11 @@
iomodel->FetchDataToInput(elements,FrictionCEnum);
iomodel->FetchDataToInput(elements,FrictionMEnum);
+ break;
+ case 4:
+ iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
+ iomodel->FetchDataToInput(elements,FrictionPEnum);
+ iomodel->FetchDataToInput(elements,FrictionQEnum);
+ iomodel->FetchDataToInput(elements,PressureEnum);
+ iomodel->FetchDataToInput(elements,TemperatureEnum);
break;
default:
Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18731)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18732)
@@ -48,6 +48,10 @@
case 2:
GetAlpha2Weertman(palpha2,gauss);
- case 3:
+ break;
+ case 3:
GetAlpha2Hydro(palpha2,gauss);
+ break;
+ case 4:
+ GetAlpha2Temp(palpha2,gauss);
break;
default:
@@ -219,4 +223,31 @@
*palpha2=alpha2;
}/*}}}*/
+void Friction::GetAlpha2Temp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+ /*Here, we want to parameterize the friction as a function of temperature
+ *
+ * alpha2 = alpha2_viscous * 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);
+
+ /*Assign output pointers:*/
+ *palpha2=alpha2;
+}/*}}}*/
void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18731)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18732)
@@ -33,4 +33,5 @@
void GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
+ void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
void GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
};
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18731)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18732)
@@ -95,4 +95,5 @@
FrictionCEnum,
FrictionLawEnum,
+ FrictionGammaEnum,
GeometryHydrostaticRatioEnum,
HydrologyModelEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18731)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18732)
@@ -103,4 +103,5 @@
case FrictionCEnum : return "FrictionC";
case FrictionLawEnum : return "FrictionLaw";
+ case FrictionGammaEnum : return "FrictionGamma";
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 18731)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18732)
@@ -103,4 +103,5 @@
else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum;
+ else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum;
else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
else if (strcmp(name,"HydrologyModel")==0) return HydrologyModelEnum;
@@ -136,9 +137,9 @@
else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
- else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
else stage=2;
}
if(stage==2){
- if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
+ if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
+ else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
@@ -259,9 +260,9 @@
else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
- else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
else stage=3;
}
if(stage==3){
- if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
+ if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
+ else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
@@ -382,9 +383,9 @@
else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
- else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
else stage=4;
}
if(stage==4){
- if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
+ if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
+ else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
@@ -505,9 +506,9 @@
else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
- else if (strcmp(name,"Air")==0) return AirEnum;
else stage=5;
}
if(stage==5){
- if (strcmp(name,"Ice")==0) return IceEnum;
+ if (strcmp(name,"Air")==0) return AirEnum;
+ else if (strcmp(name,"Ice")==0) return IceEnum;
else if (strcmp(name,"Melange")==0) return MelangeEnum;
else if (strcmp(name,"Water")==0) return WaterEnum;
@@ -628,9 +629,9 @@
else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
- else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
else stage=6;
}
if(stage==6){
- if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
+ if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
+ else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
@@ -751,9 +752,9 @@
else if (strcmp(name,"SeaiceMinConcentration")==0) return SeaiceMinConcentrationEnum;
else if (strcmp(name,"SeaiceMinThickness")==0) return SeaiceMinThicknessEnum;
- else if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
else stage=7;
}
if(stage==7){
- if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
+ if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
+ else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
else if (strcmp(name,"SeaiceCoriolisFactor")==0) return SeaiceCoriolisFactorEnum;
Index: /issm/trunk-jpl/src/m/classes/frictiontemp.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictiontemp.m (revision 18732)
+++ /issm/trunk-jpl/src/m/classes/frictiontemp.m (revision 18732)
@@ -0,0 +1,60 @@
+%FRICTIONTEMP class definition
+%
+% Usage:
+% frictiontemp=frictiontemp();
+
+classdef frictiontemp
+ properties (SetAccess=public)
+ gamma = 0;
+ coefficient = NaN;
+ p = NaN;
+ q = NaN;
+ end
+ methods
+ function obj = frictiontemp(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) % {{{
+
+ %By default gamma = 1;
+ obj.gamma = 1;
+
+
+ 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','friction.gamma','NaN',1,'numel',1,'>',0.);
+
+ %Check that temperature is provided
+ md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'size',[md.mesh.numberofvertices 1]);
+ 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, r=q/p and s=1/p)'));
+ fielddisplay(obj,'gamma','submelt sliding parameter f(T) = exp((T-Tpmp)/gamma)');
+ fielddisplay(obj,'coefficient','frictiontemp coefficient [SI]');
+ fielddisplay(obj,'p','p exponent');
+ fielddisplay(obj,'q','q exponent');
+ end % }}}
+ function marshall(obj,md,fid) % {{{
+
+ WriteData(fid,'enum',FrictionLawEnum,'data',4,'format','Integer');
+ WriteData(fid,'class','friction','object',obj,'fieldname','gamma','format','Double');
+ 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);
+ end % }}}
+ end
+end
Index: /issm/trunk-jpl/src/m/classes/groundingline.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/groundingline.m (revision 18731)
+++ /issm/trunk-jpl/src/m/classes/groundingline.m (revision 18732)
@@ -9,19 +9,19 @@
end
methods
- function createxml(obj,fid) % {{{
- fprintf(fid, '\n\n');
- fprintf(fid, '%s\n', '');
-
- % Convergence criteria
- fprintf(fid,'%s\n%s\n%s\n','','');
-
- % migration (SoftMigration, AggressiveMigration, or None)drop-down
- fprintf(fid,'%s\n%s\n%s\n%s\n','',' ',' type of grounding line migration: "SoftMigration","AggressiveMigration" or "None" ');
- fprintf(fid,'%s\n',' ');
- fprintf(fid,'%s\n',' ');
- fprintf(fid, '%s\n%s\n',' ','');
+ function createxml(obj,fid) % {{{
+ fprintf(fid, '\n\n');
+ fprintf(fid, '%s\n', '');
- fprintf(fid,'%s\n%s\n','');
- end % }}}
+ % Convergence criteria
+ fprintf(fid,'%s\n%s\n%s\n','','');
+
+ % migration (SoftMigration, AggressiveMigration, or None)drop-down
+ fprintf(fid,'%s\n%s\n%s\n%s\n','',' ',' type of grounding line migration: "SoftMigration","AggressiveMigration" or "None" ');
+ fprintf(fid,'%s\n',' ');
+ fprintf(fid,'%s\n',' ');
+ fprintf(fid, '%s\n%s\n',' ','');
+
+ fprintf(fid,'%s\n%s\n','');
+ end % }}}
function obj = groundingline(varargin) % {{{
switch nargin
Index: /issm/trunk-jpl/src/m/classes/initialization.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.m (revision 18731)
+++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 18732)
@@ -19,21 +19,21 @@
end
methods
- function createxml(obj,fid) % {{{
- fprintf(fid, '\n\n');
- fprintf(fid, '%s\n', '');
- fprintf(fid,'%s\n%s\n%s\n','','');
-
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' x component of velocity [m/yr] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' y component of velocity [m/yr] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' z component of velocity [m/yr] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' velocity norm [m/yr] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' pressure field [Pa] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' fraction of water in the ice ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' ice thickness [m] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' sediment water head of subglacial system [m] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' epl water head of subglacial system [m] ','');
- fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','',' ',' thickness of subglacial water [m] ','');
- fprintf(fid,'%s\n%s\n','');
- end % }}}
+ function createxml(obj,fid) % {{{
+ fprintf(fid, '\n\n');
+ fprintf(fid, '%s\n', '');
+ fprintf(fid,'%s\n%s\n%s\n','','');
+
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' x component of velocity [m/yr] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' y component of velocity [m/yr] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' z component of velocity [m/yr] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' velocity norm [m/yr] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' pressure field [Pa] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' fraction of water in the ice ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' ice thickness [m] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' sediment water head of subglacial system [m] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','',' ',' epl water head of subglacial system [m] ','');
+ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','',' ',' thickness of subglacial water [m] ','');
+ fprintf(fid,'%s\n%s\n','');
+ end % }}}
function obj = initialization(varargin) % {{{
switch nargin
@@ -73,4 +73,5 @@
end
md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
+ md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'size',[md.mesh.numberofvertices 1]);
end
if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy)
Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py (revision 18731)
+++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 18732)
@@ -69,4 +69,5 @@
md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices])
md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
+ md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'size',[md.mesh.numberofvertices])
if md.mesh.dimension()==3:
md = checkfield(md,'fieldname','initialization.vz','NaN',1,'size',[md.mesh.numberofvertices])
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18731)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18732)
@@ -95,4 +95,5 @@
def FrictionCEnum(): return StringToEnum("FrictionC")[0]
def FrictionLawEnum(): return StringToEnum("FrictionLaw")[0]
+def FrictionGammaEnum(): return StringToEnum("FrictionGamma")[0]
def GeometryHydrostaticRatioEnum(): return StringToEnum("GeometryHydrostaticRatio")[0]
def HydrologyModelEnum(): return StringToEnum("HydrologyModel")[0]
Index: /issm/trunk-jpl/src/m/enum/FrictionGammaEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/FrictionGammaEnum.m (revision 18732)
+++ /issm/trunk-jpl/src/m/enum/FrictionGammaEnum.m (revision 18732)
@@ -0,0 +1,11 @@
+function macro=FrictionGammaEnum()
+%FRICTIONGAMMAENUM - Enum of FrictionGamma
+%
+% WARNING: DO NOT MODIFY THIS FILE
+% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+% Please read src/c/shared/Enum/README for more information
+%
+% Usage:
+% macro=FrictionGammaEnum()
+
+macro=StringToEnum('FrictionGamma');