Changeset 16026


Ignore:
Timestamp:
08/30/13 06:48:44 (12 years ago)
Author:
jbondzio
Message:

ADD: update enthalpy module

Location:
issm/trunk-jpl/src
Files:
14 edited

Legend:

Unmodified
Added
Removed

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

    r16007 r16026  
    306306                                        ./modules/PositiveDegreeDayx/PositiveDegreeDayx.h\
    307307                                        ./modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp\
     308                                        ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\
     309                                        ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\
    308310                                        ./modules/Delta18oParameterizationx/Delta18oParameterizationx.h\
    309311                                        ./modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp\

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

    r15849 r16026  
    9292                if(isenthalpy)  InputToResultx(femmodel,WaterfractionEnum);
    9393                if(isenthalpy)  InputToResultx(femmodel,EnthalpyEnum);
    94                 if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
     94                if(isenthalpy)  InputToResultx(femmodel,WatercolumnEnum);
     95                //if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
     96                InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
    9597                femmodel->RequestedOutputsx(requested_outputs,numoutputs);
    9698        }

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

    r15983 r16026  
    111111                        else{
    112112                                enthalpy_core(femmodel);
     113                                PostprocessingEnthalpyx(femmodel);
    113114                        }
    114115                        #else

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

    r16014 r16026  
    127127                #ifdef _HAVE_THERMAL_
    128128                virtual void UpdateThermalBasalConstraints(void)=0;
    129                 //virtual void ComputeBasalMeltingrate(void)=0;
    130                 //virtual void DrainWaterfraction(void)=0;
     129                virtual void ComputeBasalMeltingrate(void)=0;
     130        virtual void DrainWaterfraction(void)=0;
    131131                #endif
    132132

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

    r16023 r16026  
    47004700        IssmDouble scalar,enthalpy,enthalpyup;
    47014701        IssmDouble pressure,pressureup;
     4702    IssmDouble watercolumn;
    47024703        IssmDouble basis[NUMVERTICES];
    47034704        Friction*  friction=NULL;
    … …  
    47244725        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    47254726        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    4726 
    4727         /*Build frictoin element, needed later: */
     4727    Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
     4728
     4729        /*Build friction element, needed later: */
    47284730        friction=new Friction("3d",inputs,matpar,analysis_type);
    47294731
    … …  
    47414743                enthalpy_input->GetInputValue(&enthalpy,gauss);
    47424744                pressure_input->GetInputValue(&pressure,gauss);
    4743 //              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
    4744 //                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    4745 //                      pressure_input->GetInputValue(&pressureup,gaussup);
    4746 //                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
    4747 //                              //do nothing, don't add heatflux
    4748 //                      }
    4749 //                      else{
    4750 //                              //need to change spcenthalpy according to Aschwanden
    4751 //                              //NEED TO UPDATE
    4752 //                      }
    4753 //              }
    4754 //              else{
    4755                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4756                         friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4757                         vx_input->GetInputValue(&vx,gauss);
    4758                         vy_input->GetInputValue(&vy,gauss);
    4759                         basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4760 
    4761                         scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
    4762                         if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    4763 
    4764                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    4765 //              }
     4745                if(enthalpy>=matpar->PureIceEnthalpy(pressure)){
     4746                        enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     4747                        pressure_input->GetInputValue(&pressureup,gaussup);
     4748            if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
     4749                                // temperate ice has positive thickness: grad enthalpy*n=0.
     4750                        }
     4751                        else{
     4752                                // only base temperate, set Dirichlet BCs in Penta::UpdateThermalBasalConstraints()
     4753                        }
     4754                }
     4755                else{
     4756            watercolumn_input->GetInputValue(&watercolumn,gauss);
     4757            if(watercolumn==0.){
     4758                /*add geothermal heat flux and basal friction*/
     4759                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4760                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4761                vx_input->GetInputValue(&vx,gauss);
     4762                vy_input->GetInputValue(&vy,gauss);
     4763                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4764
     4765                scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     4766                if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
     4767
     4768                for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4769            }
     4770            else { /*do nothing (water layer acts as insulation)*/ }
     4771                }
    47664772        }
    47674773
    … …  
    52025208                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    52035209                                break;
     5210            case LliboutryDuvalEnum:
     5211                B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,                                                   (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,                                       material->GetN());
     5212                for(i=0;i<numdof;i++) B[i]=B_average;
     5213                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
     5214                                break;
    52045215                        default:
    52055216                                _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
    … …  
    52205231        /*Intermediary*/
    52215232        bool        isenthalpy,isdynamicbasalspc,istemperatelayer;
    5222         int         numindices;
    5223         IssmDouble  h_pmp,pressure;
    5224         int        *indices = NULL;
     5233        int         numindices, numindicesup;
     5234        IssmDouble  h_pmp,pressure, pressureup;
     5235    IssmDouble  enthalpy, enthalpyup;
     5236        int        *indices = NULL, *indicesup = NULL;
    52255237
    52265238        /* Only update Constraints at the base of grounded ice*/
    5227         if(!IsOnBed() || !IsFloating()) return;
    5228 
    5229         /*Check wether dynamic basal boudary conditions are activated -> TODO: Johannes :) */
     5239        if(!IsOnBed() || IsFloating()) return;
     5240
     5241        /*Check wether dynamic basal boundary conditions are activated */
    52305242        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    52315243        if(!isenthalpy) return;
    … …  
    52345246        if(!isdynamicbasalspc) return;
    52355247
    5236 
    5237         /*Fetch indices of basal nodes for this finite element*/
     5248        /*Fetch indices of basal & surface nodes for this finite element*/
    52385249        BasalNodeIndices(&numindices,&indices,this->element_type);
     5250    SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
     5251    _assert_(numindices==numindicesup);
    52395252
    52405253        /*Get parameters and inputs: */
    52415254        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    5242 
    5243         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     5255    Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     5256
     5257        /*if there is a temperate layer of positive thickness, set enthalpy=h_pmp at that node*/
    52445258        GaussPenta* gauss=new GaussPenta();
     5259    GaussPenta* gaussup=new GaussPenta();
    52455260        for(int i=0;i<numindices;i++){
    52465261                gauss->GaussNode(this->element_type,indices[i]);
     5262        gaussup->GaussNode(this->element_type,indicesup[i]); // TODO: check: are the nodes corresponding?
    52475263
    52485264                /*Check wether there is a temperate layer at the base or not -> TODO: Johannes:) */
    5249                 istemperatelayer = false;
    5250 
    5251                 /*Add Dirichlet constraint to this node if there is a positive thickness of temperate ice*/
    5252                 if(istemperatelayer){
    5253 
     5265        /*check if node is temperate, if not, return*/
     5266        enthalpy_input->GetInputValue(&enthalpy, gauss);
     5267                pressure_input->GetInputValue(&pressure, gauss);
     5268        if (enthalpy<matpar->PureIceEnthalpy(pressure)){
     5269          // TODO: reset, if necessary, all spcs to non-valid
     5270          continue;
     5271        }
     5272        /*check if upper node is temperate. if yes, then we have a temperate layer of positive thickness. if not, continue.*/
     5273        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
     5274                pressure_input->GetInputValue(&pressureup, gaussup);   
     5275        istemperatelayer = false;
     5276        if (enthalpyup>=matpar->PureIceEnthalpy(pressureup))
     5277          istemperatelayer=true;
     5278
     5279                /*Add Dirichlet constraint to this node if there is no layer of temperate ice with positive thickness*/
     5280                if(!istemperatelayer){
    52545281                        /*Calculate enthalpy at pressure melting point */
    5255                         pressure_input->GetInputValue(&pressure,gauss);
    52565282                        h_pmp=matpar->PureIceEnthalpy(pressure);
    52575283
    5258 
    52595284                        /*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*/
    5260                         nodes[indices[i]]->ApplyConstraint(0,h_pmp);
    5261                 }
     5285                        nodes[indices[i]]->ApplyConstraint(1,h_pmp);
     5286                }
     5287        else {
     5288          /*remove spc*/
     5289          nodes[indices[i]]->DofInFSet(0);
     5290        }
    52625291        }
    52635292
    52645293        /*Free ressources:*/
    52655294        xDelete<int>(indices);
    5266 }
    5267 /*}}}*/
    5268 /*FUNCTION Penta::ComputeBasalMeltrate{{{*/
    5269 void Penta::ComputeBasalMeltrate(void){
     5295    xDelete<int>(indicesup);
     5296}
     5297/*}}}*/
     5298/*FUNCTION Penta::ComputeBasalMeltingrate{{{*/
     5299void Penta::ComputeBasalMeltingrate(void){
    52705300        /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
    52715301
    52725302        /* Intermediaries */
    5273         bool        isenthalpy;
     5303        bool        isenthalpy, checkpositivethickness, istemperatelayer;
    52745304        int         i,j,analysis_type;
    52755305        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    52795309        IssmDouble  normal_base[3], d1enthalpy[3];
    52805310        IssmDouble  kappa;
    5281     IssmDouble  meltrate[3], watercolumn[3];
    5282         IssmDouble  enthalpy, enthalpyup;
     5311    IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy;
     5312        IssmDouble  enthalpy[NUMVERTICES], enthalpyup;
    52835313        IssmDouble  pressure, pressureup;
    52845314        IssmDouble  temperature, waterfraction;
    5285         IssmDouble  latentheat;
     5315        IssmDouble  latentheat, rho_ice;
    52865316        IssmDouble  basalfriction, alpha2;
    52875317        IssmDouble  vx,vy,vz;
     5318    IssmDouble  connectivity;
    52885319        IssmDouble  dt;
    52895320        Friction   *friction  = NULL;
    … …  
    52985329        /*Fetch parameters and inputs */
    52995330        latentheat=matpar->GetLatentHeat();
     5331    rho_ice=this->matpar->GetRhoIce();
    53005332        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
     5333    Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(basalmeltingrate_input);
    53015334        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
    53025335        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    … …  
    53125345        friction=new Friction("3d",inputs,matpar,analysis_type);
    53135346
    5314         /*Ok, get meltrates now from basal conditions*/
     5347        /*Ok, get meltingrates now from basal conditions*/
    53155348        GaussPenta* gauss=new GaussPenta();
    53165349        GaussPenta* gaussup=new GaussPenta();
    … …  
    53185351                gauss->GaussVertex(iv);
    53195352                gaussup->GaussVertex(iv+3);
    5320 
    5321         // TODO: make sure that no node is computed twice (insert mask)
    5322 
     5353               
     5354        checkpositivethickness=true;
    53235355                watercolumn_input->GetInputValue(&watercolumn[iv], gauss);
    5324                 enthalpy_input->GetInputValue(&enthalpy, gauss);
     5356        basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss);
     5357                enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
    53255358                pressure_input->GetInputValue(&pressure, gauss);
    53265359
    5327                 /*Calculate basal meltrate*/
    5328                 if((watercolumn[iv]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
    5329                         enthalpy=matpar->PureIceEnthalpy(pressure);
    5330                 }
    5331                 else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
    5332                         meltrate[iv]=0.;   
    5333                         watercolumn[iv]=0.;
    5334                         continue;
    5335                 }
    5336 
    5337                 /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice) */
    5338                 enthalpy_input->GetInputValue(&enthalpyup, gaussup);
    5339                 pressure_input->GetInputValue(&pressureup, gaussup);   
    5340                 if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
    5341             for(i=0;i<3;i++) vec_heatflux[i]=0.;
    5342                 }
    5343                 else{
    5344                         enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    5345                         kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    5346                         for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    5347                 }
    5348 
    5349                 /*Get normal vector to the bed */
    5350                 BedNormal(&normal_base[0],xyz_list_tria);
    5351 
    5352                 heatflux=0.;
    5353                 for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
    5354                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    5355 
    5356                 matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure);
    5357 
    5358                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    5359                 vx_input->GetInputValue(&vx,gauss);
    5360                 vy_input->GetInputValue(&vy,gauss);
    5361                 vz_input->GetInputValue(&vz,gauss);
    5362                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
    5363                 meltrate[iv]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat;
    5364 
    5365                 /*Update water column*/
     5360                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
     5361        meltingrate_enthalpy=0.;
     5362                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure))){
     5363            /*ensure that no ice is at T<Tm(p), if water layer present*/
     5364            enthalpy[iv]=matpar->PureIceEnthalpy(pressure);
     5365            //meltingrate_enthalpy=0.;
     5366                }
     5367                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure)){
     5368            /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
     5369                        meltingrate_enthalpy=0.;
     5370            checkpositivethickness=false;
     5371                }
     5372        else {/*do nothing, go to next check*/}
     5373
     5374        if(checkpositivethickness){
     5375            /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
     5376            istemperatelayer=false;
     5377            enthalpy_input->GetInputValue(&enthalpyup, gaussup);
     5378            pressure_input->GetInputValue(&pressureup, gaussup);   
     5379            if(enthalpyup>=matpar->PureIceEnthalpy(pressureup))
     5380                istemperatelayer=true;
     5381            if(istemperatelayer)
     5382                for(i=0;i<3;i++) vec_heatflux[i]=0.;
     5383            else{
     5384                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
     5385                kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure);
     5386                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
     5387            }
     5388
     5389            /*geothermal heatflux*/
     5390            BedNormal(&normal_base[0],xyz_list_tria);
     5391            heatflux=0.;
     5392            for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
     5393            geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     5394
     5395            /*basal friction*/
     5396            friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     5397            vx_input->GetInputValue(&vx,gauss);
     5398            vy_input->GetInputValue(&vy,gauss);
     5399            vz_input->GetInputValue(&vz,gauss);
     5400            basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
     5401
     5402            matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure);
     5403            meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent
     5404        }
     5405
     5406                /*Update water column, basal meltingrate*/
     5407        connectivity=IssmDouble(vertices[iv]->Connectivity());
     5408        meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex
     5409        basalmeltingrate[iv]+=meltingrate_enthalpy;
    53665410                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    53675411                if(reCast<bool,IssmDouble>(dt))
    5368             watercolumn[iv]+=dt*meltrate[iv];
     5412            watercolumn[iv]+=dt*meltingrate_enthalpy;
    53695413                else
    5370             watercolumn[iv]=meltrate[iv];
     5414            watercolumn[iv]+=meltingrate_enthalpy;
    53715415        } 
    5372     /*update meltrate and watercolumn*/
     5416    /*feed updated variables back into model*/
     5417    this->inputs->AddInput(new PentaInput(EnthalpyEnum, enthalpy, P1Enum));
    53735418    this->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumn, P1Enum));
    5374     this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum, meltrate, P1Enum));
     5419    this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum, basalmeltingrate, P1Enum));
    53755420
    53765421        /*Clean up and return*/
    53775422        delete gauss;
    53785423        delete gaussup;
     5424    delete friction;
    53795425}
    53805426/*}}}*/
    … …  
    53825428void Penta::DrainWaterfraction(void){
    53835429   
    5384 // TODO: create vector to mark whether node has been drained already i.o. to not drain nodes multiple times
    5385 
    53865430    /*Intermediaries*/
    53875431    bool isenthalpy;
    … …  
    53895433    IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    53905434    IssmDouble latentheat, dt;
     5435    IssmDouble connectivity;
    53915436    GaussPenta* gauss;
    53925437   
    … …  
    54015446    this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    54025447    latentheat=matpar->GetLatentHeat();
     5448
    54035449    gauss=new GaussPenta();
    54045450    for(int iv=0;iv<NUMVERTICES;iv++){
    54055451        gauss->GaussVertex(iv);
    5406         /*TODO: Check whether Vertex has been drained already*/
    54075452        enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
    54085453        pressure_input->GetInputValue(&pressure[iv], gauss);
    54095454        matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
    54105455   
    5411         /*drain water fraction*/
    5412         waterfraction[iv]-=dt*DrainageFunctionWaterfraction(waterfraction[iv]);
    5413         if(waterfraction[iv]<0.) waterfraction[iv]=0.;
    5414         /*update enthalpy*/
     5456        /*drain water fraction & update enthalpy*/
     5457        // TODO: replace connectivity-wise draining by draining once per node per timestep
     5458        connectivity=IssmDouble(vertices[iv]->Connectivity());
     5459        waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity;
    54155460        matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
    54165461    }
    … …  
    54185463    this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    54195464    this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    5420     this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));   
     5465    // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    // temperature should not change during drainage...
    54215466}
    54225467/*}}}*/

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

    r16023 r16026  
    357357                void           InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong);
    358358                void           UpdateThermalBasalConstraints(void);
    359                 void           ComputeBasalMeltrate(void);
     359                void           ComputeBasalMeltingrate(void);
    360360                void           DrainWaterfraction(void);
    361361                #endif

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

    r16023 r16026  
    255255                #ifdef _HAVE_THERMAL_
    256256                void UpdateThermalBasalConstraints(void){_error_("not implemented yet");};
    257                 void ComputeBasalMeltrate(void){_error_("not implemented yet");};
     257                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
     258                void DrainWaterfraction(void){_error_("not implemented yet");};
    258259                #endif
    259260

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp

    r15986 r16026  
    4444        iomodel->FetchDataToInput(elements,WaterfractionEnum);
    4545        iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum);
     46        iomodel->FetchDataToInput(elements,WatercolumnEnum);
     47        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
    4648        iomodel->FetchDataToInput(elements,VxEnum);
    4749        iomodel->FetchDataToInput(elements,VyEnum);

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

    r15082 r16026  
    7575#include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
    7676#include "./PositiveDegreeDayx/PositiveDegreeDayx.h"
     77#include "./PostprocessingEnthalpyx/PostprocessingEnthalpyx.h"
    7778#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
    7879#include "./Reduceloadx/Reduceloadx.h"

  • issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp

    r15891 r16026  
    88
    99/*FUNCTION IssmDouble DrainageFunctionWaterfraction()*/
    10 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction){
     10IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.){
    1111    /* DrainageFunctionWaterfraction returns how much of the waterfraction is drained per year */
    12     _assert_(waterfraction>0);
     12    _assert_(waterfraction>=0.);
     13    _assert_(dt>=0.);
    1314
    1415    IssmDouble w0=0.01, w1=0.02, w2=0.03;
    15     IssmDouble D0=0, D1=0.005, D2=0.05;
     16    IssmDouble Dret, D0=0, D1=0.005, D2=0.05;
    1617
     18    /*get drainage function value*/
    1719    if((w0==w1)||(w1==w2)||(w0==w2))
    18         perror ("error in DrainageFunctionWaterfraction. Abort");
     20        _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
    1921    if(waterfraction<=w0)
    20         return D0;
     22        Dret=D0;
    2123    if((waterfraction>w0) && (waterfraction<=w1))
    22         return (D1-D0)/(w1-w0)*(waterfraction-w0)+D0;
     24        Dret=(D1-D0)/(w1-w0)*(waterfraction-w0)+D0;
    2325    if((waterfraction>w1) && (waterfraction<=w2))
    24         return (D2-D1)/(w2-w1)*(waterfraction-w1)+D1;
     26        Dret=(D2-D1)/(w2-w1)*(waterfraction-w1)+D1;
    2527    else
    26         return D2;
     28        Dret=D2;
     29   
     30    /*check if dt*Dret>waterfraction. If so, drain whole waterfraction*/
     31    if(dt==0.){
     32      if(Dret>waterfraction)
     33        return waterfraction;
     34      else
     35        return Dret;
     36    }
     37    else{
     38      if(dt*Dret>waterfraction)
     39        return waterfraction;
     40      else
     41        return dt*Dret;
     42    }
    2743}

  • issm/trunk-jpl/src/c/shared/Elements/elements.h

    r15897 r16026  
    2020                                     IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
    2121                                          IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
    22 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction);
     22IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
    2323
    2424/*Print arrays*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r16010 r16026  
    236236        SurfaceforcingsBPosEnum,
    237237        SurfaceforcingsBNegEnum,
     238        ThermalIsenthalpyEnum,
     239        ThermalIsdynamicbasalspcEnum,
    238240        ThermalMaxiterEnum,
    239241        ThermalPenaltyFactorEnum,
    … …  
    242244        ThermalSpctemperatureEnum,
    243245        ThermalStabilizationEnum,
    244         ThermalIsenthalpyEnum,
    245246        GiaMantleViscosityEnum,
    246247        GiaLithosphereThicknessEnum,

  • issm/trunk-jpl/src/m/classes/model/model.m

    r15987 r16026  
    671671                        if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end;
    672672                        if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end;
     673            if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);end;
    673674
    674675                        %bedinfo and surface info

  • issm/trunk-jpl/src/m/classes/thermal.m

    r15133 r16026  
    1313                penalty_factor    = 0;
    1414                isenthalpy        = 0;
     15                isdynamicbasalspc = 0;
    1516        end
    1617        methods
    … …  
    3940                        %Should we use cold ice (default) or enthalpy formulation
    4041                        obj.isenthalpy=0;
     42           
     43                        %will basal boundary conditions be set dynamically
     44                        obj.isdynamicbasalspc=0;
    4145                end % }}}
    4246                function md = checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5256                                md = checkfield(md,'thermal.spctemperature(find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN))','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos),'message','spctemperature should be below the adjusted melting point');
    5357                                md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0 1]);
     58                                md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0 1]);
    5459                        end
    5560                end % }}}
    … …  
    6368                        fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
    6469                        fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    65 
     70                        fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
     71           
    6672                end % }}}
    6773                function marshall(obj,md,fid) % {{{
    … …  
    7379                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
    7480                        WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
    75                 end % }}}
     81                        WriteData(fid,'object',obj,'fieldname','isdynamicbasalspc','format','Boolean');
     82                end % }}}
    7683        end
    7784end
Note: See TracChangeset for help on using the changeset viewer.