Changeset 18612


Ignore:
Timestamp:
10/10/14 06:23:20 (10 years ago)
Author:
jbondzio
Message:

CHG: new unified decision chart for basal thermal condition

Location:
issm/trunk-jpl
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r18591 r18612  
    207207        InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
    208208
    209         int solution_type;
    210         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    211         if(solution_type!=SteadystateSolutionEnum){
    212                 PostProcessing(femmodel);
    213         }
     209        IssmDouble dt;
     210        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     211        if(dt==0.) ComputeBasalMeltingrate(femmodel);
     212        else PostProcessing(femmodel);
    214213
    215214}/*}}}*/
    … …  
    563562        if(!element->IsOnBase() || element->IsFloating()) return NULL;
    564563
    565         IssmDouble  dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
    566         IssmDouble  enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
     564        int i, state;
     565        IssmDouble  dt,Jdet,scalar;
     566        IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
     567        IssmDouble      vx,vy,vz;
     568        IssmDouble  alpha2,basalfriction,geothermalflux,heatflux;
    567569        IssmDouble *xyz_list_base = NULL;
    568570
    … …  
    580582        Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
    581583        Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
    582         Input* enthalpy_input       = element->GetInput(EnthalpyPicardEnum);              _assert_(enthalpy_input);
    583         Input* pressure_input       = element->GetInput(PressureEnum);                    _assert_(pressure_input);
     584        Input* enthalpy_input            = element->GetInput(EnthalpyPicardEnum);                                        _assert_(enthalpy_input);
     585        Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
     586        Input* watercolumn_input         = element->GetInput(WatercolumnEnum);                                                   _assert_(watercolumn_input);
     587        Input* meltingrate_input         = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);                                                   _assert_(meltingrate_input);
    584588        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    585         Input* watercolumn_input    = element->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
    586         IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
    587         IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     589        IssmDouble  rho_ice                      = element->GetMaterialParameter(MaterialsRhoIceEnum);
    588590
    589591        /*Build friction element, needed later: */
    … …  
    591593
    592594        /* Start  looping on the number of gaussian points: */
    593         Gauss* gauss   = element->NewGaussBase(2);
    594         Gauss* gaussup = element->NewGaussTop(2);
     595        GaussPenta* gauss=new GaussPenta(2,2);
     596        GaussPenta* gaussup=new GaussPenta(2,2);
    595597        for(int ig=gauss->begin();ig<gauss->end();ig++){
    596598                gauss->GaussPoint(ig);
    … …  
    600602
    601603                enthalpy_input->GetInputValue(&enthalpy,gauss);
     604                enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    602605                pressure_input->GetInputValue(&pressure,gauss);
     606                pressure_input->GetInputValue(&pressureup,gaussup);
    603607                watercolumn_input->GetInputValue(&watercolumn,gauss);
    604 
    605                 if((dt==0.) || ((watercolumn<=0.) && (enthalpy<PureIceEnthalpy(element,pressure)))){
    606                         /* the above check is equivalent to
    607                          NOT [(watercolumn>0.) AND (enthalpy<PIE)] AND (enthalpy<PIE)*/
    608                         geothermalflux_input->GetInputValue(&geothermalflux,gauss);
    609 
    610                         friction->GetAlpha2(&alpha2,gauss);
    611                         vx_input->GetInputValue(&vx,gauss);
    612                         vy_input->GetInputValue(&vy,gauss);
    613                         vz_input->GetInputValue(&vz,gauss);
    614                         basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
    615                         heatflux      = (basalfriction+geothermalflux)/(rho_ice);
    616 
    617                         scalar = gauss->weight*Jdet*heatflux;
    618                         if(dt!=0.) scalar=dt*scalar;
    619 
    620                         for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    621                 }
    622                 else if(enthalpy >= PureIceEnthalpy(element,pressure)){
    623                         /* check positive thickness of temperate basal ice layer */
    624                         enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    625                         pressure_input->GetInputValue(&pressureup,gaussup);
    626                         if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
    627                                 // do nothing, set grad enthalpy*n=0.
    628                         }
    629                         else{
    630                                 // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
    631                         }
    632                 }
    633                 else{
    634                         // base cold, but watercolumn positive. Set base to pressure melting point enthalpy
     608                meltingrate_input->GetInputValue(&meltingrate,gauss);
     609
     610                state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
     611                switch (state) {
     612                        case 0:
     613                                // cold, dry base: apply basal surface forcing
     614                                geothermalflux_input->GetInputValue(&geothermalflux,gauss);
     615                                friction->GetAlpha2(&alpha2,gauss);
     616                                vx_input->GetInputValue(&vx,gauss);
     617                                vy_input->GetInputValue(&vy,gauss);
     618                                vz_input->GetInputValue(&vz,gauss);
     619                                basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
     620                                heatflux=(basalfriction+geothermalflux)/(rho_ice);
     621                                scalar=gauss->weight*Jdet*heatflux;
     622                                if(dt!=0.) scalar=dt*scalar;
     623                                for(i=0;i<numnodes;i++)
     624                                        pe->values[i]+=scalar*basis[i];
     625                                break;
     626                        case 1:
     627                                // cold, wet base: keep at pressure melting point
     628                        case 2:
     629                                // temperate, thin refreezing base: release spc
     630                        case 3:
     631                                // temperate, thin melting base: set spc
     632                        case 4:
     633                                // temperate, thick melting base: set grad H*n=0
     634                                for(i=0;i<numnodes;i++)
     635                                        pe->values[i]+=0.;
     636                                break;
     637                        default:
     638                                _printf0_("     unknown thermal basal state found!");
    635639                }
    636640        }
    … …  
    653657        if(!element->IsOnBase() || !element->IsFloating()) return NULL;
    654658
    655         IssmDouble  h_pmp,dt,Jdet,scalar_ocean,pressure;
     659        IssmDouble  Hpmp,dt,Jdet,scalar_ocean,pressure;
    656660        IssmDouble *xyz_list_base = NULL;
    657661
    … …  
    683687
    684688                pressure_input->GetInputValue(&pressure,gauss);
    685                 h_pmp=element->PureIceEnthalpy(pressure);
    686 
    687                 scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
     689                Hpmp=element->PureIceEnthalpy(pressure);
     690
     691                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*Hpmp/(heatcapacity*rho_ice);
    688692                if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
    689693
    … …  
    887891
    888892        /*Intermediaries*/
    889         int solution_type, i;
    890893        bool computebasalmeltingrates=true;
    891         bool isdrainage=true;
    892         bool updatebasalconstraints=true;
    893 
    894         if(isdrainage){
    895                 /*Drain excess water fraction in ice column: */
    896                 for(i=0;i<femmodel->elements->Size();i++){
    897                         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    898                         DrainWaterfractionIcecolumn(element);
    899                 }
    900         }
    901 
    902         if(computebasalmeltingrates){
    903                 /*Compute basal melting rates: */
    904                 for(i=0;i<femmodel->elements->Size();i++){
    905                         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    906                         ComputeBasalMeltingrate(element);
    907                 }
    908         }
    909 
    910         if(updatebasalconstraints){
    911                 /*Update basal dirichlet BCs for enthalpy in transient runs: */
    912                 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    913                 if(solution_type==TransientSolutionEnum){
    914                         for(i=0;i<femmodel->elements->Size();i++){
    915                                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    916                                 UpdateBasalConstraints(element);
    917                         }
    918                 }
     894        bool drainicecolumn=true;
     895        bool isdynamicbasalspc;
     896        IssmDouble dt;
     897
     898        femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
     899        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     900
     901        //TODO: use dt to decide what to do
     902        if(drainicecolumn)      DrainWaterfraction(femmodel);
     903        if(computebasalmeltingrates)    ComputeBasalMeltingrate(femmodel);
     904        if(isdynamicbasalspc)   UpdateBasalConstraints(femmodel);
     905
     906}/*}}}*/
     907void EnthalpyAnalysis::ComputeBasalMeltingrate(FemModel* femmodel){/*{{{*/
     908        /*Compute basal melting rates: */
     909        for(int i=0;i<femmodel->elements->Size();i++){
     910                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     911                ComputeBasalMeltingrate(element);
    919912        }
    920913}/*}}}*/
    … …  
    931924        /* Intermediaries */
    932925        const int   dim=3;
    933         int         i,is,vertexdown,vertexup,numvertices,numsegments;
    934         IssmDouble  heatflux;
    935         IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
    936         IssmDouble  basalfriction,alpha2;
     926        int         i,is,state;
     927        int                     vertexdown,vertexup,numvertices,numsegments;
     928        const int       enthalpy_enum=EnthalpyEnum;
     929        IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
     930        IssmDouble  basalfriction,alpha2,geothermalflux,heatflux;
    937931        IssmDouble  dt,yts;
    938932        IssmDouble  melting_overshoot,lambda;
    939         IssmDouble  geothermalflux;
    940933        IssmDouble  vx,vy,vz;
    941934        IssmDouble *xyz_list      = NULL;
    … …  
    943936        int        *pairindices   = NULL;
    944937
    945         /*Fetch parameters and inputs */
     938        /*Fetch parameters*/
    946939        element->GetVerticesCoordinates(&xyz_list);
    947940        element->GetVerticesCoordinatesBase(&xyz_list_base);
    … …  
    949942        IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
    950943        IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    951         Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
     944        IssmDouble beta          = element->GetMaterialParameter(MaterialsBetaEnum);
     945        IssmDouble kappa                 = EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>=0.);
     946        IssmDouble kappa_mix;
     947
     948        /*retrieve inputs*/
     949        Input* enthalpy_input         = element->GetInput(enthalpy_enum);                    _assert_(enthalpy_input);
     950        Input* pressure_input                   = element->GetInput(PressureEnum);                                                       _assert_(pressure_input);
    952951        Input* geothermalflux_input   = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    953952        Input* vx_input               = element->GetInput(VxEnum);                          _assert_(vx_input);
    954953        Input* vy_input               = element->GetInput(VyEnum);                          _assert_(vy_input);
    955954        Input* vz_input               = element->GetInput(VzEnum);                          _assert_(vz_input);
    956         IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>=0.);
     955
     956        /*Build friction element, needed later: */
     957        Friction* friction=new Friction(element,dim);
     958
     959        /******** MELTING RATES  ************************************//*{{{*/
    957960        element->NormalBase(&normal_base[0],xyz_list_base);
    958961        element->VerticalSegmentIndices(&pairindices,&numsegments);
    … …  
    960963        IssmDouble* heating = xNew<IssmDouble>(numsegments);   
    961964
    962         /*Build friction element, needed later: */
    963         Friction* friction=new Friction(element,dim);
    964 
    965         /******** MELTING RATES  ************************************/
    966965        numvertices=element->GetNumberOfVertices();
    967         IssmDouble* enthalpy = xNew<IssmDouble>(numvertices);
    968         IssmDouble* pressure = xNew<IssmDouble>(numvertices);
    969         IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
    970         IssmDouble* basalmeltingrate = xNew<IssmDouble>(numvertices);
    971         element->GetInputListOnVertices(enthalpy,EnthalpyEnum);
    972         element->GetInputListOnVertices(pressure,PressureEnum);
    973         element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
    974         element->GetInputListOnVertices(basalmeltingrate,BasalforcingsGroundediceMeltingRateEnum);
     966        IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
     967        IssmDouble* pressures = xNew<IssmDouble>(numvertices);
     968        IssmDouble* watercolumns = xNew<IssmDouble>(numvertices);
     969        IssmDouble* basalmeltingrates = xNew<IssmDouble>(numvertices);
     970        element->GetInputListOnVertices(enthalpies,enthalpy_enum);
     971        element->GetInputListOnVertices(pressures,PressureEnum);
     972        element->GetInputListOnVertices(watercolumns,WatercolumnEnum);
     973        element->GetInputListOnVertices(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
    975974
    976975        Gauss* gauss=element->NewGauss();
    977        
    978         for(int is=0;is<numsegments;is++){
     976        for(is=0;is<numsegments;is++){
    979977                vertexdown = pairindices[is*2+0];
    980978                vertexup   = pairindices[is*2+1];
    981979                gauss->GaussVertex(vertexdown);
    982                
    983                 bool checkpositivethickness=true;
    984                 _assert_(watercolumn[vertexdown]>=0.);
    985 
    986                 /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
    987                 meltingrate_enthalpy[is]=0.;
    988                 heating[is]=0.;
    989                 if((watercolumn[vertexdown]>0.) && (enthalpy[vertexdown]<PureIceEnthalpy(element,pressure[vertexdown]))){
    990                         /*ensure that no ice is at T<Tm(p), if water layer present*/
    991                         enthalpy[vertexdown]=element->PureIceEnthalpy(pressure[vertexdown]);
    992                 }
    993                 else if(enthalpy[vertexdown]<element->PureIceEnthalpy(pressure[vertexdown])){
    994                         /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
    995                         checkpositivethickness=false; // cold base, skip next test
    996                 }
    997                 else{/*we have a temperate base, go to next test*/}
    998 
    999                 if(checkpositivethickness){
    1000                         /*From here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
    1001                         bool istemperatelayer=false;
    1002                         if(enthalpy[vertexup]>=element->PureIceEnthalpy(pressure[vertexup])) istemperatelayer=true;
    1003                         if(istemperatelayer) for(i=0;i<dim;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
    1004                         else{
     980
     981                state=GetThermalBasalCondition(element, enthalpies[vertexdown], enthalpies[vertexup], pressures[vertexdown], pressures[vertexup], watercolumns[vertexdown], basalmeltingrates[vertexdown]);
     982                switch (state) {
     983                        case 0:
     984                                // cold, dry base: apply basal surface forcing
     985                                for(i=0;i<3;i++) vec_heatflux[i]=0.;
     986                                break;
     987                        case 1:
     988                                // cold, wet base: keep at pressure melting point
     989                        case 2:
     990                                // temperate, thin refreezing base: release spc
     991
     992                        case 3:
     993                                // temperate, thin melting base: set spc
     994                                // enthalpies[vertexdown]=element->PureIceEnthalpy(pressures[vertexdown]);
    1005995                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
    1006996                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    1007                         }
    1008 
     997                                break;
     998                        case 4:
     999                                // temperate, thick melting base: set grad H*n=0
     1000                                kappa_mix=GetWetIceConductivity(element, enthalpies[vertexdown], pressures[vertexdown]);
     1001                                pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
     1002                                for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
     1003                                break;
     1004                        default:
     1005                                _printf0_("     unknown thermal basal state found!");
     1006                }
     1007                if(state==0) meltingrate_enthalpy[is]=0.;
     1008                else{
    10091009                        /*heat flux along normal*/
    10101010                        heatflux=0.;
    … …  
    10131013                        /*basal friction*/
    10141014                        friction->GetAlpha2(&alpha2,gauss);
    1015                         vx_input->GetInputValue(&vx,gauss);
    1016                         vy_input->GetInputValue(&vy,gauss);
    1017                         vz_input->GetInputValue(&vz,gauss);
     1015                        vx_input->GetInputValue(&vx,gauss);             vy_input->GetInputValue(&vy,gauss);             vz_input->GetInputValue(&vz,gauss);
    10181016                        basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
    1019 
    10201017                        geothermalflux_input->GetInputValue(&geothermalflux,gauss);
    10211018                        /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
    … …  
    10231020                        meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
    10241021                }
    1025         }
    1026         /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
     1022        }/*}}}*/
     1023
     1024        /******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
    10271025        element->FindParam(&dt,TimesteppingTimeStepEnum);
     1026        element->FindParam(&yts, ConstantsYtsEnum);
    10281027        for(is=0;is<numsegments;is++){
    10291028                vertexdown = pairindices[is*2+0];
    10301029                vertexup   = pairindices[is*2+1];
    10311030                if(dt!=0.){
    1032                         if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){     // prevent too much freeze on                   
    1033                                 lambda = -watercolumn[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
    1034                                 watercolumn[vertexdown]=0.;
    1035                                 basalmeltingrate[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
    1036                                 enthalpy[vertexdown]+=(1.-lambda)*meltingrate_enthalpy[is]*dt*latentheat; // use rest of energy to cool down base
     1031                        if(watercolumns[vertexdown]+meltingrate_enthalpy[is]*dt<0.){    // prevent too much freeze on                   
     1032                                lambda = -watercolumns[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
     1033                                watercolumns[vertexdown]=0.;
     1034                                basalmeltingrates[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
     1035                                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
    10371036                        }
    10381037                        else{
    1039                                 basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
    1040                                 watercolumn[vertexdown]+=dt*meltingrate_enthalpy[is];
     1038                                basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
     1039                                watercolumns[vertexdown]+=dt*meltingrate_enthalpy[is];
    10411040                        }
    10421041                }
    10431042                else{
    1044                         basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
    1045                         if(watercolumn[vertexdown]+meltingrate_enthalpy[is]<0.)
    1046                                 watercolumn[vertexdown]=0.;
     1043                        basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
     1044                        if(watercolumns[vertexdown]+meltingrate_enthalpy[is]<0.)
     1045                                watercolumns[vertexdown]=0.;
    10471046                        else
    1048                                 watercolumn[vertexdown]+=meltingrate_enthalpy[is];
     1047                                watercolumns[vertexdown]+=meltingrate_enthalpy[is];
    10491048                }       
    1050                 basalmeltingrate[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
    1051                 _assert_(watercolumn[vertexdown]>=0.);
    1052         }
     1049                basalmeltingrates[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
     1050                _assert_(watercolumns[vertexdown]>=0.);
     1051        }/*}}}*/
    10531052
    10541053        /*feed updated variables back into model*/
    1055         element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
    1056         element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
    1057         element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrate,P1Enum);
     1054        element->AddInput(EnthalpyEnum,enthalpies,P1Enum); //TODO: distinguis for steadystate and transient run
     1055        element->AddInput(WatercolumnEnum,watercolumns,P1Enum);
     1056        element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1Enum);
    10581057
    10591058        /*Clean up and return*/
    … …  
    10611060        delete friction;
    10621061        xDelete<int>(pairindices);
    1063         xDelete<IssmDouble>(enthalpy);
    1064         xDelete<IssmDouble>(pressure);
    1065         xDelete<IssmDouble>(watercolumn);
    1066         xDelete<IssmDouble>(basalmeltingrate);
     1062        xDelete<IssmDouble>(enthalpies);
     1063        xDelete<IssmDouble>(pressures);
     1064        xDelete<IssmDouble>(watercolumns);
     1065        xDelete<IssmDouble>(basalmeltingrates);
    10671066        xDelete<IssmDouble>(meltingrate_enthalpy);
    10681067        xDelete<IssmDouble>(heating);
    10691068        xDelete<IssmDouble>(xyz_list);
    10701069        xDelete<IssmDouble>(xyz_list_base);
     1070}/*}}}*/
     1071void EnthalpyAnalysis::DrainWaterfraction(FemModel* femmodel){/*{{{*/
     1072        /*Drain excess water fraction in ice column: */
     1073        for(int i=0;i<femmodel->elements->Size();i++){
     1074                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     1075                DrainWaterfractionIcecolumn(element);
     1076        }
    10711077}/*}}}*/
    10721078void EnthalpyAnalysis::DrainWaterfractionIcecolumn(Element* element){/*{{{*/
    … …  
    11731179        xDelete<IssmDouble>(deltawaterfractions);
    11741180}/*}}}*/
     1181void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
     1182        /*Update basal dirichlet BCs for enthalpy: */
     1183        for(int i=0;i<femmodel->elements->Size();i++){
     1184                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     1185                UpdateBasalConstraints(element);
     1186        }
     1187}/*}}}*/
    11751188void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
    11761189
    … …  
    11821195
    11831196        /*Intermediary*/
    1184         bool        isdynamicbasalspc,setspc;
    1185         int         numindices, numindicesup;
    1186         IssmDouble  pressure, pressureup;
    1187         IssmDouble  h_pmp, enthalpy, enthalpyup;
    1188         IssmDouble  watercolumn;
    1189         int        *indices = NULL, *indicesup = NULL;
    1190         Node*       node = NULL;
     1197        bool        isdynamicbasalspc;
     1198        IssmDouble      dt;
    11911199
    11921200        /*Check wether dynamic basal boundary conditions are activated */
    … …  
    11941202        if(!isdynamicbasalspc) return;
    11951203
     1204        element->FindParam(&dt,TimesteppingTimeStepEnum);
     1205        if(dt==0.){
     1206                UpdateBasalConstraintsSteadystate(element);
     1207        }
     1208        else{
     1209                UpdateBasalConstraintsTransient(element);
     1210        }
     1211}/*}}}*/
     1212void EnthalpyAnalysis::UpdateBasalConstraintsTransient(Element* element){/*{{{*/
     1213
     1214        /* Check if ice in element */
     1215        if(!element->IsIceInElement()) return;
     1216
     1217        /* Only update Constraints at the base of grounded ice*/
     1218        if(!(element->IsOnBase()) || element->IsFloating()) return;
     1219
     1220        /*Intermediary*/
     1221        bool        isdynamicbasalspc,setspc;
     1222        int         numindices, numindicesup, state;
     1223        int        *indices = NULL, *indicesup = NULL;
     1224        Node*       node = NULL;
     1225        IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
     1226
     1227        /*Check wether dynamic basal boundary conditions are activated */
     1228        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
     1229        if(!isdynamicbasalspc) return;
     1230
     1231        /*Get parameters and inputs: */
     1232        Input* enthalpy_input       = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input); //TODO: check EnthalpyPicard?
     1233        Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
     1234        Input* watercolumn_input         = element->GetInput(WatercolumnEnum);                                                   _assert_(watercolumn_input);
     1235        Input* meltingrate_input         = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);                                                   _assert_(meltingrate_input);
     1236
    11961237        /*Fetch indices of basal & surface nodes for this finite element*/
    11971238        Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
    11981239        penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
    1199         penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());
    1200         _assert_(numindices==numindicesup);
     1240        penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());  _assert_(numindices==numindicesup);
     1241
     1242        GaussPenta* gauss=new GaussPenta();
     1243        GaussPenta* gaussup=new GaussPenta();
     1244
     1245        for(int i=0;i<numindices;i++){
     1246                gauss->GaussNode(element->GetElementType(),indices[i]);
     1247                gaussup->GaussNode(element->GetElementType(),indicesup[i]);
     1248               
     1249                enthalpy_input->GetInputValue(&enthalpy,gauss);
     1250                enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     1251                pressure_input->GetInputValue(&pressure,gauss);
     1252                pressure_input->GetInputValue(&pressureup,gaussup);
     1253                watercolumn_input->GetInputValue(&watercolumn,gauss);
     1254                meltingrate_input->GetInputValue(&meltingrate,gauss);
     1255
     1256                state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
     1257
     1258                setspc=false;
     1259                switch (state) {
     1260                        case 0:
     1261                                // cold, dry base: apply basal surface forcing
     1262                                break;
     1263                        case 1:
     1264                                // cold, wet base: keep at pressure melting point
     1265                                setspc=true;
     1266                                break;
     1267                        case 2:
     1268                                // temperate, thin refreezing base: release spc
     1269                                break;
     1270                        case 3:
     1271                                // temperate, thin melting base: set spc
     1272                                setspc=true;
     1273                                break;
     1274                        case 4:
     1275                                // temperate, thick melting base: set grad H*n=0
     1276                                break;
     1277                        default:
     1278                                _printf0_("     unknown thermal basal state found!");
     1279                }
     1280
     1281                /*apply or release spc*/
     1282                node=element->GetNode(indices[i]);
     1283                if(setspc){
     1284                        pressure_input->GetInputValue(&pressure, gauss);
     1285                        node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
     1286                }
     1287                else                   
     1288                        node->DofInFSet(0);
     1289        }
     1290
     1291        /*Free ressources:*/
     1292        xDelete<int>(indices);
     1293        xDelete<int>(indicesup);
     1294        delete gauss;
     1295        delete gaussup;
     1296}/*}}}*/
     1297void EnthalpyAnalysis::UpdateBasalConstraintsSteadystate(Element* element){/*{{{*/
     1298
     1299        /* Check if ice in element */
     1300        if(!element->IsIceInElement()) return;
     1301
     1302        /* Only update Constraints at the base of grounded ice*/
     1303        if(!(element->IsOnBase()) || element->IsFloating()) return;
     1304
     1305        /*Intermediary*/
     1306        bool        isdynamicbasalspc,setspc;
     1307        int         numindices, numindicesup, state;
     1308        int        *indices = NULL, *indicesup = NULL;
     1309        Node*       node = NULL;
     1310        IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
     1311
     1312        /*Check wether dynamic basal boundary conditions are activated */
     1313        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
     1314        if(!isdynamicbasalspc) return;
    12011315
    12021316        /*Get parameters and inputs: */
    1203         Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
    1204         Input* enthalpy_input=element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
    1205         Input* watercolumn_input=element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
    1206 
    1207         /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
     1317        Input* enthalpy_input            = element->GetInput(EnthalpyPicardEnum);                                        _assert_(enthalpy_input);
     1318        Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
     1319        Input* watercolumn_input         = element->GetInput(WatercolumnEnum);                                                   _assert_(watercolumn_input);
     1320        Input* meltingrate_input         = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);                                                   _assert_(meltingrate_input);
     1321
     1322        /*Fetch indices of basal & surface nodes for this finite element*/
     1323        Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
     1324        penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
     1325        penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());  _assert_(numindices==numindicesup);
     1326
    12081327        GaussPenta* gauss=new GaussPenta();
    12091328        GaussPenta* gaussup=new GaussPenta();
    … …  
    12121331                gaussup->GaussNode(element->GetElementType(),indicesup[i]);
    12131332
    1214                 /*Check wether there is a temperate layer at the base or not */
    1215                 /*check if node is temperate, else continue*/
    1216                 enthalpy_input->GetInputValue(&enthalpy, gauss);
    1217                 pressure_input->GetInputValue(&pressure, gauss);
     1333                enthalpy_input->GetInputValue(&enthalpy,gauss);
     1334                enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     1335                pressure_input->GetInputValue(&pressure,gauss);
     1336                pressure_input->GetInputValue(&pressureup,gaussup);
    12181337                watercolumn_input->GetInputValue(&watercolumn,gauss);
    1219                 h_pmp=PureIceEnthalpy(element,pressure);
    1220                 if (enthalpy>=h_pmp){
    1221                         /*check if upper node is temperate, too.
    1222                                 if yes, then we have a temperate layer of positive thickness and reset the spc.
    1223                                 if not, apply dirichlet BC.*/
    1224                         enthalpy_input->GetInputValue(&enthalpyup, gaussup);
    1225                         pressure_input->GetInputValue(&pressureup, gaussup);
    1226                         setspc=((enthalpyup<PureIceEnthalpy(element,pressureup)) && (watercolumn>=0.))?true:false;
    1227                 }
    1228                 else
    1229                         setspc = false;
    1230 
     1338                meltingrate_input->GetInputValue(&meltingrate,gauss);
     1339
     1340                state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
     1341                setspc=false;
     1342                switch (state) {
     1343                        case 0:
     1344                                // cold, dry base: apply basal surface forcing
     1345                                break;
     1346                        case 1:
     1347                                // cold, wet base: keep at pressure melting point
     1348                                setspc=true;
     1349                                break;
     1350                        case 2:
     1351                                // temperate, thin refreezing base: release spc
     1352                                break;
     1353                        case 3:
     1354                                // temperate, thin melting base: set spc
     1355                                setspc=true;
     1356                                break;
     1357                        case 4:
     1358                                // temperate, thick melting base: s
     1359                                setspc=true;
     1360                                break;
     1361                        default:
     1362                                _printf0_("     unknown thermal basal state found!");
     1363                }
     1364
     1365                /*apply or release spc*/
    12311366                node=element->GetNode(indices[i]);
    1232                 if(setspc)
    1233                         node->ApplyConstraint(0,h_pmp); /*apply spc*/
     1367                if(setspc){
     1368                        pressure_input->GetInputValue(&pressure, gauss);
     1369                        node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
     1370                }
    12341371                else                   
    1235                         node->DofInFSet(0); /*remove spc*/
     1372                        node->DofInFSet(0);
    12361373        }
    12371374
    … …  
    12411378        delete gauss;
    12421379        delete gaussup;
     1380}/*}}}*/
     1381int EnthalpyAnalysis::GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpyup, IssmDouble pressure, IssmDouble pressureup, IssmDouble watercolumn, IssmDouble meltingrate){/*{{{*/
     1382
     1383        /* Check if ice in element */
     1384        if(!element->IsIceInElement()) return -1;
     1385
     1386        /* Only update Constraints at the base of grounded ice*/
     1387        if(!(element->IsOnBase())) return -1;
     1388
     1389        /*Intermediary*/
     1390        int state=-1;
     1391        IssmDouble      dt;
     1392
     1393        /*Get parameters and inputs: */
     1394        element->FindParam(&dt,TimesteppingTimeStepEnum);
     1395
     1396        if(dt==0.){ // steadystate case
     1397                state=0; //TODO: add consistent steadystate basal condition scheme
     1398//              if(enthalpy<PureIceEnthalpy(element,pressure)){ /*is base cold?*/
     1399//                      if(watercolumn<=0.) state=0; // cold, dry base
     1400//                      else state=1; // cold, wet base (refreezing)
     1401//              }
     1402//              else{ /*base is temperate, check if upper node is temperate, too.*/
     1403//                      if(enthalpyup<PureIceEnthalpy(element,pressureup))
     1404//                              if(meltingrate<0.)      state=2; // refreezing temperate base (non-physical, only for steadystate solver)
     1405//                              else                                            state=3; // melting temperate base with no temperate layer
     1406//                      else
     1407//                              state=4; // melting temperate base with temperate layer of positive thickness
     1408//              }
     1409        }
     1410        else{ // transient case
     1411                if(enthalpy<PureIceEnthalpy(element,pressure)){
     1412                        if(watercolumn<=0.) state=0; // cold, dry base
     1413                        else state=1; // cold, wet base (refreezing)
     1414                }
     1415                else{
     1416                        if(enthalpyup<PureIceEnthalpy(element,pressureup))      state=3; // temperate base, but no temperate layer
     1417                        else state=4; // temperate layer with positive thickness
     1418                }
     1419        }
     1420
     1421        _assert_(state>=0);
     1422        return state;
     1423}/*}}}*/
     1424IssmDouble EnthalpyAnalysis::GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure){/*{{{*/
     1425
     1426        IssmDouble temperature, waterfraction;
     1427        IssmDouble kappa_w = 0.6; // thermal conductivity of water (in W/m/K)
     1428        IssmDouble kappa_i = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
     1429        element->EnthalpyToThermal(&temperature, &waterfraction, enthalpy, pressure);
     1430
     1431        return (1.-waterfraction)*kappa_i + waterfraction*kappa_w;
    12431432}/*}}}*/
    12441433

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

    r18057 r18612  
    88/*Headers*/
    99#include "./Analysis.h"
     10#include "../classes/classes.h"
    1011
    1112class EnthalpyAnalysis: public Analysis{
    … …  
    4142                /*Modules*/
    4243                static void PostProcessing(FemModel* femmodel);
     44                static void ComputeBasalMeltingrate(FemModel* femmodel);
    4345                static void ComputeBasalMeltingrate(Element* element);
     46                static void DrainWaterfraction(FemModel* femmodel);
    4447                static void DrainWaterfractionIcecolumn(Element* element);
    4548                static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
     49                static void UpdateBasalConstraints(FemModel* femmodel);
    4650                static void UpdateBasalConstraints(Element* element);
     51                static void UpdateBasalConstraintsTransient(Element* element);
     52                static void UpdateBasalConstraintsSteadystate(Element* element);
     53//              static int GetThermalBasalCondition(Element* element, Gauss* gauss, Gauss* gaussup, int enthalpy_enum);
     54//              static int GetThermalBasalCondition(Element* element, GaussPenta* gauss, GaussPenta* gaussup, int enthalpy_enum);
     55                static int GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
     56                static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);
     57
    4758
    4859                /*Intermediaries*/
Note: See TracChangeset for help on using the changeset viewer.