Changeset 21718

Timestamp:

05/16/17 16:55:45 (8 years ago)

Author:

jbondzio

Message:

CHG: adapting computation of basal melting rate (enthalpy) to higher-order elements

File:

• issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (modified) (4 diffs)

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

@@ -316 +316 @@
         const int   dim=3;
         int         i,is,state;
-        int                     vertexdown,vertexup,numvertices,numsegments;
+        int                     nodedown,nodeup,numnodes,numsegments;
         int                     enthalpy_enum;
         IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
@@ -357 +357 @@
         /******** MELTING RATES  ************************************//*{{{*/
         element->NormalBase(&normal_base[0],xyz_list_base);
-        element->VerticalSegmentIndices(&pairindices,&numsegments);
+        element->VerticalSegmentIndicesBase(&pairindices,&numsegments);
         IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
         IssmDouble* heating = xNew<IssmDouble>(numsegments);    
 
-        numvertices=element->GetNumberOfVertices();
-        IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
-        IssmDouble* pressures = xNew<IssmDouble>(numvertices);
-        IssmDouble* watercolumns = xNew<IssmDouble>(numvertices);
-        IssmDouble* basalmeltingrates = xNew<IssmDouble>(numvertices);
-        element->GetInputListOnVertices(enthalpies,enthalpy_enum);
-        element->GetInputListOnVertices(pressures,PressureEnum);
-        element->GetInputListOnVertices(watercolumns,WatercolumnEnum);
-        element->GetInputListOnVertices(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
+        numnodes=element->GetNumberOfNodes();
+        IssmDouble* enthalpies = xNew<IssmDouble>(numnodes);
+        IssmDouble* pressures = xNew<IssmDouble>(numnodes);
+        IssmDouble* watercolumns = xNew<IssmDouble>(numnodes);
+        IssmDouble* basalmeltingrates = xNew<IssmDouble>(numnodes);
+        element->GetInputListOnNodes(enthalpies,enthalpy_enum);
+        element->GetInputListOnNodes(pressures,PressureEnum);
+        element->GetInputListOnNodes(watercolumns,WatercolumnEnum);
+        element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
 
         Gauss* gauss=element->NewGauss();
         for(is=0;is<numsegments;is++){
-                vertexdown = pairindices[is*2+0];
-                vertexup   = pairindices[is*2+1];
-                gauss->GaussVertex(vertexdown);
-
-                state=GetThermalBasalCondition(element, enthalpies[vertexdown], enthalpies[vertexup], pressures[vertexdown], pressures[vertexup], watercolumns[vertexdown], basalmeltingrates[vertexdown]);
+                nodedown = pairindices[is*2+0];
+                nodeup   = pairindices[is*2+1];
+                gauss->GaussNode(element->GetElementType(),nodedown);
+
+                state=GetThermalBasalCondition(element, enthalpies[nodedown], enthalpies[nodeup], pressures[nodedown], pressures[nodeup], watercolumns[nodedown], basalmeltingrates[nodedown]);
                 switch (state) {
                         case 0:
@@ -392 +392 @@
                         case 4:
                                 // temperate, thick melting base: set grad H*n=0
-                                kappa_mix=GetWetIceConductivity(element, enthalpies[vertexdown], pressures[vertexdown]);
+                                kappa_mix=GetWetIceConductivity(element, enthalpies[nodedown], pressures[nodedown]);
                                 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
                                 for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
@@ -418 +418 @@
         /******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
         for(is=0;is<numsegments;is++){
-                vertexdown = pairindices[is*2+0];
-                vertexup   = pairindices[is*2+1];
+                nodedown = pairindices[is*2+0];
+                nodeup   = pairindices[is*2+1];
                 if(dt!=0.){
-                        if(watercolumns[vertexdown]+meltingrate_enthalpy[is]*dt<0.){    // prevent too much freeze on                   
-                                lambda = -watercolumns[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
-                                watercolumns[vertexdown]=0.;
-                                basalmeltingrates[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
-                                enthalpies[vertexdown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice; // use rest of energy to cool down base: dE=L*m, m=(1-lambda)*meltingrate*rho_ice
+                        if(watercolumns[nodedown]+meltingrate_enthalpy[is]*dt<0.){      // prevent too much freeze on                   
+                                lambda = -watercolumns[nodedown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
+                                watercolumns[nodedown]=0.;
+                                basalmeltingrates[nodedown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
+                                enthalpies[nodedown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice; // use rest of energy to cool down base: dE=L*m, m=(1-lambda)*meltingrate*rho_ice
                         }
                         else{
-                                basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
-                                watercolumns[vertexdown]+=dt*meltingrate_enthalpy[is]; 
+                                basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
+                                watercolumns[nodedown]+=dt*meltingrate_enthalpy[is]; 
                         }
                 }
                 else{
-                        basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
-                        if(watercolumns[vertexdown]+meltingrate_enthalpy[is]<0.)
-                                watercolumns[vertexdown]=0.;
+                        basalmeltingrates[nodedown]=meltingrate_enthalpy[is];
+                        if(watercolumns[nodedown]+meltingrate_enthalpy[is]<0.)
+                                watercolumns[nodedown]=0.;
                         else
-                                watercolumns[vertexdown]+=meltingrate_enthalpy[is];
+                                watercolumns[nodedown]+=meltingrate_enthalpy[is];
                 }       
-                basalmeltingrates[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
-                _assert_(watercolumns[vertexdown]>=0.);
+                basalmeltingrates[nodedown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
+                _assert_(watercolumns[nodedown]>=0.);
         }/*}}}*/
 
         /*feed updated variables back into model*/
         if(dt!=0.){
-                //element->AddInput(enthalpy_enum,enthalpies,P1Enum); //TODO: distinguis for steadystate and transient run
-                element->AddInput(WatercolumnEnum,watercolumns,P1Enum);
-        }
-        element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1Enum);
+                element->AddInput(enthalpy_enum,enthalpies,element->GetElementType()); 
+                element->AddInput(WatercolumnEnum,watercolumns,element->GetElementType());
+        }
+        element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,element->GetElementType());
 
         /*Clean up and return*/