hydrology_core.cpp File Reference

: core of the hydrology solution More...

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"

Go to the source code of this file.

Functions
void	hydrology_core (FemModel *femmodel)

Detailed Description

: core of the hydrology solution

Definition in file hydrology_core.cpp.

Function Documentation

hydrology_core()

void hydrology_core

(

FemModel *

femmodel

)

Definition at line 12 of file hydrology_core.cpp.

                                       {
 
   /*Start profiler*/
   femmodel->profiler->Start(HYDROLOGYCORE);
 
   /*intermediary*/
   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*/
   femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
   femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
   femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum);
   if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum);
 
   /*Using the Shreve based Model*/
   if (hydrology_model==HydrologyshreveEnum){
      if(VerboseSolution()) _printf0_("   computing water heads\n");
      /*first compute slopes: */
      surfaceslope_core(femmodel);
      bedslope_core(femmodel);
      /*and then go to water column*/
      if(VerboseSolution()) _printf0_("   computing water column\n");
      femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
      solutionsequence_nonlinear(femmodel,modify_loads);
      /*transfer water column thickness to old water column thickness: */
      InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
 
   }
 
   /*Using the double continuum model*/
   else if (hydrology_model==HydrologydcEnum){
      /*intermediary: */
      bool       isefficientlayer;
      /*recover parameters: */
      femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 
      /*first we exclude frozen nodes of the solved nodes*/
      femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
      femmodel->HydrologyIDSupdateDomainx(&ThawedNodes);
 
      if(ThawedNodes>0){
         /*check if we need sub steps*/
         int        dtslices;
         femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum);
 
         if(dtslices>1){
            int        substep, numaveragedinput, hydro_averaging;
            IssmDouble global_time, subtime, yts;
            IssmDouble dt, subdt;
 
            femmodel->parameters->FindParam(&global_time,TimeEnum);
            femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
            femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
            femmodel->parameters->FindParam(&hydro_averaging,HydrologyAveragingEnum);
 
            subtime=global_time-dt; //getting the time back to the start of the timestep
            subdt=dt/dtslices; //computing hydro dt from dt and a divider
            substep=0;
            femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum);
            femmodel->parameters->SetParam(subtime,TimeEnum);
 
            /*intermiedaries to deal with averaging*/
            static const int substeplist[4] = {EffectivePressureSubstepEnum,SedimentHeadSubstepEnum,EplHeadSubstepEnum,HydrologydcEplThicknessSubstepEnum};
            static const int transientlist[4] = {EffectivePressureTransientEnum,SedimentHeadTransientEnum,EplHeadTransientEnum,HydrologydcEplThicknessTransientEnum};
            static const int averagelist[4] = {EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
            std::vector<int> substepinput;
            std::vector<int> transientinput;
            std::vector<int> averagedinput;
 
            if (isefficientlayer){
               /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
               numaveragedinput = 4;
               substepinput.assign(substeplist,substeplist+4);
               transientinput.assign(transientlist,transientlist+4);
               averagedinput.assign(averagelist,averagelist+4);
            }
            else{
               numaveragedinput = 2;
               substepinput.assign(substeplist,substeplist+2);
               transientinput.assign(transientlist,transientlist+2);
               averagedinput.assign(averagelist,averagelist+2);
            }
            femmodel->InitTransientInputx(&transientinput[0],numaveragedinput);
            while(substep<dtslices){ //loop on hydro dts
               substep+=1;
               subtime+=subdt;
               /*Setting substep time as global time*/
               femmodel->parameters->SetParam(subtime,TimeEnum);
               if(VerboseSolution()) _printf0_("sub iteration " << substep << "/" << dtslices << "  time [yr]: " << setprecision(4) << subtime/yts << " (time step: " << subdt/yts << ")\n");
               if(VerboseSolution()) _printf0_("   computing water heads\n");
               /*save preceding timestep*/
               InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadOldEnum);
               if (isefficientlayer){
                  InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadOldEnum);
                  InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessOldEnum);
               }
               /*Proceed now to heads computations*/
               solutionsequence_hydro_nonlinear(femmodel);
               /*If we have a sub-timestep we store the substep inputs in a transient input here*/
               femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
            }
            /*averaging the stack*/
            femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
            /*And reseting to global time*/
            femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
            femmodel->parameters->SetParam(global_time,TimeEnum);
         }
         else{
            InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadOldEnum);
            if (isefficientlayer){
               InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadOldEnum);
               InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessOldEnum);
            }
            /*Proceed now to heads computations*/
            if(VerboseSolution()) _printf0_("   computing water heads\n");
            solutionsequence_hydro_nonlinear(femmodel);
            /*If no substeps are present we want to duplicate the results for coupling purposes*/
            InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadEnum);
            InputDuplicatex(femmodel,EffectivePressureSubstepEnum,EffectivePressureEnum);
            if (isefficientlayer){
               InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadEnum);
               InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessEnum);
            }
         }
      }
      if(VerboseSolution())_printf0_("   hydroDC done\n");
   }
 
   /*Using the SHAKTI model*/
   else if (hydrology_model==HydrologyshaktiEnum){
      femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum);
      InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);
      solutionsequence_shakti_nonlinear(femmodel);
      if(VerboseSolution()) _printf0_("   updating gap height\n");
      HydrologyShaktiAnalysis* analysis = new HydrologyShaktiAnalysis();
      analysis->UpdateGapHeight(femmodel);
      delete analysis;
   }
 
   /*Using the GlaDS model*/
   else if (hydrology_model==HydrologyGlaDSEnum){
      HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
      femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
 
      /*Set fields as old*/
      InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
      InputDuplicatex(femmodel,HydrologySheetThicknessEnum,HydrologySheetThicknessOldEnum);
      analysis->SetChannelCrossSectionOld(femmodel);
 
      /*Solve for new potential*/
      solutionsequence_glads_nonlinear(femmodel);
 
      if(VerboseSolution()) _printf0_("   updating effective pressure\n");
      analysis->UpdateEffectivePressure(femmodel);
      delete analysis;
   }
 
   /*Using the PISM hydrology model*/
   else if (hydrology_model==HydrologypismEnum){
      femmodel->SetCurrentConfiguration(HydrologyPismAnalysisEnum);
      if(VerboseSolution()) _printf0_("   updating water column\n");
      HydrologyPismAnalysis* analysis = new HydrologyPismAnalysis();
      InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
      analysis->UpdateWaterColumn(femmodel);
      delete analysis;
   }
   else{
      _error_("Hydrology model "<< EnumToStringx(hydrology_model) <<" not supported yet");
   }
   if(save_results){
      if(VerboseSolution()) _printf0_("   saving hydrology results \n");
      if(hydrology_model==HydrologydcEnum && ThawedNodes==0){
         if(VerboseSolution()) _printf0_("   No thawed node hydro is skiped \n");}
      else{
         femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
      }
   }
   /*Free ressources:*/
   if(numoutputs){
      for(int i=0;i<numoutputs;i++){
         xDelete<char>(requested_outputs[i]);
      }
      xDelete<char*>(requested_outputs);
   }
   /*End profiler*/
   femmodel->profiler->Stop(HYDROLOGYCORE);
}