Changeset 18732

Timestamp:

11/04/14 14:13:02 (10 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added new friction law based on Tpmp, request from martin.rueckamp@…

Location:

issm/trunk-jpl/src

Files:

• c/analyses/EnthalpyAnalysis.cpp (modified) (2 diffs)
• c/analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• c/classes/Loads/Friction.cpp (modified) (2 diffs)
• c/classes/Loads/Friction.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (7 diffs)
• m/classes/frictiontemp.m (added)
• m/classes/groundingline.m (modified) (1 diff)
• m/classes/initialization.m (modified) (2 diffs)
• m/classes/initialization.py (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/FrictionGammaEnum.m (added)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -26 +26 @@
         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 +100 @@
                         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:

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -130 +130 @@
         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 +266 @@
                         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:

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -48 +48 @@
                 case 2:
                         GetAlpha2Weertman(palpha2,gauss);
-          case 3:
+                        break;
+                case 3:
                         GetAlpha2Hydro(palpha2,gauss);
+                        break;
+                case 4:
+                        GetAlpha2Temp(palpha2,gauss);
                         break;
           default:
@@ -219 +223 @@
         *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){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Loads/Friction.h

@@ -33 +33 @@
                 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);
 };

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -95 +95 @@
         FrictionCEnum,
         FrictionLawEnum,
+        FrictionGammaEnum,
         GeometryHydrostaticRatioEnum,
         HydrologyModelEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -103 +103 @@
                 case FrictionCEnum : return "FrictionC";
                 case FrictionLawEnum : return "FrictionLaw";
+                case FrictionGammaEnum : return "FrictionGamma";
                 case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
                 case HydrologyModelEnum : return "HydrologyModel";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -103 +103 @@
               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 +137 @@
               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 +260 @@
               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 +383 @@
               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 +506 @@
               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 +629 @@
               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 +752 @@
               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;

issm/trunk-jpl/src/m/classes/groundingline.m

@@ -9 +9 @@
         end
         methods
-         function createxml(obj,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- groundingline -->');
-            
-            % Convergence criteria          
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Grounding line migration parameters">','<section name="groundingline" />');  
-                        
-            % migration (SoftMigration, AggressiveMigration, or None)drop-down
-            fprintf(fid,'%s\n%s\n%s\n%s\n','<parameter key ="migration" type="alternative" optional="false">','     <section name="groundingline" />','     <help> type of grounding line migration: "SoftMigration","AggressiveMigration" or "None" </help>');
-            fprintf(fid,'%s\n','       <option value="SoftMigration" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="AggressiveMigration" type="string" default="false"> </option>');
-            fprintf(fid, '%s\n%s\n','       <option value="None" type="string" default="false"></option>','</parameter>');
+                function createxml(obj,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- groundingline -->');
 
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                        % Convergence criteria          
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Grounding line migration parameters">','<section name="groundingline" />');  
+
+                        % migration (SoftMigration, AggressiveMigration, or None)drop-down
+                        fprintf(fid,'%s\n%s\n%s\n%s\n','<parameter key ="migration" type="alternative" optional="false">','     <section name="groundingline" />','     <help> type of grounding line migration: "SoftMigration","AggressiveMigration" or "None" </help>');
+                        fprintf(fid,'%s\n','       <option value="SoftMigration" type="string" default="true"> </option>');
+                        fprintf(fid,'%s\n','       <option value="AggressiveMigration" type="string" default="false"> </option>');
+                        fprintf(fid, '%s\n%s\n','       <option value="None" type="string" default="false"></option>','</parameter>');
+
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
                 function obj = groundingline(varargin) % {{{
                         switch nargin

issm/trunk-jpl/src/m/classes/initialization.m

@@ -19 +19 @@
         end
         methods
-        function createxml(obj,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- initialization -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Initial field values">','<section name="initialization" />');                    
-            
-                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vx" type="',class(obj.vx),'" default="',obj.vx,'">','     <section name="initialization" />','     <help> x component of velocity [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vy" type="',class(obj.vy),'" default="',obj.vy,'">','     <section name="initialization" />','     <help> y component of velocity [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vz" type="',class(obj.vz),'" default="',obj.vz,'">','     <section name="initialization" />','     <help> z component of velocity [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vel" type="',class(obj.vel),'" default="',obj.vel,'">','     <section name="initialization" />','     <help> velocity norm [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="pressure" type="',class(obj.pressure),'" default="',obj.pressure,'">','     <section name="initialization" />','     <help> pressure field [Pa] </help>','</parameter>');
-                 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="temperature" type="',class(obj.temperature),'" default="',obj.temperature,'">','     <section name="initialization" />','     <help> fraction of water in the ice </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="waterfraction" type="',class(obj.waterfraction),'" default="',obj.waterfraction,'">','     <section name="initialization" />','     <help> ice thickness [m] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="sediment_head" type="',class(obj.sediment_head),'" default="',obj.sediment_head,'">','     <section name="initialization" />','     <help> sediment water head of subglacial system [m] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="epl_head" type="',class(obj.epl_head),'" default="',obj.epl_head,'">','     <section name="initialization" />','     <help> epl water head of subglacial system [m] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="watercolumn" type="',class(obj.watercolumn),'" default="',obj.watercolumn,'">','     <section name="initialization" />','     <help> thickness of subglacial water [m] </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+                function createxml(obj,fid) % {{{
+                        fprintf(fid, '\n\n');
+                        fprintf(fid, '%s\n', '<!-- initialization -->');
+                        fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Initial field values">','<section name="initialization" />');                    
+
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vx" type="',class(obj.vx),'" default="',obj.vx,'">','     <section name="initialization" />','     <help> x component of velocity [m/yr] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vy" type="',class(obj.vy),'" default="',obj.vy,'">','     <section name="initialization" />','     <help> y component of velocity [m/yr] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vz" type="',class(obj.vz),'" default="',obj.vz,'">','     <section name="initialization" />','     <help> z component of velocity [m/yr] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vel" type="',class(obj.vel),'" default="',obj.vel,'">','     <section name="initialization" />','     <help> velocity norm [m/yr] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="pressure" type="',class(obj.pressure),'" default="',obj.pressure,'">','     <section name="initialization" />','     <help> pressure field [Pa] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="temperature" type="',class(obj.temperature),'" default="',obj.temperature,'">','     <section name="initialization" />','     <help> fraction of water in the ice </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="waterfraction" type="',class(obj.waterfraction),'" default="',obj.waterfraction,'">','     <section name="initialization" />','     <help> ice thickness [m] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="sediment_head" type="',class(obj.sediment_head),'" default="',obj.sediment_head,'">','     <section name="initialization" />','     <help> sediment water head of subglacial system [m] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="epl_head" type="',class(obj.epl_head),'" default="',obj.epl_head,'">','     <section name="initialization" />','     <help> epl water head of subglacial system [m] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="watercolumn" type="',class(obj.watercolumn),'" default="',obj.watercolumn,'">','     <section name="initialization" />','     <help> thickness of subglacial water [m] </help>','</parameter>');
+                        fprintf(fid,'%s\n%s\n','</frame>');
+                end % }}}
                 function obj = initialization(varargin) % {{{
                         switch nargin
@@ -73 +73 @@
                                 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)

issm/trunk-jpl/src/m/classes/initialization.py

@@ -69 +69 @@
                         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])

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -95 +95 @@
 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]