source: issm/trunk/src/c/analyses/HydrologySommersAnalysis.cpp@ 22758

#include "./HydrologySommersAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

/*Define 2 hardcoded parameters*/
#define OMEGA 0.001    // parameter controlling transition to nonlinear resistance in basal system (dimensionless)
#define NU    1.787e-6 //kinematic water viscosity m^2/s
#define CT    7.5e-8  // Clapeyron slope (K/Pa) 
#define CW    4.22e3   // specific heat capacity of water (J/kg/K)

/*Model processing*/
void HydrologySommersAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/

        /*retrieve some parameters: */
        int hydrology_model;
        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");

        if(hydrology_model!=HydrologysommersEnum) return;

        IoModelToConstraintsx(constraints,iomodel,"md.hydrology.spchead",HydrologySommersAnalysisEnum,P1Enum);

}/*}}}*/
void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/

        /*Fetch parameters: */
        int  hydrology_model;
        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");

        /*Now, do we really want Sommers?*/
        if(hydrology_model!=HydrologysommersEnum) 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(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
                        }
                }
                else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
                        if(iomodel->my_vertices[i]){
                                loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
                        }       
                }
        }
        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(iomodel->loadcounter+i+1,i,iomodel,segments,HydrologySommersAnalysisEnum));
                }
        }

        xDelete<int>(segments);

}/*}}}*/
void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/

        /*Fetch parameters: */
        int  hydrology_model;
        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");

        /*Now, do we really want Sommers?*/
        if(hydrology_model!=HydrologysommersEnum) return;

        if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
        ::CreateNodes(nodes,iomodel,HydrologySommersAnalysisEnum,P1Enum);
        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
}/*}}}*/
int  HydrologySommersAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
        return 1;
}/*}}}*/
void HydrologySommersAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

        /*Fetch data needed: */
        int    hydrology_model,frictionlaw;
        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");

        /*Now, do we really want Sommers?*/
        if(hydrology_model!=HydrologysommersEnum) 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(i,iomodel,analysis_counter,analysis_type,P1Enum);
                        counter++;
                }
        }

        iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
        iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
        if(iomodel->domaintype!=Domain2DhorizontalEnum){
                iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
                iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
        }
        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
        iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
        iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
        iomodel->FetchDataToInput(elements,"md.basalforcings.geothermalflux",BasalforcingsGeothermalfluxEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.head",HydrologyHeadEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.gap_height",HydrologyGapHeightEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.englacial_input",HydrologyEnglacialInputEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.bump_spacing",HydrologyBumpSpacingEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.reynolds",HydrologyReynoldsEnum);
        iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
        iomodel->FindConstant(&frictionlaw,"md.friction.law");

        /*Friction law variables*/
        switch(frictionlaw){
                case 1:
                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
                        break;
                case 8:
                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
                        break;
                default:
                        _error_("Friction law "<< frictionlaw <<" not supported");
        }
}/*}}}*/
void HydrologySommersAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

        /*retrieve some parameters: */
        int    hydrology_model;
        int    numoutputs;
        char** requestedoutputs = NULL;
        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");

        /*Now, do we really want Sommers?*/
        if(hydrology_model!=HydrologysommersEnum) 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");
}/*}}}*/

/*Finite Element Analysis*/
void           HydrologySommersAnalysis::Core(FemModel* femmodel){/*{{{*/
        _error_("not implemented");
}/*}}}*/
ElementVector* HydrologySommersAnalysis::CreateDVector(Element* element){/*{{{*/
        /*Default, return NULL*/
        return NULL;
}/*}}}*/
ElementMatrix* HydrologySommersAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
_error_("Not implemented");
}/*}}}*/
ElementMatrix* HydrologySommersAnalysis::CreateKMatrix(Element* element){/*{{{*/

        /*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;
}/*}}}*/
ElementVector* HydrologySommersAnalysis::CreatePVector(Element* element){/*{{{*/

        /*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->GetMaterialParameter(MaterialsLatentheatEnum);
        IssmDouble  g               = element->GetMaterialParameter(ConstantsGEnum);
        IssmDouble  rho_ice         = element->GetMaterialParameter(MaterialsRhoIceEnum);
        IssmDouble  rho_water       = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
        Input* head_input           = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
        Input* gap_input            = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
        Input* thickness_input      = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
        Input* base_input           = element->GetInput(BaseEnum);                       _assert_(base_input);
        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
        Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
        Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
        Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
   Input* headold_input        = element->GetInput(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;
}/*}}}*/
void           HydrologySommersAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
        element->GetSolutionFromInputsOneDof(solution,HydrologyHeadEnum);
}/*}}}*/
void           HydrologySommersAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
        _error_("Not implemented yet");
}/*}}}*/
void           HydrologySommersAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

        /*Intermediary*/
        IssmDouble dh[3];
        int* doflist = NULL;
        IssmDouble* xyz_list = NULL;

        /*Get gravity from parameters*/
           IssmDouble  g = element->GetMaterialParameter(ConstantsGEnum);

        /*Fetch number of nodes for this finite element*/
        int numnodes = element->GetNumberOfNodes();

        /*Fetch dof list and allocate solution vector*/
        element->GetDofList(&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->GetMaterialParameter(MaterialsRhoIceEnum);
        IssmDouble  rho_water = element->GetMaterialParameter(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->AddInput(HydrologyHeadEnum,values,element->GetElementType());
   element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum);

        /*Update reynolds number according to new solution*/
        element->GetVerticesCoordinates(&xyz_list);
        Input* head_input = element->GetInput(HydrologyHeadEnum);_assert_(head_input);
        head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list);
        IssmDouble conductivity = GetConductivity(element);

        IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/NU;
        element->AddInput(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);
}/*}}}*/
void           HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
        /*Default, do nothing*/
        return;
}/*}}}*/

/*Additional methods*/
IssmDouble HydrologySommersAnalysis::GetConductivity(Element* element){/*{{{*/

        /*Intermediaries */
        IssmDouble gap,reynolds;

        /*Get gravity from parameters*/
        IssmDouble  g = element->GetMaterialParameter(ConstantsGEnum);

        /*Get Reynolds and gap average values*/
        Input* reynolds_input = element->GetInput(HydrologyReynoldsEnum);  _assert_(reynolds_input);
        Input* gap_input      = element->GetInput(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;
}/*}}}*/
void HydrologySommersAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/


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

}/*}}}*/
void HydrologySommersAnalysis::UpdateGapHeight(Element* element){/*{{{*/

        /*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->GetMaterialParameter(MaterialsLatentheatEnum);
        IssmDouble  g               = element->GetMaterialParameter(ConstantsGEnum);
        IssmDouble  rho_ice         = element->GetMaterialParameter(MaterialsRhoIceEnum);
        IssmDouble  rho_water       = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(geothermalflux_input);
        Input* head_input           = element->GetInput(HydrologyHeadEnum);              _assert_(head_input);
        Input* gap_input            = element->GetInput(HydrologyGapHeightEnum);         _assert_(gap_input);
        Input* thickness_input      = element->GetInput(ThicknessEnum);                  _assert_(thickness_input);
        Input* base_input           = element->GetInput(BaseEnum);                       _assert_(base_input);
        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
        Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
        Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
        Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
        Input* br_input             = element->GetInput(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->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
 
        /*Divide by connectivity, add basal flux as an input*/
        q = q/totalweights;
        element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);

        /* Divide by connectivity, add degree of channelization as an input */
        channelization = channelization/totalweights;
        element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum);

        /*Clean up and return*/
        xDelete<IssmDouble>(xyz_list);
        delete friction;
        delete gauss;
}/*}}}*/