Changeset 23093

Timestamp:

08/16/18 13:32:45 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing PISM hydro model

Location:

issm/trunk-jpl

Files:

• src/c/analyses/HydrologyPismAnalysis.cpp (modified) (4 diffs)
• src/c/classes/FemModel.cpp (modified) (1 diff)
• src/c/classes/Materials/Matpar.cpp (modified) (1 diff)
• src/c/cores/hydrology_core.cpp (modified) (2 diffs)
• src/c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• src/c/toolkits/petsc/objects/PetscSolver.cpp (modified) (1 diff)
• src/m/classes/hydrologypism.m (modified) (1 diff)
• src/m/classes/hydrologyshreve.m (modified) (1 diff)
• test/Archives/Archive610.arch (modified) ( previous)
• test/NightlyRun/test610.m (added)

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

@@ -32 +32 @@
 
         /*Add input to elements*/
+        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
+        iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
+        iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
         iomodel->FetchDataToInput(elements,"md.hydrology.drainage_rate",HydrologyDrainageRateEnum);
-        iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn",WatercolumnEnum,0.);
+        iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
 }/*}}}*/
 void HydrologyPismAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
@@ -46 +49 @@
         if(hydrology_model!=HydrologypismEnum) return;
         parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
+
+        /*Requested outputs*/
+        iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
+        parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
+        if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
+        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
 
         /*Nothing else to add for now*/
@@ -98 +107 @@
         /*Retrieve all inputs and parameters*/
         element->FindParam(&dt,TimesteppingTimeStepEnum);
+        IssmDouble  rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
 
         /*Get water column and drainage rate*/
@@ -109 +120 @@
 
         /*Add water*/
-        for(int i=0;i<numvertices;i++) watercolumn[i] += (meltingrate[i]-drainagerate[i])*dt;
+        for(int i=0;i<numvertices;i++) watercolumn[i] += (meltingrate[i]/rho_ice*rho_water-drainagerate[i])*dt;
 
         /* Divide by connectivity, add degree of channelization as an input */

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -771 +771 @@
                                 analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum;
                                 analyses_temp[numanalyses++]=HydrologyShaktiAnalysisEnum;
+                                analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum;
                                 analyses_temp[numanalyses++]=HydrologyDCInefficientAnalysisEnum;
                                 analyses_temp[numanalyses++]=HydrologyDCEfficientAnalysisEnum;

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -150 +150 @@
                         }
                         else if(hydrology_model==HydrologyshaktiEnum){
+                                /*Nothing to add*/
+                        }
+                        else if(hydrology_model==HydrologypismEnum){
                                 /*Nothing to add*/
                         }

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

@@ -13 +13 @@
 
         /*intermediary*/
-        int                                                     hydrology_model;
-        int                                                     solution_type;
-        int                                                     numoutputs                              =       0;
-        bool                                            save_results;
-        bool                                            modify_loads                    =       true;
-        char                                    **requested_outputs = NULL;
-        IssmDouble                      ThawedNodes;
+        int          hydrology_model;
+        int          solution_type;
+        int          numoutputs        = 0;
+        bool         save_results;
+        bool         modify_loads      = true;
+        char       **requested_outputs = NULL;
+        IssmDouble   ThawedNodes;
 
         /*first recover parameters common to all solutions*/
@@ -46 +46 @@
         else if (hydrology_model==HydrologydcEnum){
                 /*intermediary: */
-                bool                            isefficientlayer;
-                bool                            isthermal;
-                int                                     step,hydroslices;
-                IssmDouble      time,init_time,hydrotime,yts;
-                IssmDouble      dt,hydrodt;
+                bool       isefficientlayer;
+                bool       isthermal;
+                int        step,hydroslices;
+                IssmDouble time,init_time,hydrotime,yts;
+                IssmDouble dt,hydrodt;
 
                 femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

@@ -204 +204 @@
                 case 2: return HydrologyshreveEnum;
                 case 3: return HydrologyshaktiEnum;
+                case 4: return HydrologypismEnum;
                 default: _error_("Marshalled hydrology code \""<<enum_in<<"\" not supported yet"); 
         }

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp

@@ -115 +115 @@
         KSPGetPC(ksp,&pc);
         if (solver_type==MUMPSPACKAGE_LU){
-                #if defined(_HAVE_PETSCDEV_)
+                #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=9)
                 PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
                 #else

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

@@ -47 +47 @@
                         WriteData(fid,prefix,'name','md.hydrology.model','data',4,'format','Integer');
                         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s
+                        WriteData(fid,prefix,'data',{'Watercolumn'},'name','md.hydrology.requested_outputs','format','StringArray');
                 end % }}}
         end

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

@@ -8 +8 @@
                 spcwatercolumn     = NaN;
                 stabilization      = 0;
-    requested_outputs  = {};
+                requested_outputs  = {};
         end
         methods