HydrologyShaktiAnalysis Class Reference

#include <HydrologyShaktiAnalysis.h>

Inheritance diagram for HydrologyShaktiAnalysis:

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)
IssmDouble	GetConductivity (Element *element)
void	UpdateGapHeight (FemModel *femmodel)
void	UpdateGapHeight (Element *element)
Public Member Functions inherited from Analysis
virtual	~Analysis ()

Detailed Description

Definition at line 11 of file HydrologyShaktiAnalysis.h.

Member Function Documentation

CreateConstraints()

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

virtual

Implements Analysis.

Definition at line 14 of file HydrologyShaktiAnalysis.cpp.

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

CreateLoads()

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

virtual

Implements Analysis.

Definition at line 25 of file HydrologyShaktiAnalysis.cpp.

                                                                       {/*{{{*/
 
   /*Fetch parameters: */
   int  hydrology_model;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want Shakti?*/
   if(hydrology_model!=HydrologyshaktiEnum) return;
 
   /*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 HydrologyShaktiAnalysis::CreateNodes ( Nodes *  nodes, IoModel *  iomodel, bool  isamr = false  )

virtual

Implements Analysis.

Definition at line 66 of file HydrologyShaktiAnalysis.cpp.

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

DofsPerNode()

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

virtual

Implements Analysis.

Definition at line 79 of file HydrologyShaktiAnalysis.cpp.

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

UpdateElements()

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

virtual

Implements Analysis.

Definition at line 82 of file HydrologyShaktiAnalysis.cpp.

                                                                                                                                       {/*{{{*/
 
   /*Fetch data needed: */
   int    hydrology_model,frictionlaw;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want Shakti?*/
   if(hydrology_model!=HydrologyshaktiEnum) 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);
   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.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.head",HydrologyHeadEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.gap_height",HydrologyGapHeightEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.reynolds",HydrologyReynoldsEnum);
   iomodel->FetchDataToInput(inputs2,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
   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 HydrologyShaktiAnalysis::UpdateParameters ( Parameters *  parameters, IoModel *  iomodel, int  solution_enum, int  analysis_enum  )

virtual

Implements Analysis.

Definition at line 137 of file HydrologyShaktiAnalysis.cpp.

                                                                                                                         {/*{{{*/
 
   /*retrieve some parameters: */
   int    hydrology_model;
   int    numoutputs;
   char** requestedoutputs = NULL;
   iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
 
   /*Now, do we really want Shakti?*/
   if(hydrology_model!=HydrologyshaktiEnum) return;
 
   parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
   parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
   parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
 
  /*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 HydrologyShaktiAnalysis::Core ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 161 of file HydrologyShaktiAnalysis.cpp.

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

CreateDVector()

ElementVector * HydrologyShaktiAnalysis::CreateDVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 164 of file HydrologyShaktiAnalysis.cpp.

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

CreateJacobianMatrix()

ElementMatrix * HydrologyShaktiAnalysis::CreateJacobianMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 168 of file HydrologyShaktiAnalysis.cpp.

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

CreateKMatrix()

ElementMatrix * HydrologyShaktiAnalysis::CreateKMatrix ( Element *  element )

virtual

Implements Analysis.

Definition at line 171 of file HydrologyShaktiAnalysis.cpp.

                                                                     {/*{{{*/
 
   /*Intermediaries */
   IssmDouble Jdet;
   IssmDouble* xyz_list = NULL;
 
   /*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 conductivity from inputs*/
   IssmDouble conductivity = GetConductivity(element);
 
   /*Get englacial storage coefficient*/
   IssmDouble storage,dt;
   element->FindParam(&storage,HydrologyStorageEnum);
   element->FindParam(&dt,TimesteppingTimeStepEnum);
 
   /* Start  looping on the number of gaussian points: */
   Gauss* gauss=element->NewGauss(1);
   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);
 
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
              + gauss->weight*Jdet*storage/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 * HydrologyShaktiAnalysis::CreatePVector ( Element *  element )

virtual

Implements Analysis.

Definition at line 220 of file HydrologyShaktiAnalysis.cpp.

                                                                     {/*{{{*/
 
   /*Skip if water or ice shelf element*/
   if(element->IsFloating()) return NULL;
 
   /*Intermediaries */
   IssmDouble  Jdet,meltrate,G,dh[2],B,A,n;
   IssmDouble  gap,bed,thickness,head,ieb,head_old;
   IssmDouble  lr,br,vx,vy,beta;
   IssmDouble  alpha2,frictionheat;
   IssmDouble  PMPheat,dpressure_water[2],dbed[2]; 
   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);
   IssmDouble  latentheat      = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble  g               = element->FindParam(ConstantsGEnum);
   IssmDouble  rho_ice         = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble  rho_water       = element->FindParam(MaterialsRhoFreshwaterEnum);
   Input2* geothermalflux_input = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
   Input2* head_input           = element->GetInput2(HydrologyHeadEnum);              _assert_(head_input);
   Input2* gap_input            = element->GetInput2(HydrologyGapHeightEnum);         _assert_(gap_input);
   Input2* thickness_input      = element->GetInput2(ThicknessEnum);                  _assert_(thickness_input);
   Input2* base_input           = element->GetInput2(BaseEnum);                       _assert_(base_input);
   Input2* B_input              = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input              = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* englacial_input      = element->GetInput2(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
   Input2* vx_input             = element->GetInput2(VxEnum);                         _assert_(vx_input);
   Input2* vy_input             = element->GetInput2(VyEnum);                         _assert_(vy_input);
   Input2* lr_input             = element->GetInput2(HydrologyBumpSpacingEnum);       _assert_(lr_input);
   Input2* br_input             = element->GetInput2(HydrologyBumpHeightEnum);        _assert_(br_input);
   Input2* headold_input        = element->GetInput2(HydrologyHeadOldEnum);           _assert_(headold_input);
 
   /*Get conductivity from inputs*/
   IssmDouble conductivity = GetConductivity(element);
 
   /*Get englacial storage coefficient*/
   IssmDouble storage,dt;
   element->FindParam(&storage,HydrologyStorageEnum);
   element->FindParam(&dt,TimesteppingTimeStepEnum);
 
   /*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);
      geothermalflux_input->GetInputValue(&G,gauss);
      base_input->GetInputValue(&bed,gauss);
      base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
      thickness_input->GetInputValue(&thickness,gauss);
      gap_input->GetInputValue(&gap,gauss);
      head_input->GetInputValue(&head,gauss);
      head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
      englacial_input->GetInputValue(&ieb,gauss);
      lr_input->GetInputValue(&lr,gauss);
      br_input->GetInputValue(&br,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      headold_input->GetInputValue(&head_old,gauss);
 
      /*Get ice A parameter*/
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      A=pow(B,-n);
 
      /*Compute beta term*/
      if(gap<br)
       beta = (br-gap)/lr;
      else
       beta = 0.;
 
      /*Compute frictional heat flux*/
      friction->GetAlpha2(&alpha2,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      frictionheat=alpha2*(vx*vx+vy*vy);
 
      /*Get water and ice pressures*/
      IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
      IssmDouble pressure_water = rho_water*g*(head-bed);
      if(pressure_water>pressure_ice) pressure_water = pressure_ice;
 
      /*Get water pressure from previous time step to use in lagged creep term*/
      IssmDouble pressure_water_old = rho_water*g*(head_old-bed);
      if(pressure_water_old>pressure_ice) pressure_water_old = pressure_ice;
 
      /*Compute change in sensible heat due to changes in pressure melting point*/
      dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
      dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
      PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
 
      meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
      _assert_(meltrate>0.);
 
      for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
       (
        meltrate*(1/rho_water-1/rho_ice)
        +A*pow(fabs(pressure_ice - pressure_water),n-1)*(pressure_ice - pressure_water)*gap
        -beta*sqrt(vx*vx+vy*vy)
        +ieb
        +storage*head_old/dt
        )*basis[i];        
   }
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(basis);
   delete friction;
   delete gauss;
   return pe;
}/*}}}*/

GetSolutionFromInputs()

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

virtual

Implements Analysis.

Definition at line 342 of file HydrologyShaktiAnalysis.cpp.

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

GradientJ()

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

virtual

Implements Analysis.

Definition at line 345 of file HydrologyShaktiAnalysis.cpp.

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

InputUpdateFromSolution()

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

virtual

Implements Analysis.

Definition at line 348 of file HydrologyShaktiAnalysis.cpp.

                                                                                                    {/*{{{*/
 
   /*Intermediary*/
   IssmDouble dh[3];
   int* doflist = NULL;
   IssmDouble* xyz_list = NULL;
 
   /*Get gravity from parameters*/
   IssmDouble  g = element->FindParam(ConstantsGEnum);
 
   /*Fetch number of nodes for this finite element*/
   int numnodes = element->GetNumberOfNodes();
 
   /*Fetch dof list and allocate solution vector*/
   element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
   IssmDouble* values = xNew<IssmDouble>(numnodes);
 
   /*Get thickness and base on nodes to apply cap on water head*/
   IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
   IssmDouble* thickness = xNew<IssmDouble>(numnodes);
   IssmDouble* bed       = xNew<IssmDouble>(numnodes);
   IssmDouble  rho_ice   = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble  rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
   element->GetInputListOnNodes(&thickness[0],ThicknessEnum);
   element->GetInputListOnNodes(&bed[0],BaseEnum);
 
   /*Get head from previous time-step and under-relaxation coefficient to use in under-relaxation for nonlinear convergence*/
   IssmDouble* head_old  = xNew<IssmDouble>(numnodes); 
   element->GetInputListOnNodes(&head_old[0],HydrologyHeadEnum);
   IssmDouble relaxation; 
   element->FindParam(&relaxation,HydrologyRelaxationEnum);
 
   /*Use the dof list to index into the solution vector: */
   for(int i=0;i<numnodes;i++){
      values[i]=solution[doflist[i]];
 
      /*make sure that p_water<p_ice ->  h<rho_i H/rho_w + zb*/
      if(values[i]>rho_ice*thickness[i]/rho_water+bed[i]){
         values[i] = rho_ice*thickness[i]/rho_water+bed[i];
      }
 
      /*Make sure that negative pressure is not allowed*/
  //    if(values[i]<bed[i]){
   //    values[i] = bed[i];
   // }
 
      /*Under-relaxation*/
      values[i] = head_old[i] - relaxation*(head_old[i]-values[i]);
 
      /*Calculate effective pressure*/
      eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
 
      if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
      if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
   }
 
   /*Add input to the element: */
   element->AddInput2(HydrologyHeadEnum,values,element->GetElementType());
   element->AddInput2(EffectivePressureEnum,eff_pressure,P1Enum);
 
   /*Update reynolds number according to new solution*/
   element->GetVerticesCoordinates(&xyz_list);
   Input2* head_input = element->GetInput2(HydrologyHeadEnum);_assert_(head_input);
   IssmDouble conductivity = GetConductivity(element);
 
   /*Get gap height derivatives at the center of the element*/
   Gauss* gauss=element->NewGauss(1);
   head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
   delete gauss;
 
   IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
   element->AddInput2(HydrologyReynoldsEnum,&reynolds,P0Enum);
 
   /*Free resources:*/
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(thickness);
   xDelete<IssmDouble>(bed);
   xDelete<IssmDouble>(xyz_list);
   xDelete<int>(doflist);
   xDelete<IssmDouble>(eff_pressure);
   xDelete<IssmDouble>(head_old);
}/*}}}*/

UpdateConstraints()

void HydrologyShaktiAnalysis::UpdateConstraints ( FemModel *  femmodel )

virtual

Implements Analysis.

Definition at line 430 of file HydrologyShaktiAnalysis.cpp.

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

GetConductivity()

IssmDouble HydrologyShaktiAnalysis::GetConductivity

(

Element *

element

)

Definition at line 436 of file HydrologyShaktiAnalysis.cpp.

                                                                   {/*{{{*/
 
   /*Intermediaries */
   IssmDouble gap,reynolds;
 
   /*Get gravity from parameters*/
   IssmDouble  g = element->FindParam(ConstantsGEnum);
 
   /*Get Reynolds and gap average values*/
   Input2* reynolds_input = element->GetInput2(HydrologyReynoldsEnum);  _assert_(reynolds_input);
   Input2* gap_input      = element->GetInput2(HydrologyGapHeightEnum); _assert_(gap_input);
   reynolds_input->GetInputAverage(&reynolds);
   gap_input->GetInputAverage(&gap);
 
   /*Compute conductivity*/
   IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
   _assert_(conductivity>0);
 
   /*Clean up and return*/
   return conductivity;
}/*}}}*/

UpdateGapHeight() [1/2]

void HydrologyShaktiAnalysis::UpdateGapHeight

(

FemModel *

femmodel

)

Definition at line 457 of file HydrologyShaktiAnalysis.cpp.

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

UpdateGapHeight() [2/2]

void HydrologyShaktiAnalysis::UpdateGapHeight

(

Element *

element

)

Definition at line 465 of file HydrologyShaktiAnalysis.cpp.

                                                             {/*{{{*/
 
   /*Skip if water or ice shelf element*/
   if(element->IsFloating()) return;
 
   /*Intermediaries */
   IssmDouble newgap = 0.;
   IssmDouble  Jdet,meltrate,G,dh[2],B,A,n,dt;
   IssmDouble  gap,bed,thickness,head,ieb;
   IssmDouble  lr,br,vx,vy,beta;
   IssmDouble  alpha2,frictionheat;
   IssmDouble* xyz_list = NULL;
   IssmDouble  dpressure_water[2],dbed[2],PMPheat;
   IssmDouble q = 0.;
   IssmDouble channelization = 0.;
 
   /*Retrieve all inputs and parameters*/
   element->GetVerticesCoordinates(&xyz_list);
   element->FindParam(&dt,TimesteppingTimeStepEnum);
   IssmDouble  latentheat      = element->FindParam(MaterialsLatentheatEnum);
   IssmDouble  g               = element->FindParam(ConstantsGEnum);
   IssmDouble  rho_ice         = element->FindParam(MaterialsRhoIceEnum);
   IssmDouble  rho_water       = element->FindParam(MaterialsRhoFreshwaterEnum);
   Input2* geothermalflux_input = element->GetInput2(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
   Input2* head_input           = element->GetInput2(HydrologyHeadEnum);              _assert_(head_input);
   Input2* gap_input            = element->GetInput2(HydrologyGapHeightEnum);         _assert_(gap_input);
   Input2* thickness_input      = element->GetInput2(ThicknessEnum);                  _assert_(thickness_input);
   Input2* base_input           = element->GetInput2(BaseEnum);                       _assert_(base_input);
   Input2* B_input              = element->GetInput2(MaterialsRheologyBEnum);         _assert_(B_input);
   Input2* n_input              = element->GetInput2(MaterialsRheologyNEnum);         _assert_(n_input);
   Input2* englacial_input      = element->GetInput2(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
   Input2* vx_input             = element->GetInput2(VxEnum);                         _assert_(vx_input);
   Input2* vy_input             = element->GetInput2(VyEnum);                         _assert_(vy_input);
   Input2* lr_input             = element->GetInput2(HydrologyBumpSpacingEnum);       _assert_(lr_input);
   Input2* br_input             = element->GetInput2(HydrologyBumpHeightEnum);        _assert_(br_input);
 
   /*Get conductivity from inputs*/
   IssmDouble conductivity = GetConductivity(element);
 
   /*Build friction element, needed later: */
   Friction* friction=new Friction(element,2);
 
   /*Keep track of weights*/
   IssmDouble totalweights=0.;
 
   /* 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);
 
      geothermalflux_input->GetInputValue(&G,gauss);
      base_input->GetInputValue(&bed,gauss);
      base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
      thickness_input->GetInputValue(&thickness,gauss);
      gap_input->GetInputValue(&gap,gauss);
      head_input->GetInputValue(&head,gauss);
      head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss);
      englacial_input->GetInputValue(&ieb,gauss);
      lr_input->GetInputValue(&lr,gauss);
      br_input->GetInputValue(&br,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
 
      /*Get ice A parameter*/
      B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      A=pow(B,-n);
 
      /*Compute beta term*/
      if(gap<br)
       beta = (br-gap)/lr;
      else
       beta = 0.;
 
      /*Compute frictional heat flux*/
      friction->GetAlpha2(&alpha2,gauss);
      vx_input->GetInputValue(&vx,gauss);
      vy_input->GetInputValue(&vy,gauss);
      frictionheat=alpha2*(vx*vx+vy*vy);
 
      /*Get water and ice pressures*/
      IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.); 
      IssmDouble pressure_water = rho_water*g*(head-bed);
      if(pressure_water>pressure_ice) pressure_water = pressure_ice;
 
      /* Compute change in sensible heat due to changes in pressure melting point*/
      dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
      dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
      PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
 
      meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
      _assert_(meltrate>0.);
 
      newgap += gauss->weight*Jdet*(gap+dt*(
               meltrate/rho_ice
               -A*pow(fabs(pressure_ice-pressure_water),n-1)*(pressure_ice-pressure_water)*gap
               +beta*sqrt(vx*vx+vy*vy)
               ));
      totalweights +=gauss->weight*Jdet;
 
      /* Compute basal water flux */
      q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
 
      /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */
      channelization += gauss->weight*Jdet*(meltrate/rho_ice/(meltrate/rho_ice+beta*sqrt(vx*vx+vy*vy)));
   }
 
   /*Divide by connectivity*/
   newgap = newgap/totalweights;
   IssmDouble mingap = 1e-6;
   if(newgap<mingap) newgap=mingap;
 
   /*Limit gap height to grow to surface*/
   if(newgap>thickness)
    newgap = thickness;
 
   /*Add new gap as an input*/
   element->AddInput2(HydrologyGapHeightEnum,&newgap,P0Enum);
 
   /*Divide by connectivity, add basal flux as an input*/
   q = q/totalweights;
   element->AddInput2(HydrologyBasalFluxEnum,&q,P0Enum);
 
   /* Divide by connectivity, add degree of channelization as an input */
   channelization = channelization/totalweights;
   element->AddInput2(DegreeOfChannelizationEnum,&channelization,P0Enum);
 
   /*Clean up and return*/
   xDelete<IssmDouble>(xyz_list);
   delete friction;
   delete gauss;
}/*}}}*/

---

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

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