Changeset 17014

Timestamp:

12/09/13 07:48:15 (11 years ago)

Author:

jbondzio

Message:

CHG: filled EnthalpyAnalysis::ComputeBasalMeltingrate with content. ADD: GetUpperElement, ThermalToEnthalpy for parent class Element*.

Location:

issm/trunk-jpl/src/c

Files:

• analyses/EnthalpyAnalysis.cpp (modified) (12 diffs)
• analyses/EnthalpyAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (2 diffs)
• classes/Elements/Penta.cpp (modified) (14 diffs)
• classes/Elements/Penta.h (modified) (2 diffs)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -403 +403 @@
 
         /*Intermediaries*/
-        int         stabilization;
+        int         i, stabilization;
         IssmDouble  Jdet,phi,dt;
-        IssmDouble  enthalpy;
-        IssmDouble  kappa,tau_parameter,diameter;
+        IssmDouble  enthalpy, Hpmp;
+        IssmDouble  enthalpypicard, d1enthalpypicard[3];
+        IssmDouble  pressure, d1pressure[3], d2pressure;
+        IssmDouble  waterfractionpicard;
+        IssmDouble  kappa,tau_parameter,diameter,kappa_w;
         IssmDouble  u,v,w;
-        IssmDouble  scalar_def,scalar_transient;
+        IssmDouble  scalar_def, scalar_sens ,scalar_transient;
         IssmDouble* xyz_list = NULL;
+        IssmDouble  d1H_d1P, d1P2;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -423 +427 @@
         element->GetVerticesCoordinates(&xyz_list);
         IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
         IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+        IssmDouble  temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
+        IssmDouble  beta                = element->GetMaterialParameter(MaterialsBetaEnum);
+        IssmDouble  latentheat          = element->GetMaterialParameter(MaterialsLatentheatEnum);
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         element->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -429 +437 @@
         Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
         Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
-        Input* enthalpy_input = NULL;
+        Input* enthalpypicard_input=element->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
+        Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+        Input* enthalpy_input=NULL;
         if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
         if(stabilization==2){
@@ -443 +453 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
+                
+                /*viscous dissipation*/
                 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
 
@@ -448 +460 @@
                 if(dt!=0.) scalar_def=scalar_def*dt;
 
-                /*TODO: add -beta*laplace T_m(p)*/
-                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
+                for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
+
+                /*sensible heat flux in temperate ice*/
+                enthalpypicard_input->GetInputValue(&enthalpypicard,gauss);
+                pressure_input->GetInputValue(&pressure,gauss);
+                Hpmp=this->PureIceEnthalpy(element, pressure);
+
+                if(enthalpypicard>=Hpmp){
+                        enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
+                        pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
+                        d2pressure=0.; // for linear elements, 2nd derivative is zero
+                        
+                        d1H_d1P=0.;
+                        for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
+                        d1P2=0.;
+                        for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
+
+                        scalar_sens=-beta*((temperateiceconductivity - thermalconductivity)/latentheat*(d1H_d1P + beta*heatcapacity*d1P2))/rho_ice;
+                        if(dt!=0.) scalar_sens=scalar_sens*dt;
+                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_sens*basis[i];
+                }               
 
                 /* Build transient now */
@@ -455 +486 @@
                         enthalpy_input->GetInputValue(&enthalpy, gauss);
                         scalar_transient=enthalpy*Jdet*gauss->weight;
-                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
+                        for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
                 }
 
@@ -466 +497 @@
                         tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
 
-                        for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+                        for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 
                         if(dt!=0.){
-                                for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+                                for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
                         }
                 }
@@ -484 +515 @@
 ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
 
-        /* Geothermal flux on ice sheet base and basal friction */
+        /* implementation of the basal condition decision chart of Aschwanden 2012, Fig.5 */
         if(!element->IsOnBed() || element->IsFloating()) return NULL;
 
@@ -527 +558 @@
                 watercolumn_input->GetInputValue(&watercolumn,gauss);
 
-                if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
+                if((watercolumn<=0.) && (enthalpy<PureIceEnthalpy(element,pressure))){
                         /* the above check is equivalent to 
-                         NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
+                         NOT [(watercolumn>0.) AND (enthalpy<PIE)] AND (enthalpy<PIE)*/
                         geothermalflux_input->GetInputValue(&geothermalflux,gauss);
 
@@ -549 +580 @@
                         pressure_input->GetInputValue(&pressureup,gaussup);
                         if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
-                                // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
+                                // do nothing, set grad enthalpy*n=0.
                         }
                         else{
@@ -556 +587 @@
                 }
                 else{
-                        // base cold, but watercolumn positive. Set base to h_pmp.
+                        // base cold, but watercolumn positive. Set base to pressure melting point enthalpy
                 }
         }
@@ -810 +841 @@
         }
 }/*}}}*/
+
+
+
+
+
+
 void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
-
-}/*}}}*/
-void EnthalpyAnalysis::DrainWaterfraction(Element* element){/*{{{*/
-
-}/*}}}*/
+        /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
+        /* melting rate is positive when melting, negative when refreezing*/
+
+        /* Intermediaries */
+        int         i,is,vertexdown,vertexup,dim=3,numvertices,numsegments;
+        IssmDouble  heatflux;
+        IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
+        IssmDouble  temperature, waterfraction;
+        IssmDouble  basalfriction,alpha2;
+        IssmDouble  dt,yts;
+        IssmDouble  melting_overshoot,lambda;
+        IssmDouble  geothermalflux;
+        IssmDouble  vx,vy,vz;
+        IssmDouble *xyz_list      = NULL;
+        IssmDouble *xyz_list_base = NULL;
+        int        *pairindices   = NULL;
+
+        /* Only compute melt rates at the base of grounded ice*/
+        if(!element->IsOnBed() || element->IsFloating()) return;
+
+        /*Fetch parameters and inputs */
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
+        IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        //Input* watercolumn_input      = inputs->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
+        Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
+        //  Input* pressure_input         = inputs->GetInput(PressureEnum);                    _assert_(pressure_input);
+        //      Input* basalmeltingrate_input = inputs->GetInput(BasalforcingsMeltingRateEnum);    _assert_(basalmeltingrate_input);
+        Input* geothermalflux_input   = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+        Input* vx_input               = element->GetInput(VxEnum);                          _assert_(vx_input);
+        Input* vy_input               = element->GetInput(VyEnum);                          _assert_(vy_input);
+        Input* vz_input               = element->GetInput(VzEnum);                          _assert_(vz_input);
+        IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>0.);
+        element->NormalBase(&normal_base[0],xyz_list_base);
+        element->VerticalSegmentIndices(&pairindices,&numsegments);
+        IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
+        IssmDouble* heating = xNew<IssmDouble>(numsegments);    
+
+        /*Build friction element, needed later: */
+        Friction* friction=new Friction(element,dim);
+
+        /******** MELTING RATES  ************************************/
+        numvertices=element->GetNumberOfVertices();
+        IssmDouble* enthalpy = xNew<IssmDouble>(numvertices);
+        IssmDouble* pressure = xNew<IssmDouble>(numvertices);
+        IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
+        IssmDouble* basalmeltingrate = xNew<IssmDouble>(numvertices);
+        element->GetInputListOnVertices(enthalpy,EnthalpyEnum);
+        element->GetInputListOnVertices(pressure,PressureEnum);
+        element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
+        element->GetInputListOnVertices(basalmeltingrate,BasalforcingsMeltingRateEnum);
+
+        Gauss* gauss=element->NewGauss();
+        
+        for(int is=0;is<numsegments;is++){
+                vertexdown = pairindices[is*2+0];
+                vertexup   = pairindices[is*2+1];
+                gauss->GaussVertex(vertexdown);
+                
+                bool checkpositivethickness=true;
+                _assert_(watercolumn>=0.);
+
+                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
+                meltingrate_enthalpy[is]=0.;
+                heating[is]=0.;
+                if((watercolumn[vertexdown]>0.) && (enthalpy[vertexdown]<PureIceEnthalpy(element,pressure[vertexdown]))){
+                        /*ensure that no ice is at T<Tm(p), if water layer present*/
+                        enthalpy[vertexdown]=element->PureIceEnthalpy(pressure[vertexdown]); 
+                }
+                else if(enthalpy[vertexdown]<element->PureIceEnthalpy(pressure[vertexdown])){
+                        /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
+                        checkpositivethickness=false; // cold base, skip next test
+                }
+                else{/*we have a temperate base, go to next test*/}
+
+                if(checkpositivethickness){
+                        /*From here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
+                        bool istemperatelayer=false;
+                        if(enthalpy[vertexup]>=element->PureIceEnthalpy(pressure[vertexup])) istemperatelayer=true;
+                        if(istemperatelayer) for(i=0;i<dim;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
+                        else{
+                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
+                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
+                        }
+
+                        /*heat flux along normal*/
+                        heatflux=0.;
+                        for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
+
+                        /*basal friction*/
+                        friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+                        vx_input->GetInputValue(&vx,gauss);
+                        vy_input->GetInputValue(&vy,gauss);
+                        vz_input->GetInputValue(&vz,gauss);
+                        basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
+
+                        element->EnthalpyToThermal(&temperature,&waterfraction,enthalpy[vertexdown],pressure[vertexdown]);
+                        geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+                        // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
+                        heating[is]=(heatflux+basalfriction+geothermalflux);
+                        meltingrate_enthalpy[is]=heating[is]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent //????
+                }
+        }
+        // enthalpy might have been changed, update 
+        //element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
+
+        /******** DRAINAGE *****************************************/
+        IssmDouble* drainrate_column = xNew<IssmDouble>(numsegments); //TODO: xDelete?
+        IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
+        for(is=0;is<numsegments;is++)   drainrate_column[is]=0.;
+        Element* elementi = element;
+        for(;;){
+                for(is=0;is<numsegments;is++)   drainrate_element[is]=0.;
+                DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
+                for(is=0;is<numsegments;is++)   drainrate_column[is]+=drainrate_element[is];
+
+                if(elementi->IsOnSurface()) break;
+                elementi=elementi->GetUpperElement();                   
+        }
+        // add drained water to melting rate
+        for(is=0;is<numsegments;is++) meltingrate_enthalpy[is]+=drainrate_column[is];
+
+        /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        for(is=0;is<numsegments;is++){
+                vertexdown = pairindices[is*2+0];
+                vertexup   = pairindices[is*2+1];
+                if(reCast<bool,IssmDouble>(dt)){
+                        if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){     // prevent too much freeze on                   
+                                melting_overshoot=watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt;
+                                lambda=melting_overshoot/(meltingrate_enthalpy[is]*dt); _assert_(lambda>0); _assert_(lambda<1);
+                                basalmeltingrate[vertexdown]=(1.-lambda)*meltingrate_enthalpy[is];
+                                watercolumn[vertexdown]=0.;
+                                yts=365*24*60*60;
+                                enthalpy[vertexdown]+=dt/yts*lambda*heating[is];
+                        }
+                        else{
+                                basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
+                                watercolumn[vertexdown]+=dt*meltingrate_enthalpy[is]; 
+                        }
+                }
+                else{
+                        basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
+                        watercolumn[vertexdown]+=meltingrate_enthalpy[is];
+                }               
+                _assert_(watercolumn[vertexdown]>=0.);
+        }
+
+        /*feed updated variables back into model*/
+        element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
+        element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
+        element->AddInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum);
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+        xDelete<IssmDouble>(enthalpy);
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(watercolumn);
+        xDelete<IssmDouble>(basalmeltingrate);
+        xDelete<IssmDouble>(meltingrate_enthalpy);
+        xDelete<IssmDouble>(heating);
+        xDelete<IssmDouble>(drainrate_column);
+        xDelete<IssmDouble>(drainrate_element);
+}/*}}}*/
+
+void EnthalpyAnalysis::DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element){/*{{{*/
+
+        /*Intermediaries*/
+        int iv,is,vertexdown,vertexup,numsegments;      
+        IssmDouble dt, height_element;
+        IssmDouble rho_water, rho_ice;
+        int numvertices = element->GetNumberOfVertices();
+
+        IssmDouble* xyz_list = NULL;
+        IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
+        IssmDouble* pressures = xNew<IssmDouble>(numvertices);
+        IssmDouble* temperatures = xNew<IssmDouble>(numvertices);
+        IssmDouble* waterfractions = xNew<IssmDouble>(numvertices);
+        IssmDouble* deltawaterfractions = xNew<IssmDouble>(numvertices);
+        int        *pairindices   = NULL;
+        
+        rho_ice=element->GetMaterialParameter(MaterialsRhoIceEnum);
+        rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
+
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
+        element->GetInputListOnVertices(pressures,PressureEnum);
+
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        for(iv=0;iv<numvertices;iv++){ 
+                element->EnthalpyToThermal(&temperatures[iv],&waterfractions[iv], enthalpies[iv],pressures[iv]); 
+                deltawaterfractions[iv]=DrainageFunctionWaterfraction(waterfractions[iv], dt);
+        }
+        
+        /*drain waterfraction, feed updated variables back into model*/
+        for(iv=0;iv<numvertices;iv++){
+                if(reCast<bool,IssmDouble>(dt))
+                        waterfractions[iv]-=deltawaterfractions[iv]*dt;
+                else
+                        waterfractions[iv]-=deltawaterfractions[iv];
+                element->ThermalToEnthalpy(&enthalpies[iv], temperatures[iv], waterfractions[iv], pressures[iv]);
+        }
+        element->AddInput(EnthalpyEnum,enthalpies,P1Enum);
+        element->AddInput(WaterfractionEnum,waterfractions,P1Enum);
+
+        /*return meltwater column equivalent to drained water*/
+        element->VerticalSegmentIndices(&pairindices,&numsegments);
+        for(is=0;is<numsegments;is++){
+                vertexdown = pairindices[is*2+0];
+                vertexup   = pairindices[is*2+1];
+                height_element=fabs(xyz_list[vertexup*3+2]-xyz_list[vertexdown*3+2]);
+                pdrainrate_element[is]=(deltawaterfractions[vertexdown]+deltawaterfractions[vertexup])/2.*rho_water/rho_ice*height_element;
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(enthalpies);
+        xDelete<IssmDouble>(pressures);
+        xDelete<IssmDouble>(temperatures);
+        xDelete<IssmDouble>(waterfractions);
+        xDelete<IssmDouble>(deltawaterfractions);
+}/*}}}*/
+
+
+
+
+
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
 
-}/*}}}*/
+        /* /\*Intermediary*\/ */
+        /* bool        isdynamicbasalspc,setspc; */
+        /* int         numindices, numindicesup; */
+        /* IssmDouble  pressure, pressureup; */
+        /* IssmDouble  h_pmp, enthalpy, enthalpyup; */
+        /* IssmDouble  watercolumn; */
+        /* int        *indices = NULL, *indicesup = NULL; */
+        
+        /* /\* Only update Constraints at the base of grounded ice*\/ */
+        /* if(!IsOnBed() || IsFloating()) return; */
+
+        /* /\*Check wether dynamic basal boundary conditions are activated *\/ */
+        /* element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); */
+        /* if(!isdynamicbasalspc) return; */
+
+        /* /\*Fetch indices of basal & surface nodes for this finite element*\/ */
+        /* // BasalNodeIndices(&numindices,&indices,this->element_type); */
+        /* // SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); */
+        /* //   _assert_(numindices==numindicesup); */
+
+        /* /\*Get parameters and inputs: *\/ */
+        /* Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); */
+        /* Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); */
+        /* Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); */
+
+        /* /\*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*\/ */
+        /* GaussPenta* gauss=new GaussPenta(); */
+        /* GaussPenta* gaussup=new GaussPenta(); */
+        /* for(int i=0;i<numindices;i++){ */
+        /*      gauss->GaussNode(this->element_type,indices[i]); */
+        /*      gaussup->GaussNode(this->element_type,indicesup[i]); */
+
+        /*      /\*Check wether there is a temperate layer at the base or not *\/ */
+        /*      /\*check if node is temperate, if not, continue*\/ */
+        /*      enthalpy_input->GetInputValue(&enthalpy, gauss); */
+        /*      pressure_input->GetInputValue(&pressure, gauss); */
+        /*      watercolumn_input->GetInputValue(&watercolumn,gauss); */
+        /*      setspc = false; */
+        /*      // TODO: add case H<Hpmp && watercolumn>0.; */
+        /*      if (enthalpy>=element->PureIceEnthalpy(pressure)){ */
+        /*              /\*check if upper node is temperate, too. */
+        /*                      if yes, then we have a temperate layer of positive thickness and reset the spc. */
+        /*                      if not, apply dirichlet BC.*\/ */
+        /*              enthalpy_input->GetInputValue(&enthalpyup, gaussup); */
+        /*              pressure_input->GetInputValue(&pressureup, gaussup); */
+        /*              setspc=((enthalpyup<element->PureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; */
+        /*      } */
+
+        /*      if (setspc) { */
+        /*              /\*Calculate enthalpy at pressure melting point *\/ */
+        /*              h_pmp=matpar->PureIceEnthalpy(pressure); */
+        /*              /\*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*\/ */
+        /*              //element->ApplyConstraint(indices[i],1,h_pmp); */
+        /*              //nodes[indices[i]]->ApplyConstraint(1,h_pmp); */
+        /*      } */
+        /*      else { */
+        /*              /\*remove spc*\/ */
+        /*              //element->RemoveConstraint(indices[i],0); */
+        /*              //nodes[indices[i]]->DofInFSet(0); */
+        /*      } */
+        /* } */
+
+        /* /\*Free ressources:*\/ */
+        /* xDelete<int>(indices); */
+        /* xDelete<int>(indicesup); */
+        /* delete gauss; */
+        /* delete gaussup; */
+}/*}}}*/
+
+
+
+
 
 /*Intermediaries*/

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

@@ -40 +40 @@
                 static void PostProcessing(FemModel* femmodel);
                 static void ComputeBasalMeltingrate(Element* element);
-                static void DrainWaterfraction(Element* element);
+                static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
                 static void UpdateBasalConstraints(Element* element);
 
                 /*Intermediaries*/
-                IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
-                IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
-                IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
-                IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
+                static IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
+                static IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
+                static IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
+                static IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
 };
 #endif

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -105 +105 @@
                 virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
                 virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
+                virtual void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0;
                 virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
                 virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
@@ -123 +124 @@
                 virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
 
+                virtual Element* GetUpperElement(void)=0;
+                virtual Element* GetLowerElement(void)=0;
+                virtual Element* GetSurfaceElement(void)=0;
                 virtual Element* GetBasalElement(void)=0;
                 virtual void    GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -134 +134 @@
                                 penta->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
                                 if (penta->IsOnSurface()) break;
-                                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+                                penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
                         }
                 }
@@ -148 +148 @@
 }
 /*}}}*/
-
 /*FUNCTION Penta::BasalFrictionCreateInput {{{*/
 void Penta::BasalFrictionCreateInput(void){
@@ -479 +478 @@
 }
 /*}}}*/
+/*FUNCTION Penta::ThermalToEnthalpy{{{*/
+void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){
+        matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
+}
+/*}}}*/
 /*FUNCTION Penta::EnthalpyToThermal{{{*/
 void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
@@ -542 +546 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetBasalElement{{{*/
-Element* Penta::GetBasalElement(void){
+/*FUNCTION Penta::GetUpperPenta{{{*/
+Penta* Penta::GetUpperPenta(void){
+
+        Penta* upper_penta=NULL;
+
+        upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
+
+        return upper_penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetLowerPenta{{{*/
+Penta* Penta::GetLowerPenta(void){
+
+        Penta* lower_penta=NULL;
+
+        lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
+
+        return lower_penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetSurfacePenta{{{*/
+Penta* Penta::GetSurfacePenta(void){
 
         /*Output*/
         Penta* penta=NULL;
 
-        /*Go through all elements till the bed is reached*/
+        /*Go through all pentas till the surface is reached*/
+        penta=this;
+        for(;;){
+                /*Stop if we have reached the surface, else, take upper penta*/
+                if (penta->IsOnSurface()) break;
+
+                /* get upper Penta*/
+                penta=penta->GetUpperPenta();
+                _assert_(penta->Id()!=this->id);
+        }
+
+        /*return output*/
+        return penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetBasalPenta{{{*/
+Penta* Penta::GetBasalPenta(void){
+
+        /*Output*/
+        Penta* penta=NULL;
+
+        /*Go through all pentas till the bed is reached*/
         penta=this;
         for(;;){
@@ -555 +600 @@
 
                 /* get lower Penta*/
-                penta=penta->GetLowerElement();
+                penta=penta->GetLowerPenta();
                 _assert_(penta->Id()!=this->id);
         }
@@ -561 +606 @@
         /*return output*/
         return penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetUpperElement{{{*/
+Element* Penta::GetUpperElement(void){
+
+        /*Output*/
+        Element* upper_element=this->GetUpperPenta();
+        return upper_element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetLowerElement{{{*/
+Element* Penta::GetLowerElement(void){
+
+        /*Output*/
+        Element* lower_element=this->GetLowerPenta();
+        return lower_element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetSurfaceElement{{{*/
+Element* Penta::GetSurfaceElement(void){
+
+        /*Output*/
+        Element* element=this->GetSurfacePenta();
+        return element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetBasalElement{{{*/
+Element* Penta::GetBasalElement(void){
+
+        /*Output*/
+        Element* element=this->GetBasalPenta();
+        return element;
 }
 /*}}}*/
@@ -833 +910 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetLowerElement{{{*/
-Penta* Penta::GetLowerElement(void){
-
-        Penta* upper_penta=NULL;
-
-        upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
-
-        return upper_penta;
-}
-/*}}}*/
 /*FUNCTION Penta::GetNodeIndex {{{*/
 int Penta::GetNodeIndex(Node* node){
@@ -1121 +1188 @@
 
         return tau_parameter;
-}
-/*}}}*/
-/*FUNCTION Penta::GetUpperElement{{{*/
-Penta* Penta::GetUpperElement(void){
-
-        Penta* upper_penta=NULL;
-
-        upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
-
-        return upper_penta;
 }
 /*}}}*/
@@ -1486 +1543 @@
 
                 /* get upper Penta*/
-                penta=penta->GetUpperElement();
+                penta=penta->GetUpperPenta();
                 _assert_(penta->Id()!=this->id);
 
@@ -1555 +1612 @@
         for(;;){
                 /* get upper Penta*/
-                penta=penta->GetUpperElement();
+                penta=penta->GetUpperPenta();
                 _assert_(penta->Id()!=this->id);
 
@@ -1778 +1835 @@
 
                 /* get upper Penta*/
-                penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+                penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
         }
 
@@ -3416 +3473 @@
 
                 if(penta->IsOnSurface()) break;
-                penta=penta->GetUpperElement();                 
+                penta=penta->GetUpperPenta();                   
         }
         // add drained water to melting rate
@@ -3460 +3517 @@
 void Penta::DrainWaterfraction(IssmDouble* drainrate_element){
 
-    /*Intermediaries*/
+  /*Intermediaries*/
         bool isenthalpy;
         int iv, index0;
@@ -4609 +4666 @@
                         
                         /* get upper Penta*/
-                        penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+                        penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
                 }
         }

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -72 +72 @@
                 void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
                 void   Delta18oParameterization(void);
+                void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
                 void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
+                Penta* GetUpperPenta(void);
+                Penta* GetLowerPenta(void);
+                Penta* GetSurfacePenta(void);
+                Penta* GetBasalPenta(void);
+                Element* GetUpperElement(void);
+                Element* GetLowerElement(void);
+                Element* GetSurfaceElement(void);
                 Element* GetBasalElement(void);
                 void     GetDofList(int** pdoflist,int approximation_enum,int setenum);
@@ -211 +219 @@
                 void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
                 Node*          GetNode(int node_number);
-                Penta*         GetUpperElement(void);
-                Penta*         GetLowerElement(void);
                 void           InputChangeName(int input_enum, int enum_type_old);
                 void             InputExtrude(int enum_type,int object_type);

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -71 +71 @@
                 void        Delta18oParameterization(void){_error_("not implemented yet");};
                 void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+                void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
                 IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
                 IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
                 int         FiniteElement(void);
+                Element*    GetUpperElement(void){_error_("not implemented yet");};
+                Element*    GetLowerElement(void){_error_("not implemented yet");};
+                Element*    GetSurfaceElement(void){_error_("not implemented yet");};
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
                 void        GetDofList(int** pdoflist,int approximation_enum,int setenum){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -69 +69 @@
                 void        Delta18oParameterization(void);
                 void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+                void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
                 void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
                 int         FiniteElement(void);
+                Element*    GetUpperElement(void){_error_("not implemented yet");};
+                Element*    GetLowerElement(void){_error_("not implemented yet");};
+                Element*    GetSurfaceElement(void){_error_("not implemented yet");};
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
                 void          GetDofList(int** pdoflist,int approximation_enum,int setenum);