Changeset 19744 for issm/trunk-jpl/src

Timestamp:

11/17/15 16:10:39 (9 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added discrete moulins to hydrology model

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/analyses/HydrologySommersAnalysis.cpp (modified) (2 diffs)
• c/classes/Loads/Moulin.cpp (added)
• c/classes/Loads/Moulin.h (added)
• c/classes/classes.h (modified) (1 diff)
• c/cores/hydrology_core.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (8 diffs)
• m/classes/hydrologysommers.m (modified) (4 diffs)
• m/consistency/ismodelselfconsistent.m (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/HydrologyMoulinInputEnum.m (added)
• m/enum/MoulinEnum.m (added)
• m/plot/plot_transient_movie.m (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -69 +69 @@
                                         ./classes/Loads/Penpair.cpp\
                                         ./classes/Loads/Pengrid.cpp\
+                                        ./classes/Loads/Moulin.cpp\
                                         ./classes/Loads/Numericalflux.cpp\
                                         ./classes/matrix/ElementMatrix.cpp\

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

@@ -22 +22 @@
 }/*}}}*/
 void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-        /*No loads*/
+
+        /*Fetch parameters: */
+        int  hydrology_model;
+        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
+
+        /*Now, do we really want Sommers?*/
+        if(hydrology_model!=HydrologysommersEnum) return;
+
+        //create penalties for nodes: no node can have water above the max
+        CreateSingleNodeToElementConnectivity(iomodel);
+        for(int i=0;i<iomodel->numberofvertices;i++){
+                if (iomodel->domaintype!=Domain3DEnum){
+                        /*keep only this partition's nodes:*/
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
+                        }
+                }
+                else if(reCast<int>(iomodel->Data(MeshVertexonbaseEnum)[i])){
+                        if(iomodel->my_vertices[i]){
+                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
+                        }       
+                }
+        }
+        iomodel->DeleteData(1,MeshVertexonbaseEnum);
+
 }/*}}}*/
 void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -72 +96 @@
         iomodel->FetchDataToInput(elements,HydrologyGapHeightEnum);
         iomodel->FetchDataToInput(elements,HydrologyEnglacialInputEnum);
+        iomodel->FetchDataToInput(elements,HydrologyMoulinInputEnum);
         iomodel->FetchDataToInput(elements,HydrologyBumpSpacingEnum);
         iomodel->FetchDataToInput(elements,HydrologyBumpHeightEnum);

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

@@ -37 +37 @@
 #include "./Loads/Penpair.h"
 #include "./Loads/Pengrid.h"
+#include "./Loads/Moulin.h"
 
 /*Elements: */

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

@@ -60 +60 @@
                 
                 /*Proceed now to heads computations*/
-                if(VerboseSolution()) _printf0_("   computing water head\n");
                 solutionsequence_hydro_nonlinear(femmodel);
 

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

@@ -162 +162 @@
         HydrologyBumpHeightEnum,
         HydrologyEnglacialInputEnum,
+        HydrologyMoulinInputEnum,
         HydrologyReynoldsEnum,
         HydrologySpcheadEnum,
@@ -579 +580 @@
         NumericalfluxTypeEnum,
         ParamEnum,
+        MoulinEnum,
         PengridEnum,
         PenpairEnum,

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

@@ -168 +168 @@
                 case HydrologyBumpHeightEnum : return "HydrologyBumpHeight";
                 case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput";
+                case HydrologyMoulinInputEnum : return "HydrologyMoulinInput";
                 case HydrologyReynoldsEnum : return "HydrologyReynolds";
                 case HydrologySpcheadEnum : return "HydrologySpchead";
@@ -571 +572 @@
                 case NumericalfluxTypeEnum : return "NumericalfluxType";
                 case ParamEnum : return "Param";
+                case MoulinEnum : return "Moulin";
                 case PengridEnum : return "Pengrid";
                 case PenpairEnum : return "Penpair";

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

@@ -171 +171 @@
               else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum;
               else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
+              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
               else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
               else if (strcmp(name,"HydrologySpchead")==0) return HydrologySpcheadEnum;
@@ -259 +260 @@
               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
               else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
-              else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
+              if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
+              else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
               else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
@@ -382 +383 @@
               else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
               else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
-              else if (strcmp(name,"SmbP")==0) return SmbPEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"SmbSwf")==0) return SmbSwfEnum;
+              if (strcmp(name,"SmbP")==0) return SmbPEnum;
+              else if (strcmp(name,"SmbSwf")==0) return SmbSwfEnum;
               else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
               else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum;
@@ -505 +506 @@
               else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
               else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
-              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
+              if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
+              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
               else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
               else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
@@ -583 +584 @@
               else if (strcmp(name,"NumericalfluxType")==0) return NumericalfluxTypeEnum;
               else if (strcmp(name,"Param")==0) return ParamEnum;
+              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
               else if (strcmp(name,"Pengrid")==0) return PengridEnum;
               else if (strcmp(name,"Penpair")==0) return PenpairEnum;
@@ -627 +629 @@
               else if (strcmp(name,"Free")==0) return FreeEnum;
               else if (strcmp(name,"Open")==0) return OpenEnum;
-              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
-              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
+              if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
+              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
+              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
               else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
               else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
@@ -750 +752 @@
               else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
               else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
-              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
-              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
+              if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
+              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
+              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
               else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
               else if (strcmp(name,"J")==0) return JEnum;
@@ -873 +875 @@
               else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
               else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
-              else if (strcmp(name,"MisfitLocal")==0) return MisfitLocalEnum;
-              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
+              if (strcmp(name,"MisfitLocal")==0) return MisfitLocalEnum;
+              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
               else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
               else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;

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

@@ -11 +11 @@
                 bump_height     = NaN;
                 englacial_input = NaN;
+                moulin_input    = NaN;
                 reynolds        = NaN;
                 spchead         = NaN;
@@ -40 +41 @@
                         md = checkfield(md,'fieldname','hydrology.bump_spacing','>',0,'size',[md.mesh.numberofelements 1],'NaN',1);
                         md = checkfield(md,'fieldname','hydrology.bump_height','>=',0,'size',[md.mesh.numberofelements 1],'NaN',1);
-                        md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'size',[md.mesh.numberofvertices 1],'NaN',1);
+                        md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'NaN',1,'timeseries',1);
+                        md = checkfield(md,'fieldname','hydrology.moulin_input','>=',0,'NaN',1,'timeseries',1);
                         md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements 1],'NaN',1);
                         md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);
@@ -51 +53 @@
                         fielddisplay(self,'bump_height','characteristic bedrock bump height (m)');
                         fielddisplay(self,'englacial_input','liquid water input from englacial to subglacial system (m/yr)');
+                        fielddisplay(self,'moulin_input','liquid water input from moulins (at the vertices) to subglacial system (m^3/s)');
                         fielddisplay(self,'reynolds','Reynolds'' number');
                         fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)');
@@ -63 +66 @@
                         WriteData(fid,'object',self,'class','hydrology','fieldname','bump_spacing','format','DoubleMat','mattype',2);
                         WriteData(fid,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',2);
-                        WriteData(fid,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts);
+                        WriteData(fid,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1);
+                        WriteData(fid,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
                         WriteData(fid,'object',self,'class','hydrology','fieldname','reynolds','format','DoubleMat','mattype',2);
                         WriteData(fid,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1);

issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m

@@ -80 +80 @@
 
         case TransientSolutionEnum(),
-                analyses=[StressbalanceAnalysisEnum();StressbalanceVerticalAnalysisEnum();StressbalanceSIAAnalysisEnum();L2ProjectionBaseAnalysisEnum();ThermalAnalysisEnum();MeltingAnalysisEnum();EnthalpyAnalysisEnum();MasstransportAnalysisEnum()];
+                analyses=[StressbalanceAnalysisEnum();StressbalanceVerticalAnalysisEnum();StressbalanceSIAAnalysisEnum();L2ProjectionBaseAnalysisEnum();ThermalAnalysisEnum();MeltingAnalysisEnum();EnthalpyAnalysisEnum();MasstransportAnalysisEnum();HydrologySommersAnalysisEnum];
 
         case FlaimSolutionEnum(),

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

@@ -160 +160 @@
 def HydrologyBumpHeightEnum(): return StringToEnum("HydrologyBumpHeight")[0]
 def HydrologyEnglacialInputEnum(): return StringToEnum("HydrologyEnglacialInput")[0]
+def HydrologyMoulinInputEnum(): return StringToEnum("HydrologyMoulinInput")[0]
 def HydrologyReynoldsEnum(): return StringToEnum("HydrologyReynolds")[0]
 def HydrologySpcheadEnum(): return StringToEnum("HydrologySpchead")[0]
@@ -563 +564 @@
 def NumericalfluxTypeEnum(): return StringToEnum("NumericalfluxType")[0]
 def ParamEnum(): return StringToEnum("Param")[0]
+def MoulinEnum(): return StringToEnum("Moulin")[0]
 def PengridEnum(): return StringToEnum("Pengrid")[0]
 def PenpairEnum(): return StringToEnum("Penpair")[0]

issm/trunk-jpl/src/m/plot/plot_transient_movie.m

@@ -62 +62 @@
 
                         clf;
-                        titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+                        titlestring=[field ' at time ' num2str(results(i).time) ' year'];
                         plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
                         apply_options_movie(md,options,titlestring);