HydrologyDCInefficientAnalysis Class Reference

#include <HydrologyDCInefficientAnalysis.h>

Inheritance diagram for HydrologyDCInefficientAnalysis:

Public Member Functions
int	DofsPerNode (int **doflist, int domaintype, int approximation)
void	UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
void	UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
void	CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false)
void	CreateConstraints (Constraints *constraints, IoModel *iomodel)
void	CreateLoads (Loads *loads, IoModel *iomodel)
void	Core (FemModel *femmodel)
ElementVector *	CreateDVector (Element *element)
ElementMatrix *	CreateJacobianMatrix (Element *element)
ElementMatrix *	CreateKMatrix (Element *element)
ElementVector *	CreatePVector (Element *element)
void	GetSolutionFromInputs (Vector< IssmDouble > *solution, Element *element)
void	GradientJ (Vector< IssmDouble > *gradient, Element *element, int control_type, int control_index)
void	InputUpdateFromSolution (IssmDouble *solution, Element *element)
void	UpdateConstraints (FemModel *femmodel)
void	GetB (IssmDouble *B, Element *element, IssmDouble *xyz_list, Gauss *gauss)
IssmDouble	SedimentStoring (Element *element, Gauss *gauss, Input2 *sed_head_input, Input2 *base_input)
IssmDouble	SedimentTransmitivity (Element *element, Gauss *gauss, Input2 *sed_head_input, Input2 *base_input, Input2 *SedTrans_input)
void	GetHydrologyDCInefficientHmax (IssmDouble *ph_max, Element *element, Node *innode)
IssmDouble	GetHydrologyKMatrixTransfer (Element *element)
IssmDouble	GetHydrologyPVectorTransfer (Element *element, Gauss *gauss, Input2 *epl_head_input)
void	ElementizeEplMask (FemModel *femmodel)
void	HydrologyIDSGetMask (Vector< IssmDouble > *vec_mask, Element *element)
void	ElementizeIdsMask (FemModel *femmodel)
void	HydrologyIdsGetActive (Vector< IssmDouble > *active_vec, Element *element)
Public Member Functions inherited from Analysis
virtual	~Analysis ()

Detailed Description

Definition at line 12 of file HydrologyDCInefficientAnalysis.h.

Member Function Documentation

DofsPerNode()

int HydrologyDCInefficientAnalysis::DofsPerNode ( int **  doflist, int  domaintype, int  approximation  )

virtual

Implements Analysis.

Definition at line 10 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                              {/*{{{*/
   return 1;
}/*}}}*/

UpdateParameters()

void HydrologyDCInefficientAnalysis::UpdateParameters ( Parameters *  parameters, IoModel *  iomodel, int  solution_enum, int  analysis_enum  )

virtual

Implements Analysis.

Definition at line 13 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                                {/*{{{*/
 
   int         hydrology_model;
   int         sedimentlimit_flag;
   int         transfer_flag;
   int         unconfined_flag;
   int         penalty_lock;
   int         hydro_maxiter;
   int         hydroslices;
   int         averaging_method;
   int         numoutputs;
   bool        isefficientlayer;
   IssmDouble  penalty_factor;
   IssmDouble  rel_tol;
   IssmDouble  leakagefactor;
   IssmDouble  sedimentlimit;
   char**      requestedoutputs = NULL;
 
   /*retrieve some parameters: */
   bool   issmb;
   iomodel->FindConstant(&issmb,"md.transient.issmb");
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want DC?*/
   if(hydrology_model!=HydrologydcEnum) return;
 
   iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
   iomodel->FetchData(&transfer_flag,      "md.hydrology.transfer_flag" );
   iomodel->FetchData(&unconfined_flag,    "md.hydrology.unconfined_flag" );
   iomodel->FetchData(&penalty_lock,       "md.hydrology.penalty_lock" );
   iomodel->FetchData(&hydro_maxiter,      "md.hydrology.max_iter" );
   iomodel->FetchData(&hydroslices,        "md.hydrology.steps_per_step");
   iomodel->FetchData(&averaging_method,   "md.hydrology.averaging");
   iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
   iomodel->FetchData(&penalty_factor,     "md.hydrology.penalty_factor" );
   iomodel->FetchData(&rel_tol,            "md.hydrology.rel_tol" );
 
   parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
   parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
   parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
   parameters->AddObject(new IntParam(HydrologydcUnconfinedFlagEnum,unconfined_flag));
   parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
   parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
   parameters->AddObject(new IntParam(HydrologyStepsPerStepEnum,hydroslices));
   parameters->AddObject(new IntParam(HydrologyAveragingEnum,averaging_method));
 
   parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
   parameters->AddObject(new DoubleParam(HydrologydcPenaltyFactorEnum,penalty_factor));
   parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
   if(transfer_flag==1){
      iomodel->FetchData(&leakagefactor,"md.hydrology.leakage_factor");
      parameters->AddObject(new DoubleParam(HydrologydcLeakageFactorEnum,leakagefactor));
   }
   if(sedimentlimit_flag==1){
      iomodel->FetchData(&sedimentlimit,"md.hydrology.sedimentlimit");
      parameters->AddObject(new DoubleParam(HydrologydcSedimentlimitEnum,sedimentlimit));
   }
   if(!issmb){
      parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
      parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
   }
 
  /*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");
}/*}}}*/

UpdateElements()

void HydrologyDCInefficientAnalysis::UpdateElements ( Elements *  elements, Inputs2 *  inputs2, IoModel *  iomodel, int  analysis_counter, int  analysis_type  )

virtual

Implements Analysis.

Definition at line 81 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                                              {/*{{{*/
 
   bool   isefficientlayer;
   int    hydrology_model;
 
   /*Fetch data needed: */
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want DC?*/
   if(hydrology_model!=HydrologydcEnum) return;
 
   /*Fetch data needed: */
   iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
 
   /*Update elements: */
   int counter=0;
   for(int i=0;i<iomodel->numberofelements;i++){
      if(iomodel->my_elements[i]){
         Element* element=(Element*)elements->GetObjectByOffset(counter);
         element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,P1Enum);
         counter++;
      }
   }
   iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.basal_moulin_input",HydrologydcBasalMoulinInputEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.sediment_head",SedimentHeadSubstepEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum);
   if(iomodel->domaintype!=Domain2DhorizontalEnum){
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
   }
   if(isefficientlayer){
      iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.mask_eplactive_node",HydrologydcMaskEplactiveNodeEnum);
      iomodel->FetchDataToInput(inputs2,elements,"md.initialization.epl_head",EplHeadSubstepEnum);
   }
 
}/*}}}*/

CreateNodes()

void HydrologyDCInefficientAnalysis::CreateNodes ( Nodes *  nodes, IoModel *  iomodel, bool  isamr = false  )

virtual

Implements Analysis.

Definition at line 122 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                        {/*{{{*/
 
   /*Fetch parameters: */
   int  hydrology_model;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want DC?*/
   if(hydrology_model!=HydrologydcEnum) return;
 
   if(iomodel->domaintype!=Domain2DhorizontalEnum){
      iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
   }
   ::CreateNodes(nodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum);
   iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
}/*}}}*/

CreateConstraints()

void HydrologyDCInefficientAnalysis::CreateConstraints ( Constraints *  constraints, IoModel *  iomodel  )

virtual

Implements Analysis.

Definition at line 137 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                               {/*{{{*/
 
   /*retrieve some parameters: */
   int hydrology_model;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
   if(hydrology_model!=HydrologydcEnum) return;
 
   IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spcsediment_head",HydrologyDCInefficientAnalysisEnum,P1Enum);
}/*}}}*/

CreateLoads()

void HydrologyDCInefficientAnalysis::CreateLoads ( Loads *  loads, IoModel *  iomodel  )

virtual

Implements Analysis.

Definition at line 146 of file HydrologyDCInefficientAnalysis.cpp.

                                                                              {/*{{{*/
 
   /*Fetch parameters: */
   int hydrology_model;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
   if(hydrology_model!=HydrologydcEnum) return;
 
   if(iomodel->domaintype==Domain3DEnum){
      iomodel->FetchData(1,"md.mesh.vertexonbase");
   }
   //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 Pengrid(i+1,i,iomodel));
            loads->AddObject(new Moulin(i+1,i,iomodel));
         }
      }
      else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
         if(iomodel->my_vertices[i]){
            loads->AddObject(new Pengrid(i+1,i,iomodel));
            loads->AddObject(new Moulin(i+1,i,iomodel));
         }
      }
   }
   iomodel->DeleteData(1,"md.mesh.vertexonbase");
}/*}}}*/

Core()

void HydrologyDCInefficientAnalysis::Core ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 177 of file HydrologyDCInefficientAnalysis.cpp.

                                                                     {/*{{{*/
   _error_("not implemented");
}/*}}}*/

CreateDVector()

ElementVector * HydrologyDCInefficientAnalysis::CreateDVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 180 of file HydrologyDCInefficientAnalysis.cpp.

                                                                            {/*{{{*/
   /*Default, return NULL*/
   return NULL;
}/*}}}*/

CreateJacobianMatrix()

ElementMatrix * HydrologyDCInefficientAnalysis::CreateJacobianMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 184 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                   {/*{{{*/
_error_("Not implemented");
}/*}}}*/

CreateKMatrix()

ElementMatrix * HydrologyDCInefficientAnalysis::CreateKMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 187 of file HydrologyDCInefficientAnalysis.cpp.

                                                                            {/*{{{*/
 
   /*Intermediaries*/
   bool     thawed_element;
   int      domaintype;
   Element* basalelement;
 
   /*Get basal element*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         basalelement = element;
         break;
      case Domain3DEnum:
         if(!element->IsOnBase()) return NULL;
         basalelement = element->SpawnBasalElement();
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
 
   basalelement->GetInput2Value(&thawed_element,HydrologydcMaskThawedEltEnum);
 
   /*Check that all nodes are active, else return empty matrix*/
   if(!thawed_element) {
      if(domaintype!=Domain2DhorizontalEnum){
         basalelement->DeleteMaterials();
         delete basalelement;
      }
      return NULL;
   }
 
   /*Intermediaries */
   bool        active_element,isefficientlayer;
   IssmDouble  D_scalar,Jdet,dt;
   IssmDouble  sediment_transmitivity;
   IssmDouble  transfer,sediment_storing;
   IssmDouble *xyz_list  = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = basalelement->GetNumberOfNodes();
 
   /*Initialize Element vector*/
   ElementMatrix* Ke     = basalelement->NewElementMatrix();
   IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
   IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
   IssmDouble     D[2][2]= {0.};
 
   /*Retrieve all inputs and parameters*/
   basalelement ->GetVerticesCoordinates(&xyz_list);
   basalelement ->FindParam(&dt,TimesteppingTimeStepEnum);
   basalelement ->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
   Input2* SedTrans_input = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
   Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
   Input2* base_input     = basalelement->GetInput2(BaseEnum);
   Input2* old_wh_input   = basalelement->GetInput2(SedimentHeadOldEnum);                  _assert_(old_wh_input);
 
   /*Transfer related Inputs*/
   if(isefficientlayer){
      basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
   }
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=basalelement->NewGauss(2);
 
   for(int ig=gauss -> begin();ig<gauss->end();ig++){
      gauss          -> GaussPoint(ig);
      basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
 
      sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
      sediment_storing       = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
 
      /*Diffusivity*/
      D_scalar=sediment_transmitivity*gauss->weight*Jdet;
      //D_scalar=gauss->weight*Jdet;
      if(dt!=0.) D_scalar=D_scalar*dt;
      D[0][0]=D_scalar;
      D[1][1]=D_scalar;
      GetB(B,basalelement,xyz_list,gauss);
      TripleMultiply(B,2,numnodes,1,
                            &D[0][0],2,2,0,
                            B,2,numnodes,0,
                            &Ke->values[0],1);
 
      /*Transient*/
      if(dt!=0.){
         basalelement->NodalFunctions(&basis[0],gauss);
         D_scalar=sediment_storing*gauss->weight*Jdet;
         //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
         TripleMultiply(basis,numnodes,1,0,
                               &D_scalar,1,1,0,
                               basis,1,numnodes,0,
                               &Ke->values[0],1);
 
         /*Transfer EPL part*/
         if(isefficientlayer){
            if(active_element){
               transfer=GetHydrologyKMatrixTransfer(basalelement);
               basalelement->NodalFunctions(&basis[0],gauss);
               D_scalar=dt*transfer*gauss->weight*Jdet;
               //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
               TripleMultiply(basis,numnodes,1,0,
                                     &D_scalar,1,1,0,
                                     basis,1,numnodes,0,
                                     &Ke->values[0],1);
            }
         }
      }
   }
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(B);
   xDelete<IssmDouble>(basis);
   delete gauss;
   if(domaintype!=Domain2DhorizontalEnum){
      basalelement->DeleteMaterials();
      delete basalelement;
   }
   return Ke;
}/*}}}*/

CreatePVector()

ElementVector * HydrologyDCInefficientAnalysis::CreatePVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 306 of file HydrologyDCInefficientAnalysis.cpp.

                                                                            {/*{{{*/
 
   /*Intermediaries*/
   bool      thawed_element;
   int          domaintype;
   Element* basalelement;
 
   /*Get basal element*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         basalelement = element;
         break;
      case Domain3DEnum:
         if(!element->IsOnBase()) return NULL;
         basalelement = element->SpawnBasalElement();
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
 
   basalelement->GetInput2Value(&thawed_element,HydrologydcMaskThawedEltEnum);
 
   /*Check that all nodes are active, else return empty matrix*/
   if(!thawed_element) {
   if(domaintype!=Domain2DhorizontalEnum){
         basalelement->DeleteMaterials();
         delete basalelement;
      }
      return NULL;
   }
 
   /*Intermediaries */
   bool       active_element,isefficientlayer;
   int        smb_model;
   int        smbsubstepping;
   int        hydrologysubstepping;
   int        smb_averaging;
   IssmDouble dt,scalar,sediment_storing;
   IssmDouble water_head,sediment_transmitivity;
   IssmDouble water_load,runoff_value,transfer;
   IssmDouble Jdet;
 
   IssmDouble *xyz_list             = NULL;
   Input2     *active_element_input = NULL;
   Input2     *old_wh_input         = NULL;
   Input2     *dummy_input          = NULL;
   Input2     *surface_runoff_input = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = basalelement->GetNumberOfNodes();
 
   /*Initialize Element vector*/
   ElementVector* pe    = basalelement->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   basalelement->GetVerticesCoordinates(&xyz_list);
   basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
   basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
   basalelement->FindParam(&smb_model,SmbEnum);
   basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
 
   Input2*  sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
   Input2*  epl_head_input   = basalelement->GetInput2(EplHeadSubstepEnum);
   Input2*  base_input       = basalelement->GetInput2(BaseEnum);
   Input2*  basal_melt_input = basalelement->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input);
   Input2*  SedTrans_input   = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
 
   IssmDouble time;
   basalelement->FindParam(&time,TimeEnum);
 
   if(dt!= 0.){
      old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
   }
   if(smb_model==SMBgradientscomponentsEnum){
      basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
      basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
 
      if(smbsubstepping==1){
         //no substeping for the smb we take the result from there
         dummy_input = basalelement->GetInput2(SmbRunoffEnum); _assert_(dummy_input);
      }
      else if(smbsubstepping>1 && smbsubstepping<=hydrologysubstepping){
         //finer hydro stepping, we take the value at the needed time
         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
      }
      else{
         //finer stepping in smb, we average the runoff from transient input
         dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
      }
      surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
   }
 
   /*Transfer related Inputs*/
   if(isefficientlayer){
      basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
   }
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=basalelement->NewGauss(2);
 
   IssmDouble yts;
   basalelement->FindParam(&yts,ConstantsYtsEnum);
 
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
      basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
      basalelement->NodalFunctions(basis,gauss);
      sediment_transmitivity = SedimentTransmitivity(basalelement,gauss,sed_head_input,base_input,SedTrans_input);
 
      /*Loading term*/
      if(!isefficientlayer){
         basal_melt_input->GetInputValue(&water_load,gauss);
         if(surface_runoff_input) surface_runoff_input->GetInputValue(&runoff_value,gauss);
         else                     runoff_value = 0.;
         scalar = Jdet*gauss->weight*(water_load+runoff_value);
         //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
         if(dt!=0.) scalar = scalar*dt;
         for(int i=0;i<numnodes;i++){
            pe->values[i]+=scalar*basis[i];
         }
      }
      else{
         /*if EPL is present and active input is there not here*/
         if(!active_element){
            basal_melt_input->GetInputValue(&water_load,gauss);
            if(surface_runoff_input)surface_runoff_input->GetInputValue(&runoff_value,gauss);
            else runoff_value = 0.;
            scalar = Jdet*gauss->weight*(water_load+runoff_value);
            //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
            if(dt!=0.) scalar = scalar*dt;
            for(int i=0;i<numnodes;i++){
               pe->values[i]+=scalar*basis[i];
            }
         }
      }
 
      /*Transient and transfer terms*/
      if(dt!=0.){
         old_wh_input->GetInputValue(&water_head,gauss);
         sediment_storing = SedimentStoring(basalelement,gauss,sed_head_input,base_input);
         if(isefficientlayer){
            /*Dealing with the sediment part of the transfer term*/
            if(active_element){
               transfer=GetHydrologyPVectorTransfer(basalelement,gauss,epl_head_input);
            }
            else{
               transfer=0.0;
            }
            scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer));
            //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity;
            for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
         }
         else{
            scalar = Jdet*gauss->weight*(water_head*sediment_storing);
            //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity;
            for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
         }
      }
   }
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   delete gauss;
   if(domaintype!=Domain2DhorizontalEnum){
      basalelement->DeleteMaterials();
      delete basalelement;
   }
   return pe;
}/*}}}*/

GetSolutionFromInputs()

void HydrologyDCInefficientAnalysis::GetSolutionFromInputs ( Vector< IssmDouble > *  solution, Element *  element  )

virtual

Implements Analysis.

Definition at line 503 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                       {/*{{{*/
   element->GetSolutionFromInputsOneDof(solution,SedimentHeadSubstepEnum);
}/*}}}*/

GradientJ()

void HydrologyDCInefficientAnalysis::GradientJ ( Vector< IssmDouble > *  gradient, Element *  element, int  control_type, int  control_index  )

virtual

Implements Analysis.

Definition at line 506 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                              {/*{{{*/
   _error_("Not implemented yet");
}/*}}}*/

InputUpdateFromSolution()

void HydrologyDCInefficientAnalysis::InputUpdateFromSolution ( IssmDouble *  solution, Element *  element  )

virtual

Implements Analysis.

Definition at line 509 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                 {/*{{{*/
 
   /*Intermediaries*/
   int    domaintype;
   Element* basalelement=NULL;
   bool   converged;
   int*     doflist = NULL;
 
   /*Get basal element*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         basalelement = element;
         break;
      case Domain3DEnum:
         if(!element->IsOnBase()) return;
         basalelement = element->SpawnBasalElement();
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
   /*Fetch number of nodes for this finite element*/
   int numnodes = basalelement->GetNumberOfNodes();
 
   /*Fetch dof list and allocate solution vector*/
   basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
   IssmDouble* values   = xNew<IssmDouble>(numnodes);
   IssmDouble* pressure = xNew<IssmDouble>(numnodes);
   IssmDouble* residual = xNew<IssmDouble>(numnodes);
 
   for(int i=0;i<numnodes;i++){
      values[i] =solution[doflist[i]];
      if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
 
   }
 
   /*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;
      IssmDouble* thickness = xNew<IssmDouble>(numnodes);
      IssmDouble* base = xNew<IssmDouble>(numnodes);
      IssmDouble* transmitivity = xNew<IssmDouble>(numnodes);
 
      basalelement->FindParam(&kmax,HydrologySedimentKmaxEnum);
      basalelement->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
      IssmDouble rho_freshwater = basalelement->FindParam(MaterialsRhoFreshwaterEnum);
      IssmDouble rho_ice        = basalelement->FindParam(MaterialsRhoIceEnum);
      IssmDouble g              = basalelement->FindParam(ConstantsGEnum);
 
      basalelement->GetInputListOnVertices(&thickness[0],ThicknessEnum);
      basalelement->GetInputListOnVertices(&base[0],BaseEnum);
      basalelement->GetInputListOnVertices(&transmitivity[0],HydrologydcSedimentTransmitivityEnum);
 
      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);
            //residual[i] = kappa*(values[i]-h_max)*transmitivity[i];
         }
         else{
            residual[i] = 0.;
         }
         pressure[i]=(rho_ice*g*thickness[i])-(rho_freshwater*g*(values[i]-base[i]));
      }
      xDelete<IssmDouble>(thickness);
      xDelete<IssmDouble>(base);
      xDelete<IssmDouble>(transmitivity);
   }
 
   /*Add input to the element: */
   element->AddBasalInput2(SedimentHeadSubstepEnum,values,P1Enum);
   element->AddBasalInput2(EffectivePressureSubstepEnum,pressure,P1Enum);
   element->AddBasalInput2(SedimentHeadResidualEnum,residual,P1Enum);
 
   /*Free ressources:*/
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(residual);
   xDelete<IssmDouble>(pressure);
   xDelete<int>(doflist);
   if(domaintype!=Domain2DhorizontalEnum){
      basalelement->DeleteMaterials();
      delete basalelement;
   }
}/*}}}*/

UpdateConstraints()

void HydrologyDCInefficientAnalysis::UpdateConstraints ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 598 of file HydrologyDCInefficientAnalysis.cpp.

                                                                        {/*{{{*/
   /*Default, do nothing*/
   return;
}/*}}}*/

GetB()

void HydrologyDCInefficientAnalysis::GetB	(	IssmDouble *	B,
		Element *	element,
		IssmDouble *	xyz_list,
		Gauss *	gauss
	)

Definition at line 476 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                         {/*{{{*/
   /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2.
    * For node i, Bi can be expressed in the actual coordinate system
    * by:
    *       Bi=[ dN/dx ]
    *          [ dN/dy ]
    * where N is the finiteelement function for node i.
    *
    * We assume B has been allocated already, of size: 3x(2*numnodes)
    */
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Get nodal functions derivatives*/
   IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
   element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
   /*Build B: */
   for(int i=0;i<numnodes;i++){
      B[numnodes*0+i] = dbasis[0*numnodes+i];
      B[numnodes*1+i] = dbasis[1*numnodes+i];
   }
 
   /*Clean-up*/
   xDelete<IssmDouble>(dbasis);
}/*}}}*/

SedimentStoring()

IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring	(	Element *	element,
		Gauss *	gauss,
		Input2 *	sed_head_input,
		Input2 *	base_input
	)

Definition at line 603 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                                  {/*{{{*/
   int unconf_scheme;
   IssmDouble expfac;
   IssmDouble sediment_storing;
   IssmDouble storing,yield;
   IssmDouble base_elev,prestep_head,water_sheet;
   IssmDouble rho_freshwater           = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble g                        = element->FindParam(ConstantsGEnum);
   IssmDouble sediment_porosity        = element->FindParam(HydrologydcSedimentPorosityEnum);
   IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);
   IssmDouble sediment_compressibility = element->FindParam(HydrologydcSedimentCompressibilityEnum);
   IssmDouble water_compressibility    = element->FindParam(HydrologydcWaterCompressibilityEnum);
   element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
   switch(unconf_scheme){
   case 0:
      sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
      break;
   case 1:
      base_input->GetInputValue(&base_elev,gauss);
      sed_head_input->GetInputValue(&prestep_head,gauss);
      water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
 
      /* if (water_sheet<sediment_thickness){ */
      /*    sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); */
      /* } */
      /* else{ */
      /*    sediment_storing=sediment_porosity; */
      /* } */
      storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
      //using logistic function for heavyside approximation
      expfac=10.;
      yield=sediment_porosity;
      sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
      break;
   default:
      _error_("UnconfinedFlag is 0 or 1");
   }
   return sediment_storing;
}/*}}}*/

SedimentTransmitivity()

IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity	(	Element *	element,
		Gauss *	gauss,
		Input2 *	sed_head_input,
		Input2 *	base_input,
		Input2 *	SedTrans_input
	)

Definition at line 642 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                                                               {/*{{{*/
   int unconf_scheme;
   IssmDouble ratio,expfac;
   IssmDouble sediment_transmitivity;
   IssmDouble FullLayer_transmitivity;
   IssmDouble meltingrate;
   IssmDouble groundedice;
   IssmDouble base_elev,prestep_head,water_sheet;
   IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);
 
   element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
   SedTrans_input->GetInputValue(&FullLayer_transmitivity,gauss);
 
   switch(unconf_scheme){
   case 0:
      sediment_transmitivity=FullLayer_transmitivity;
      break;
   case 1:
      base_input->GetInputValue(&base_elev,gauss);
      sed_head_input->GetInputValue(&prestep_head,gauss);
      water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
 
      //min definition of the if test
      sediment_transmitivity=FullLayer_transmitivity*min(water_sheet/sediment_thickness,1.);
      if (sediment_transmitivity==0){
         sediment_transmitivity=1.0e-20;
      }
 
      break;
   default:
      _error_("UnconfinedFlag is 0 or 1");
   }
   return sediment_transmitivity;
}/*}}}*/

GetHydrologyDCInefficientHmax()

void HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax	(	IssmDouble *	ph_max,
		Element *	element,
		Node *	innode
	)

Definition at line 676 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                    {/*{{{*/
 
   int        hmax_flag;
   IssmDouble h_max;
   IssmDouble rho_ice,rho_water;
   IssmDouble thickness,bed;
   /*Get the flag to the limitation method*/
   element->FindParam(&hmax_flag,HydrologydcSedimentlimitFlagEnum);
 
   /*Switch between the different cases*/
   switch(hmax_flag){
   case 0:
      h_max=1.0e+10;
      break;
   case 1:
      element->FindParam(&h_max,HydrologydcSedimentlimitEnum);
      break;
   case 2:
      /*Compute max*/
      rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
      rho_ice   = element->FindParam(MaterialsRhoIceEnum);
      element->GetInputValue(&thickness,innode,ThicknessEnum);
      element->GetInputValue(&bed,innode,BaseEnum);
      h_max=((rho_ice*thickness)/rho_water)+bed;
      break;
   case 3:
      _error_("Using normal stress  not supported yet");
      break;
   default:
      _error_("no case higher than 3 for SedimentlimitFlag");
   }
   /*Assign output pointer*/
   *ph_max=h_max;
}

GetHydrologyKMatrixTransfer()

IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyKMatrixTransfer

(

Element *

element

)

Definition at line 711 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                      {/*{{{*/
 
   int transfermethod;
   IssmDouble leakage,transfer;
   element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
 
   /*Switch between the different transfer methods cases*/
   switch(transfermethod){
   case 0:
      /*Just keepping the transfer to zero*/
      transfer=0.0;
      break;
   case 1:
      element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
      transfer=+leakage;
      break;
   default:
      _error_("no case higher than 1 for the Transfer method");
   }
   return transfer;
}/*}}}*/

GetHydrologyPVectorTransfer()

IssmDouble HydrologyDCInefficientAnalysis::GetHydrologyPVectorTransfer	(	Element *	element,
		Gauss *	gauss,
		Input2 *	epl_head_input
	)

Definition at line 732 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                                            {/*{{{*/
 
   int transfermethod;
   IssmDouble epl_head;
   IssmDouble leakage,transfer;
   element->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
 
   /*Switch between the different transfer methods cases*/
   switch(transfermethod){
   case 0:
      /*Just keepping the transfer to zero*/
      transfer=0.0;
      break;
   case 1:
      _assert_(epl_head_input);
      epl_head_input->GetInputValue(&epl_head,gauss);
      element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
      transfer=+epl_head*leakage;
      break;
   default:
      _error_("no case higher than 1 for the Transfer method");
   }
   return transfer;
}/*}}}*/

ElementizeEplMask()

void HydrologyDCInefficientAnalysis::ElementizeEplMask

(

FemModel *

femmodel

)

Definition at line 756 of file HydrologyDCInefficientAnalysis.cpp.

                                                                        {/*{{{*/
 
   bool     element_active;
   Element* element=NULL;
 
   for(int i=0;i<femmodel->elements->Size();i++){
      element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 
      Input2* input=element->GetInput2(HydrologydcMaskEplactiveNodeEnum); _assert_(input);
      if(input->GetInputMax()>0.){
         element_active = true;
      }
      else{
         element_active = false;
      }
      element->SetBoolInput(element->inputs2,HydrologydcMaskEplactiveEltEnum,element_active);
   }
}/*}}}*/

HydrologyIDSGetMask()

void HydrologyDCInefficientAnalysis::HydrologyIDSGetMask	(	Vector< IssmDouble > *	vec_mask,
		Element *	element
	)

Definition at line 774 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                       {/*{{{*/
   bool        active_element;
   int         domaintype;
   Element*    basalelement=NULL;
 
   /*Get basal element*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         basalelement = element;
         break;
      case Domain3DEnum:
         if(!element->IsOnBase()) return;
         basalelement = element->SpawnBasalElement();
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
   /*Intermediaries*/
   int            numnodes    =  basalelement->GetNumberOfNodes();
   IssmDouble*    meltingrate =  xNew<IssmDouble>(numnodes);
   IssmDouble*    groundedice =  xNew<IssmDouble>(numnodes);
 
   basalelement->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum);
   basalelement->GetInputListOnVertices(&groundedice[0],MaskOceanLevelsetEnum);
 
   /*if melting rate is not positive and node is not floating, deactivate*/
   for(int i=0;i<numnodes;i++){
      if ((meltingrate[i]<=0.0) && (groundedice[i]>0)){
         vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
      }
      else{
         vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
      }
   }
 
   if(domaintype!=Domain2DhorizontalEnum){
      basalelement->DeleteMaterials();
      delete basalelement;
   }
   xDelete<IssmDouble>(meltingrate);
   xDelete<IssmDouble>(groundedice);
}/*}}}*/

ElementizeIdsMask()

void HydrologyDCInefficientAnalysis::ElementizeIdsMask

(

FemModel *

femmodel

)

Definition at line 816 of file HydrologyDCInefficientAnalysis.cpp.

                                                                        {/*{{{*/
 
   bool     element_active;
   Element* element=NULL;
 
   for(int i=0;i<femmodel->elements->Size();i++){
      element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
 
      Input2* input=element->GetInput2(HydrologydcMaskThawedNodeEnum); _assert_(input);
      if(input->GetInputMax()>0.){
         element_active = true;
      }
      else{
         element_active = false;
      }
      element->SetBoolInput(element->inputs2,HydrologydcMaskThawedEltEnum,element_active);
   }
}/*}}}*/

HydrologyIdsGetActive()

void HydrologyDCInefficientAnalysis::HydrologyIdsGetActive	(	Vector< IssmDouble > *	active_vec,
		Element *	element
	)

Definition at line 834 of file HydrologyDCInefficientAnalysis.cpp.

                                                                                                          {/*{{{*/
   /*Constants*/
   int      domaintype;
   Element*   basalelement=NULL;
 
   /*Get basal element*/
   element->FindParam(&domaintype,DomainTypeEnum);
   switch(domaintype){
      case Domain2DhorizontalEnum:
         basalelement = element;
         break;
      case Domain3DEnum:
         if(!element->IsOnBase()) return;
         basalelement = element->SpawnBasalElement();
         break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
   }
 
   const int   numnodes = basalelement->GetNumberOfNodes();
   IssmDouble  flag     = 0.;
   IssmDouble* active   = xNew<IssmDouble>(numnodes);
   bool active_element;
 
   /*Pass the activity mask from elements to nodes*/
   basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskThawedNodeEnum);
   basalelement->GetInput2Value(&active_element,HydrologydcMaskThawedEltEnum);
 
   for(int i=0;i<numnodes;i++) flag+=active[i];
 
   /*If any node is active all the node in the element are active*/
   if(flag>0.){
      for(int i=0;i<numnodes;i++){
         active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
      }
   }
   /*If the element is active all its nodes are active*/
   else if(active_element){
      for(int i=0;i<numnodes;i++){
         active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
      }
   }
   else{
      /*Do not do anything: at least one node is active for this element but this element is not solved for*/
   }
   if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
   xDelete<IssmDouble>(active);
}

---

The documentation for this class was generated from the following files:

• src/c/analyses/HydrologyDCInefficientAnalysis.h
• src/c/analyses/HydrologyDCInefficientAnalysis.cpp