Changeset 17742

Timestamp:

04/15/14 17:28:59 (11 years ago)

Author:

bdef

Message:

NEW: adding effective pressure computation within the Hydro

Location:

issm/trunk-jpl/src

Files:

• c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (4 diffs)
• c/cores/hydrology_core.cpp (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) (6 diffs)
• m/classes/transient.m (modified) (1 diff)
• m/classes/transient.py (modified) (1 diff)
• m/enum/EffectivePressureEnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)

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

@@ -303 +303 @@
         Input* sed_trans_input   = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum);
         Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
-        Input* base_input         = basalelement->GetInput(BaseEnum);
+        Input* base_input        = basalelement->GetInput(BaseEnum);
         Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
         if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
@@ -411 +411 @@
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* values   = xNew<IssmDouble>(numnodes);
+        IssmDouble* pressure = xNew<IssmDouble>(numnodes);
         IssmDouble* residual = xNew<IssmDouble>(numnodes);
 
@@ -419 +420 @@
         }
 
-        /*If converged keep the residual in mind*/
-        element->GetInputValue(&converged,ConvergedEnum);
+        /*If converged keep the residual in mind, also compute effective pressure*/
+        basalelement->GetInputValue(&converged,ConvergedEnum);
 
         /*Get inputs*/
         if(converged){
                 IssmDouble penalty_factor,kmax,kappa,h_max;
-                element->FindParam(&kmax,HydrologySedimentKmaxEnum);
-                element->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
+                IssmDouble rho_freshwater,rho_ice;
+                IssmDouble g;
+                IssmDouble* thickness = xNew<IssmDouble>(numnodes);
+                IssmDouble* base      = xNew<IssmDouble>(numnodes);
+
+                basalelement->FindParam(&kmax,HydrologySedimentKmaxEnum);
+                basalelement->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
+                basalelement->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
+                basalelement->GetMaterialParameter(MaterialsRhoIceEnum);
+                basalelement->GetMaterialParameter(ConstantsGEnum);
+                
+                basalelement->GetInputListOnVertices(thickness,ThicknessEnum);
+                basalelement->GetInputListOnVertices(base,BaseEnum);
 
                 kappa=kmax*pow(10.,penalty_factor);
 
                 for(int i=0;i<numnodes;i++){
+
                         GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->GetNode(i));
-                        if(values[i]>h_max) residual[i] = kappa*(values[i]-h_max);
-                        else                residual[i] = 0.;
-                }
+                        if(values[i]>h_max) {
+                                residual[i] = kappa*(values[i]-h_max);
+                                pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(values[i]-base[i]));
+                        }
+                        else{
+                                residual[i] = 0.;
+                                pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(h_max-base[i]));
+                        }
+                }
+                xDelete<IssmDouble>(thickness);
+                xDelete<IssmDouble>(base);
         }
 
@@ -440 +461 @@
         element->AddBasalInput(SedimentHeadEnum,values,P1Enum);
         element->AddBasalInput(SedimentHeadResidualEnum,residual,P1Enum);
+        element->AddBasalInput(EffectivePressureEnum,pressure,P1Enum);
 
         /*Free ressources:*/
         xDelete<IssmDouble>(values);
         xDelete<IssmDouble>(residual);
+        xDelete<IssmDouble>(pressure);
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

issm/trunk-jpl/src/c/cores/hydrology_core.cpp

@@ -64 +64 @@
                         if(VerboseSolution()) _printf0_("   saving results \n");
                         if(isefficientlayer){
-                                int outputs[8] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum};
-                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],8);
+                                int outputs[9] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum,EffectivePressureEnum};
+                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],9);
                         }
                         else{
-                                int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};
-                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+                                int outputs[3] = {SedimentHeadEnum,SedimentHeadResidualEnum,EffectivePressureEnum};
+                                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
                         }
                         /*unload results*/

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

@@ -98 +98 @@
         SedimentHeadOldEnum,
         SedimentHeadResidualEnum,
+        EffectivePressureEnum,
         EplHeadEnum,
         EplHeadOldEnum,

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

@@ -106 +106 @@
                 case SedimentHeadOldEnum : return "SedimentHeadOld";
                 case SedimentHeadResidualEnum : return "SedimentHeadResidual";
+                case EffectivePressureEnum : return "EffectivePressure";
                 case EplHeadEnum : return "EplHead";
                 case EplHeadOldEnum : return "EplHeadOld";

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

@@ -106 +106 @@
               else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
               else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
+              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
               else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
               else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
@@ -136 +137 @@
               else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
               else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
-              else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
+              if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
+              else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
               else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
               else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
@@ -259 +260 @@
               else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
               else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
-              else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
+              if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
+              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
               else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum;
               else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
@@ -382 +383 @@
               else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
               else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
-              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
+              if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
+              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
               else if (strcmp(name,"Loads")==0) return LoadsEnum;
               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
@@ -505 +506 @@
               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
               else if (strcmp(name,"Vx")==0) return VxEnum;
-              else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
+              if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
+              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
               else if (strcmp(name,"Vy")==0) return VyEnum;
               else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
@@ -628 +629 @@
               else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
               else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum;
-              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
+              if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
+              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
               else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
               else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;

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

@@ -54 +54 @@
 
                         %default output
-                        obj.requested_outputs={'default'};
+                        obj.requested_outputs={'none'};
                 end % }}}
                 function obj = setdefaultparameters(obj) % {{{

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

@@ -47 +47 @@
                         return []
 
+        #}}}
+        def setallnullparameters(self): # {{{
+                
+                #Nothing done
+                self.ismasstransport = False
+                self.isstressbalance = False
+                self.isthermal       = False
+                self.isgroundingline = False
+                self.isgia           = False
+                self.isdamageevolution = False
+                self.islevelset      = False
+                self.ishydrology     = False
+
+                #default output
+                self.requested_outputs=['none']
+                return self
         #}}}
         def setdefaultparameters(self): # {{{

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

@@ -98 +98 @@
 def SedimentHeadOldEnum(): return StringToEnum("SedimentHeadOld")[0]
 def SedimentHeadResidualEnum(): return StringToEnum("SedimentHeadResidual")[0]
+def EffectivePressureEnum(): return StringToEnum("EffectivePressure")[0]
 def EplHeadEnum(): return StringToEnum("EplHead")[0]
 def EplHeadOldEnum(): return StringToEnum("EplHeadOld")[0]