Changeset 17028

Timestamp:

12/19/13 05:36:17 (11 years ago)

Author:

jbondzio

Message:

DEL: UpdateBasalConstraintsEnthalpy, ComputeBasalMeltingrate and DrainWaterfraction from Element and its children. DEL: module PostprocessingEnthalpyx

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• modules/PostprocessingEnthalpyx (deleted)
• modules/modules.h (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -391 +391 @@
 #}}}
 #Thermal sources  {{{
-thermal_sources = ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\
-                                           ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\
-                                                ./analyses/ThermalAnalysis.h\
+thermal_sources = ./analyses/ThermalAnalysis.h\
                                                 ./analyses/ThermalAnalysis.cpp\
                                                 ./analyses/EnthalpyAnalysis.h\

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

@@ -271 +271 @@
                 virtual void UpdateConstraintsExtrudeFromTop(void)=0;
 
-                #ifdef _HAVE_THERMAL_
-                virtual void UpdateBasalConstraintsEnthalpy(void)=0;
-                virtual void ComputeBasalMeltingrate(void)=0;
-                virtual void DrainWaterfraction(IssmDouble* drainrate_element)=0;
-                #endif
-
                 #ifdef _HAVE_HYDROLOGY_
                 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0;

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

@@ -3282 +3282 @@
 #endif
 
-#ifdef _HAVE_THERMAL_
-/*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/
-void  Penta::UpdateBasalConstraintsEnthalpy(void){
-
-        /*Intermediary*/
-        bool        isenthalpy,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 */
-        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
-        if(!isenthalpy) return;
-        parameters->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;
-                if (enthalpy>=matpar->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<matpar->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)*/
-                        nodes[indices[i]]->ApplyConstraint(1,h_pmp);
-                }
-                else {
-                        /*remove spc*/
-                        nodes[indices[i]]->DofInFSet(0);
-                }
-        }
-
-        /*Free ressources:*/
-        xDelete<int>(indices);
-        xDelete<int>(indicesup);
-        delete gauss;
-        delete gaussup;
-}
-/*}}}*/
-/*FUNCTION Penta::ComputeBasalMeltingrate{{{*/
-void Penta::ComputeBasalMeltingrate(void){
-        /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
-        /* melting rate is positive when melting, negative when refreezing*/
-
-        /* Intermediaries */
-        bool        isenthalpy, checkpositivethickness, istemperatelayer;
-        int         i,j,iv,analysis_type;
-        IssmDouble  xyz_list[NUMVERTICES][3];
-        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
-        IssmDouble  heatflux,kappa;
-        IssmDouble  vec_heatflux[3];
-        IssmDouble  normal_base[3], d1enthalpy[3];
-        IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES];
-        IssmDouble  enthalpy[NUMVERTICES],pressure[NUMVERTICES];
-        IssmDouble  temperature, waterfraction;
-        IssmDouble  latentheat, rho_ice;
-        IssmDouble  basalfriction, alpha2;
-        IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
-        IssmDouble  geothermalflux[NUMVERTICES];
-        IssmDouble  dt, yts;
-        IssmDouble  melting_overshoot,meltingrate_enthalpy[NUMVERTICES2D];
-        IssmDouble  drainrate_element[NUMVERTICES2D],drainrate_column[NUMVERTICES2D];
-        IssmDouble  lambda,heating[NUMVERTICES2D];
-        Friction   *friction  = NULL;
-        Penta      *penta = NULL;
-
-        /* Only compute melt rates at the base of grounded ice*/
-        if(!IsOnBed() || IsFloating()) return;
-
-        /*Check wether enthalpy is activated*/
-        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
-        if(!isenthalpy) return;
-
-        /*Fetch parameters and inputs */
-        latentheat=matpar->GetLatentHeat();
-        rho_ice=matpar->GetRhoIce();
-        GetInputListOnVertices(&vx[0],VxEnum);
-        GetInputListOnVertices(&vy[0],VyEnum);
-        GetInputListOnVertices(&vz[0],VzEnum);
-        GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
-        GetInputListOnVertices(&pressure[0],PressureEnum);
-        GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
-        GetInputListOnVertices(&basalmeltingrate[0],BasalforcingsMeltingRateEnum);
-        GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum);
-        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
-        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
-
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
-
-        /*Build friction element, needed later: */
-        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        friction=new Friction(this,3);
-
-        /******** MELTING RATES  ************************************/
-        GaussPenta* gauss=new GaussPenta();
-        for(iv=0;iv<NUMVERTICES2D;iv++){
-
-                gauss->GaussVertex(iv);
-                checkpositivethickness=true;
-
-                _assert_(watercolumn[iv]>=0.);
-
-                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
-                meltingrate_enthalpy[iv]=0.;
-                heating[iv]=0.;
-                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
-                        /*ensure that no ice is at T<Tm(p), if water layer present*/
-                        enthalpy[iv]=matpar->PureIceEnthalpy(pressure[iv]); 
-                }
-                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv])){
-                        /*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){
-                        /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
-                        istemperatelayer=false;
-                        if(enthalpy[iv+NUMVERTICES2D]>=matpar->PureIceEnthalpy(pressure[iv+NUMVERTICES2D])) istemperatelayer=true;
-                        if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
-                        else{
-                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
-                                kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,&enthalpy[0],&pressure[0]); _assert_(kappa>0.);
-                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
-                        }
-
-                        /*geothermal heatflux*/
-                        NormalBase(&normal_base[0],&xyz_list_tria[0][0]);
-                        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);
-                        basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0));
-
-                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
-                        // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
-                        heating[iv]=(heatflux+basalfriction+geothermalflux[iv]);
-                        meltingrate_enthalpy[iv]=heating[iv]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent 
-                }               
-        }
-        // enthalpy might have been changed, update
-        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
-
-        /******** DRAINAGE *****************************************/
-        for(iv=0; iv<NUMVERTICES2D; iv++)       
-                drainrate_column[iv]=0.;
-        penta=this;
-
-        for(;;){
-                for(iv=0; iv<NUMVERTICES2D; iv++)       drainrate_element[iv]=0.;
-                penta->DrainWaterfraction(&drainrate_element[0]); // TODO: make sure every vertex is only drained once
-                for(iv=0; iv<NUMVERTICES2D; iv++)       drainrate_column[iv]+=drainrate_element[iv];
-
-                if(penta->IsOnSurface()) break;
-                penta=penta->GetUpperPenta();                   
-        }
-        // add drained water to melting rate
-        for(iv=0; iv<NUMVERTICES2D;iv++)
-                meltingrate_enthalpy[iv]+=drainrate_column[iv];
-
-        /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
-        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        for(iv=0;iv<NUMVERTICES2D;iv++){
-                if(reCast<bool,IssmDouble>(dt)){
-                        if(watercolumn[iv]+meltingrate_enthalpy[iv]*dt<0.){                             
-                                melting_overshoot=watercolumn[iv]+meltingrate_enthalpy[iv]*dt;
-                                lambda=melting_overshoot/(meltingrate_enthalpy[iv]*dt); _assert_(lambda>0); _assert_(lambda<1);
-                                basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy[iv];
-                                watercolumn[iv]=0.;
-                                yts=365*24*60*60;
-                                enthalpy[iv]+=dt/yts*lambda*heating[iv];
-                        }
-                        else{
-                                basalmeltingrate[iv]=meltingrate_enthalpy[iv];
-                                watercolumn[iv]+=dt*meltingrate_enthalpy[iv]; 
-                        }
-                }
-                else{
-                        basalmeltingrate[iv]=meltingrate_enthalpy[iv];
-                        watercolumn[iv]+=meltingrate_enthalpy[iv];
-                }         
-                
-                _assert_(watercolumn[iv]>=0.);
-        }
-
-        /*feed updated variables back into model*/
-        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
-        this->inputs->AddInput(new PentaInput(WatercolumnEnum,watercolumn,P1Enum));
-        this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum));
-
-        /*Clean up and return*/
-        delete gauss;
-        delete friction;
-}
-/*}}}*/
-/*FUNCTION Penta::DrainWaterfraction{{{*/
-void Penta::DrainWaterfraction(IssmDouble* drainrate_element){
-
-  /*Intermediaries*/
-        bool isenthalpy;
-        int iv, index0;
-        IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
-        IssmDouble dw[NUMVERTICES];
-        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; 
-        IssmDouble dt, height_element;
-        IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble rho_water, rho_ice;
-
-        /*Check wether enthalpy is activated*/
-        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
-        if(!isenthalpy) return;       
-        
-        rho_ice=matpar->GetRhoIce();
-        rho_water=matpar->GetRhoFreshwater();
-
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
-        this->GetInputListOnVertices(&pressure[0],PressureEnum);
-
-        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        for(iv=0;iv<NUMVERTICES;iv++){ 
-                matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); 
-                dw[iv]=DrainageFunctionWaterfraction(waterfraction[iv], dt);
-        }
-        
-        /*drain waterfraction, feed updated variables back into model*/
-        for(iv=0;iv<NUMVERTICES;iv++){
-                if(reCast<bool,IssmDouble>(dt))
-                        waterfraction[iv]-=dw[iv]*dt;
-                else
-                        waterfraction[iv]-=dw[iv];
-                matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);
-        }
-        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
-        this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
-
-        /*return meltwater column equivalent to drained water*/
-        for(iv=0;iv<NUMVERTICES2D;iv++){
-                index0=(iv+NUMVERTICES2D);
-                height_element=fabs(xyz_list[index0][2]-xyz_list[iv][2]);
-                drainrate_element[iv]=(dw[iv]+dw[index0])/2.*rho_water/rho_ice*height_element;
-        }
-}
-/*}}}*/
-#endif
-
 #ifdef _HAVE_CONTROL_
 /*FUNCTION Penta::ControlInputGetGradient{{{*/

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

@@ -269 +269 @@
                 void           UpdateConstraintsExtrudeFromBase(void){_error_("not implemented yet");};
                 void           UpdateConstraintsExtrudeFromTop(void){_error_("not implemented yet");};
-                #ifdef _HAVE_THERMAL_
-                void           UpdateBasalConstraintsEnthalpy(void);
-                void           ComputeBasalMeltingrate(void);
-                void           DrainWaterfraction(IssmDouble* drainrate_element);
-                #endif
                 /*}}}*/
 };

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

@@ -159 +159 @@
                 void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
 
-                #ifdef _HAVE_THERMAL_
-                void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
-                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
-                void DrainWaterfraction(IssmDouble* drainrate_element){_error_("not implemented yet");};
-                #endif
                 #ifdef _HAVE_HYDROLOGY_
                 void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};

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

@@ -266 +266 @@
                 void UpdateConstraintsExtrudeFromBase(void);
                 void UpdateConstraintsExtrudeFromTop(void);
-                #ifdef _HAVE_THERMAL_
-                void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
-                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
-                void DrainWaterfraction(IssmDouble* drainrate_element){_error_("not implemented yet");};
-                #endif
 
                 #ifdef _HAVE_HYDROLOGY_

issm/trunk-jpl/src/c/modules/modules.h

@@ -76 +76 @@
 #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
 #include "./PositiveDegreeDayx/PositiveDegreeDayx.h"
-#include "./PostprocessingEnthalpyx/PostprocessingEnthalpyx.h"
 #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
 #include "./Reduceloadx/Reduceloadx.h"