HydrologyGlaDSAnalysis Class Reference

#include <HydrologyGlaDSAnalysis.h>

Inheritance diagram for HydrologyGlaDSAnalysis:

Public Member Functions
void	CreateConstraints (Constraints *constraints, IoModel *iomodel)
void	CreateLoads (Loads *loads, IoModel *iomodel)
void	CreateNodes (Nodes *nodes, IoModel *iomodel, bool isamr=false)
int	DofsPerNode (int **doflist, int domaintype, int approximation)
void	UpdateElements (Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)
void	UpdateParameters (Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)
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	SetChannelCrossSectionOld (FemModel *femmodel)
void	UpdateSheetThickness (FemModel *femmodel)
void	UpdateSheetThickness (Element *element)
void	UpdateChannelCrossSection (FemModel *femmodel)
void	UpdateEffectivePressure (FemModel *femmodel)
void	UpdateEffectivePressure (Element *element)
Public Member Functions inherited from Analysis
virtual	~Analysis ()

Detailed Description

Definition at line 11 of file HydrologyGlaDSAnalysis.h.

Member Function Documentation

CreateConstraints()

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

virtual

Implements Analysis.

Definition at line 10 of file HydrologyGlaDSAnalysis.cpp.

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

CreateLoads()

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

virtual

Implements Analysis.

Definition at line 21 of file HydrologyGlaDSAnalysis.cpp.

                                                                      {/*{{{*/
 
   /*Fetch parameters: */
   int  hydrology_model;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want GlaDS?*/
   if(hydrology_model!=HydrologyGlaDSEnum) return;
 
   /*Add channels?*/
   int K,L;
   bool ischannels;
   IssmDouble* channelarea;
   iomodel->FindConstant(&ischannels,"md.hydrology.ischannels");
   iomodel->FetchData(&channelarea,&K,&L,"md.initialization.channelarea");
   if(ischannels){
      /*Get faces (edges in 2d)*/
      CreateFaces(iomodel);
      for(int i=0;i<iomodel->numberoffaces;i++){
 
         /*Get left and right elements*/
         int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
 
         /*Now, if this element is not in the partition*/
         if(!iomodel->my_elements[element]) continue;
 
         /*Add channelarea from initialization if exists*/
         if(K!=0 && K!=iomodel->numberoffaces){
            _error_("Unknown dimension for channel area initialization.");
         }
         if(K==0){
            loads->AddObject(new Channel(i+1,0.,i,iomodel));
         }
         else{
            loads->AddObject(new Channel(i+1,channelarea[i],i,iomodel));
         }
         iomodel->DeleteData(1,"md.initialization.channelarea");
      }
   }
 
   /*Create discrete loads for Moulins*/
   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(i+1,i,iomodel));
         }
      }
      else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
         if(iomodel->my_vertices[i]){
            loads->AddObject(new Moulin(i+1,i,iomodel));
         }  
      }
   }
   iomodel->DeleteData(1,"md.mesh.vertexonbase");
 
   /*Deal with Neumann BC*/
   int M,N;
   int *segments = NULL;
   iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
 
   /*Check that the size seem right*/
   _assert_(N==3); _assert_(M>=3);
   for(int i=0;i<M;i++){
      if(iomodel->my_elements[segments[i*3+2]-1]){
         loads->AddObject(new Neumannflux(i+1,i,iomodel,segments));
      }
   }
   xDelete<int>(segments);
 
}/*}}}*/

CreateNodes()

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

virtual

Implements Analysis.

Definition at line 93 of file HydrologyGlaDSAnalysis.cpp.

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

DofsPerNode()

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

virtual

Implements Analysis.

Definition at line 106 of file HydrologyGlaDSAnalysis.cpp.

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

UpdateElements()

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

virtual

Implements Analysis.

Definition at line 109 of file HydrologyGlaDSAnalysis.cpp.

                                                                                                                                      {/*{{{*/
 
   /*Fetch data needed: */
   int    hydrology_model,frictionlaw;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want GlaDS?*/
   if(hydrology_model!=HydrologyGlaDSEnum) return;
 
   /*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.geometry.bed",BedEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
   if(iomodel->domaintype!=Domain2DhorizontalEnum){
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
      iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
   }
   iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
   iomodel->FindConstant(&frictionlaw,"md.friction.law");
 
   /*Friction law variables*/
   switch(frictionlaw){
      case 1:
         iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
         iomodel->FetchDataToInput(inputs2,elements,"md.friction.p",FrictionPEnum);
         iomodel->FetchDataToInput(inputs2,elements,"md.friction.q",FrictionQEnum);
         break;
      case 8:
         iomodel->FetchDataToInput(inputs2,elements,"md.friction.coefficient",FrictionCoefficientEnum);
         break;
      default:
         _error_("Friction law "<< frictionlaw <<" not supported");
   }
}/*}}}*/

UpdateParameters()

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

virtual

Implements Analysis.

Definition at line 163 of file HydrologyGlaDSAnalysis.cpp.

                                                                                                                        {/*{{{*/
 
   /*retrieve some parameters: */
   int    hydrology_model;
   int    numoutputs;
   char** requestedoutputs = NULL;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want GlaDS?*/
   if(hydrology_model!=HydrologyGlaDSEnum) return;
 
   parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
   parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.pressure_melt_coefficient",HydrologyPressureMeltCoefficientEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.cavity_spacing",HydrologyCavitySpacingEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ischannels",HydrologyIschannelsEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.melt_flag",HydrologyMeltFlagEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
 
   /*Deal with friction parameters*/
   int frictionlaw;
   iomodel->FindConstant(&frictionlaw,"md.friction.law");
   if(frictionlaw==6){
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
   }
   if(frictionlaw==4){
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
   }
   if(frictionlaw==1 || frictionlaw==3 || frictionlaw==7){
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
   }
   if(frictionlaw==9){
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
      parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
      parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
   }
 
   /*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");
}/*}}}*/

Core()

void HydrologyGlaDSAnalysis::Core ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 213 of file HydrologyGlaDSAnalysis.cpp.

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

CreateDVector()

ElementVector * HydrologyGlaDSAnalysis::CreateDVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 216 of file HydrologyGlaDSAnalysis.cpp.

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

CreateJacobianMatrix()

ElementMatrix * HydrologyGlaDSAnalysis::CreateJacobianMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 220 of file HydrologyGlaDSAnalysis.cpp.

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

CreateKMatrix()

ElementMatrix * HydrologyGlaDSAnalysis::CreateKMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 223 of file HydrologyGlaDSAnalysis.cpp.

                                                                    {/*{{{*/
 
   /*Intermediaries */
   IssmDouble  Jdet,dphi[3],h,k;
   IssmDouble  A,B,n,phi_old,phi,phi_0,H,b,v1;
   IssmDouble* xyz_list = NULL;
 
   /*Hard coded coefficients*/
   const IssmPDouble alpha = 5./4.;
   const IssmPDouble beta  = 3./2.;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementMatrix* Ke     = element->NewElementMatrix();
   IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
   IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
 
   /*Get all inputs and parameters*/
   IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
   Input2* k_input   = element->GetInput2(HydrologySheetConductivityEnum);_assert_(k_input);
   Input2* phi_input = element->GetInput2(HydraulicPotentialEnum);      _assert_(phi_input);
   Input2* h_input   = element->GetInput2(HydrologySheetThicknessEnum); _assert_(h_input);
   Input2* H_input      = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input      = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* B_input      = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input      = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
 
      phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
      phi_input->GetInputValue(&phi,gauss);
      h_input->GetInputValue(&h,gauss);
      k_input->GetInputValue(&k,gauss);
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      b_input->GetInputValue(&b,gauss);
      H_input->GetInputValue(&H,gauss);
 
      /*Get norm of gradient of hydraulic potential and make sure it is >0*/
      IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
      if(normgradphi < AEPS) normgradphi = AEPS;
 
      IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);
 
      /*Diffusive term*/
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
                     coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
                     + coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
         }
      }
 
      /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
      phi_0 = rho_water*g*b + rho_ice*g*H;
      A=pow(B,-n);
      v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
         }
      }
 
      /*Transient term if dt>0*/
      if(dt>0.){
         /*Diffusive term*/
         for(int i=0;i<numnodes;i++){
            for(int j=0;j<numnodes;j++){
               Ke->values[i*numnodes+j] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*basis[i]*basis[j];
            }
         }
      }
 
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   xDelete<IssmDouble>(dbasis);
   delete gauss;
   return Ke;
}/*}}}*/

CreatePVector()

ElementVector * HydrologyGlaDSAnalysis::CreatePVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 321 of file HydrologyGlaDSAnalysis.cpp.

                                                                    {/*{{{*/
 
   /*Skip if water or ice shelf element*/
   if(element->IsFloating()) return NULL;
 
   /*Intermediaries */
   bool        meltflag;
   IssmDouble  Jdet,w,v2,vx,vy,ub,h,h_r;
   IssmDouble  G,m,frictionheat,alpha2;
   IssmDouble  A,B,n,phi_old,phi,phi_0;
   IssmDouble  H,b;
   IssmDouble* xyz_list = NULL;
 
   /*Fetch number of nodes and dof for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Initialize Element vector and other vectors*/
   ElementVector* pe    = element->NewElementVector();
   IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->FindParam(&meltflag,HydrologyMeltFlagEnum);
   IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble l_r       = element->FindParam(HydrologyCavitySpacingEnum);
   IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
   Input2* hr_input     = element->GetInput2(HydrologyBumpHeightEnum);_assert_(hr_input);
   Input2* vx_input     = element->GetInput2(VxEnum);_assert_(vx_input);
   Input2* vy_input     = element->GetInput2(VyEnum);_assert_(vy_input);
   Input2* h_input      = element->GetInput2(HydrologySheetThicknessEnum);_assert_(h_input);
   Input2* H_input      = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input      = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* G_input      = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
   Input2* m_input      = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);_assert_(m_input);
   Input2* B_input      = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input      = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* phiold_input = element->GetInput2(HydraulicPotentialOldEnum);      _assert_(phiold_input);
   Input2* phi_input    = element->GetInput2(HydraulicPotentialEnum);         _assert_(phi_input);
 
   /*Build friction element, needed later: */
   Friction* friction=new Friction(element,2);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss(2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
 
      element->JacobianDeterminant(&Jdet,xyz_list,gauss);
      element->NodalFunctions(basis,gauss);
 
      /*Get input values at gauss points*/
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      h_input->GetInputValue(&h,gauss);
      G_input->GetInputValue(&G,gauss);
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      hr_input->GetInputValue(&h_r,gauss);
      phi_input->GetInputValue(&phi,gauss);
      b_input->GetInputValue(&b,gauss);
      H_input->GetInputValue(&H,gauss);
 
      /*Get basal velocity*/
      ub = sqrt(vx*vx + vy*vy);
 
      /*Compute cavity opening w*/
      w  = 0.;
      if(h<h_r) w = ub*(h_r-h)/l_r;
 
      /*Compute frictional heat flux*/
      friction->GetAlpha2(&alpha2,gauss);
      frictionheat=alpha2*ub*ub;
 
      /*Compute melt (if necessary)*/
      if(!meltflag){
         m = (G + frictionheat)/(rho_ice*L);
      }
      else{
         m_input->GetInputValue(&m,gauss);
      }
 
      /*Compute closing rate*/
      phi_0 = rho_water*g*b + rho_ice*g*H;
      A=pow(B,-n);
      v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
 
      for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
 
      /*Transient term if dt>0*/
      if(dt>0.){
         phiold_input->GetInputValue(&phi_old,gauss);
         for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
      }
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   delete gauss;
   delete friction;
   return pe;
}/*}}}*/

GetSolutionFromInputs()

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

virtual

Implements Analysis.

Definition at line 427 of file HydrologyGlaDSAnalysis.cpp.

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

GradientJ()

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

virtual

Implements Analysis.

Definition at line 432 of file HydrologyGlaDSAnalysis.cpp.

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

InputUpdateFromSolution()

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

virtual

Implements Analysis.

Definition at line 435 of file HydrologyGlaDSAnalysis.cpp.

                                                                                                   {/*{{{*/
   element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
}/*}}}*/

UpdateConstraints()

void HydrologyGlaDSAnalysis::UpdateConstraints ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 438 of file HydrologyGlaDSAnalysis.cpp.

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

SetChannelCrossSectionOld()

void HydrologyGlaDSAnalysis::SetChannelCrossSectionOld

(

FemModel *

femmodel

)

Definition at line 444 of file HydrologyGlaDSAnalysis.cpp.

                                                                        {/*{{{*/
 
   bool ischannels;
   femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
   if(!ischannels) return;
 
   for(int i=0;i<femmodel->loads->Size();i++){
      if(femmodel->loads->GetEnum(i)==ChannelEnum){
         Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
         channel->SetChannelCrossSectionOld();
      }
   }
 
}/*}}}*/

UpdateSheetThickness() [1/2]

void HydrologyGlaDSAnalysis::UpdateSheetThickness

(

FemModel *

femmodel

)

Definition at line 458 of file HydrologyGlaDSAnalysis.cpp.

                                                                   {/*{{{*/
 
   for(int j=0;j<femmodel->elements->Size();j++){
      Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
      UpdateSheetThickness(element);
   }
 
}/*}}}*/

UpdateSheetThickness() [2/2]

void HydrologyGlaDSAnalysis::UpdateSheetThickness

(

Element *

element

)

Definition at line 466 of file HydrologyGlaDSAnalysis.cpp.

                                                                 {/*{{{*/
 
   /*Skip if water or ice shelf element*/
   if(element->IsFloating()) return;
 
   /*Intermediaries */
   IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r,H,b;
   IssmDouble  A,B,n,phi,phi_0;
   IssmDouble  alpha,beta;
 
   /*Fetch number vertices for this element*/
   int numvertices = element->GetNumberOfVertices();
 
   /*Initialize new sheet thickness*/
   IssmDouble* h_new = xNew<IssmDouble>(numvertices);
 
   /*Retrieve all inputs and parameters*/
   IssmDouble  dt       = element->FindParam(TimesteppingTimeStepEnum);
   IssmDouble  l_r      = element->FindParam(HydrologyCavitySpacingEnum);
   IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble g         = element->FindParam(ConstantsGEnum);
   Input2* hr_input = element->GetInput2(HydrologyBumpHeightEnum);_assert_(hr_input);
   Input2* vx_input = element->GetInput2(VxEnum);_assert_(vx_input);
   Input2* vy_input = element->GetInput2(VyEnum);_assert_(vy_input);
   Input2* H_input  = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input  = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* hold_input  = element->GetInput2(HydrologySheetThicknessOldEnum);_assert_(hold_input);
   Input2* B_input  = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input  = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* phi_input = element->GetInput2(HydraulicPotentialEnum);         _assert_(phi_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss();
   for(int iv=0;iv<numvertices;iv++){
      gauss->GaussVertex(iv);
 
      /*Get input values at gauss points*/
      phi_input->GetInputValue(&phi,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      hold_input->GetInputValue(&h_old,gauss);
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      hr_input->GetInputValue(&h_r,gauss);
      b_input->GetInputValue(&b,gauss);
      H_input->GetInputValue(&H,gauss);
 
      /*Get values for a few potentials*/
      phi_0 = rho_water*g*b + rho_ice*g*H;
      N = phi_0 - phi;
 
      /*Get basal velocity*/
      ub = sqrt(vx*vx + vy*vy);
 
      /*Get A from B and n*/
      A=pow(B,-n);
 
      /*Define alpha and beta*/
      if(h_old<h_r){
         alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
         beta  = ub*h_r/l_r;
      }
      else{
         alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
         beta  = 0.;
      }
 
      /*Get new sheet thickness*/
      h_new[iv] = ODE1(alpha,beta,h_old,dt,1);
 
      /*Make sure it is positive*/
      if(h_new[iv]<AEPS) h_new[iv] = AEPS;
   }
 
   element->AddInput2(HydrologySheetThicknessEnum,h_new,P1Enum);
 
   /*Clean up and return*/
   xDelete<IssmDouble>(h_new);
   delete gauss;
}/*}}}*/

UpdateChannelCrossSection()

void HydrologyGlaDSAnalysis::UpdateChannelCrossSection

(

FemModel *

femmodel

)

Definition at line 607 of file HydrologyGlaDSAnalysis.cpp.

                                                                        {/*{{{*/
 
   bool ischannels;
   femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum);
   if(!ischannels) return;
 
   for(int i=0;i<femmodel->loads->Size();i++){
      if(femmodel->loads->GetEnum(i)==ChannelEnum){
         Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i);
         channel->UpdateChannelCrossSection();
      }
   }
 
}/*}}}*/

UpdateEffectivePressure() [1/2]

void HydrologyGlaDSAnalysis::UpdateEffectivePressure

(

FemModel *

femmodel

)

Definition at line 547 of file HydrologyGlaDSAnalysis.cpp.

                                                                      {/*{{{*/
 
   for(int j=0;j<femmodel->elements->Size();j++){
      Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
      UpdateEffectivePressure(element);
   }
 
}/*}}}*/

UpdateEffectivePressure() [2/2]

void HydrologyGlaDSAnalysis::UpdateEffectivePressure

(

Element *

element

)

Definition at line 555 of file HydrologyGlaDSAnalysis.cpp.

                                                                    {/*{{{*/
 
   /*Skip if water or ice shelf element*/
   if(element->IsFloating()) return;
   
   /*Intermediary*/
   IssmDouble phi_0, phi_m, p_i;
   IssmDouble H,b,phi;
 
   int numnodes = element->GetNumberOfNodes();
 
   /*Get thickness and base on nodes to apply cap on water head*/
   IssmDouble* N = xNew<IssmDouble>(numnodes);
   IssmDouble  rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble  rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   IssmDouble  g         = element->FindParam(ConstantsGEnum);
   Input2* H_input   = element->GetInput2(ThicknessEnum); _assert_(H_input);
   Input2* b_input   = element->GetInput2(BedEnum); _assert_(b_input);
   Input2* phi_input = element->GetInput2(HydraulicPotentialEnum); _assert_(phi_input);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss();
   for(int iv=0;iv<numnodes;iv++){
      gauss->GaussNode(element->FiniteElement(),iv);
 
      /*Get input values at gauss points*/
      H_input->GetInputValue(&H,gauss);
      b_input->GetInputValue(&b,gauss);
      phi_input->GetInputValue(&phi,gauss);
 
      /*Elevation potential*/
      phi_m = rho_water*g*b;
 
      /*Compute overburden pressure*/
      p_i = rho_ice*g*H;
 
      /*Copmute overburden potential*/
      phi_0 = phi_m + p_i;
 
      /*Calculate effective pressure*/
      N[iv] = phi_0 - phi;
 
      if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(N[iv])) _error_("Inf found in solution vector");
   }
 
   element->AddInput2(EffectivePressureEnum,N,element->FiniteElement());
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(N);
}/*}}}*/

---

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

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