Changeset 16359


Ignore:
Timestamp:
10/10/13 06:55:53 (11 years ago)
Author:
jbondzio
Message:

BUG: Enthalpy: basal melting and refreezing works now for zero velocities

Location:
issm/trunk-jpl/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r16219 r16359  
    3030                InputToResultx(femmodel,EnthalpyEnum);
    3131                InputToResultx(femmodel,WaterfractionEnum);
     32                InputToResultx(femmodel,WatercolumnEnum);
     33                InputToResultx(femmodel,BasalforcingsMeltingRateEnum);
    3234        }
    3335}

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

    r16357 r16359  
    45264526        IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    45274527        IssmDouble Jdet2d,dt;
    4528         IssmDouble rho_ice,heatcapacity,geothermalflux_value;
     4528        IssmDouble rho_ice,heatcapacity;
     4529        IssmDouble geothermalflux_value, heatflux_value;
    45294530        IssmDouble basalfriction,alpha2,vx,vy,vz;
    45304531        IssmDouble scalar,enthalpy,enthalpyup;
    … …  
    45634564        gauss=new GaussPenta(0,1,2,2);
    45644565        gaussup=new GaussPenta(3,4,5,2);
     4566
    45654567        for(int ig=gauss->begin();ig<gauss->end();ig++){
    45664568
    … …  
    45754577                watercolumn_input->GetInputValue(&watercolumn,gauss);
    45764578
    4577                 if((watercolumn==0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
     4579                if((watercolumn<=0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
    45784580                        /* the above check is equivalent to
    45794581                         NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
    … …  
    45864588                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
    45874589
    4588                         scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     4590                        heatflux_value=(basalfriction+geothermalflux_value)/(rho_ice);
     4591                        scalar=gauss->weight*Jdet2d*heatflux_value;
    45894592                        if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    45904593                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    … …  
    50145017        bool        isenthalpy,isdynamicbasalspc,setspc;
    50155018        int         numindices, numindicesup;
    5016         IssmDouble  h_pmp,pressure, pressureup;
    5017         IssmDouble  enthalpy, enthalpyup;
     5019        IssmDouble  pressure, pressureup;
     5020        IssmDouble  h_pmp, enthalpy, enthalpyup;
     5021        IssmDouble  watercolumn;
    50185022        int        *indices = NULL, *indicesup = NULL;
    50195023
    … …  
    50355039        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    50365040        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     5041        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); //_assert_(watercolumn_input);
    50375042
    50385043        /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
    … …  
    50475052                enthalpy_input->GetInputValue(&enthalpy, gauss);
    50485053                pressure_input->GetInputValue(&pressure, gauss);
     5054                watercolumn_input->GetInputValue(&watercolumn,gauss);
    50495055                setspc = false;
    50505056                if (enthalpy>=matpar->PureIceEnthalpy(pressure)){               
    … …  
    50545060                        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
    50555061                        pressure_input->GetInputValue(&pressureup, gaussup);   
    5056                         setspc=(enthalpyup<matpar->PureIceEnthalpy(pressureup))?true:false;
     5062                        setspc=((enthalpyup<matpar->PureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false;
    50575063                }
    50585064
    … …  
    50965102        IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
    50975103        IssmDouble  geothermalflux[NUMVERTICES];
    5098         IssmDouble  dt,meltingrate_enthalpy;
     5104        IssmDouble  dt;
     5105        IssmDouble  meltingrate_enthalpy;
    50995106        Friction   *friction  = NULL;
    51005107
    … …  
    51085115        /*Fetch parameters and inputs */
    51095116        latentheat=matpar->GetLatentHeat();
    5110         rho_ice=this->matpar->GetRhoIce();
     5117        rho_ice=matpar->GetRhoIce();
    51115118        GetInputListOnVertices(&vx[0],VxEnum);
    51125119        GetInputListOnVertices(&vy[0],VyEnum);
    … …  
    51495156                        istemperatelayer=false;
    51505157                        if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true;
    5151                         if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.;
     5158                        if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
    51525159                        else{
    51535160                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    5154                                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]); _assert_(kappa>0.);
     5161                                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy[iv],enthalpy[iv+NUMVERTICES2D], pressure[iv],pressure[iv+NUMVERTICES2D]); _assert_(kappa>0.);
    51555162                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    51565163                        }
    … …  
    51665173
    51675174                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
    5168                         meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux[iv]))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent
    5169                 }
    5170 
     5175                        // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
     5176                        meltingrate_enthalpy=(heatflux+basalfriction+geothermalflux[iv])/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent
     5177                }
    51715178                /*Update water column, basal meltingrate*/
    5172                 basalmeltingrate[iv]+=meltingrate_enthalpy;
     5179                //basalmeltingrate[iv]+=meltingrate_enthalpy;
     5180                basalmeltingrate[iv]=meltingrate_enthalpy;
    51735181                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    51745182                if(reCast<bool,IssmDouble>(dt))
    5175                  watercolumn[iv]+=dt*meltingrate_enthalpy;
     5183                        watercolumn[iv]+=dt*meltingrate_enthalpy;
    51765184                else
    5177                  watercolumn[iv]+=meltingrate_enthalpy;
     5185                        watercolumn[iv]+=meltingrate_enthalpy;
    51785186        } 
    51795187        /*feed updated variables back into model*/
    … …  
    51935201        bool isenthalpy;
    51945202        IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
     5203        IssmDouble watercolumnbase[NUMVERTICES];
    51955204        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    51965205        IssmDouble latentheat, dt;
     5206        IssmDouble dw, dwc;
     5207        Penta *pentabase = NULL;
    51975208
    51985209        /*Check wether enthalpy is activated*/
    51995210        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    52005211        if(!isenthalpy) return;       
    5201 
     5212       
     5213        /*get basal element, needed for basal watercolumn*/
     5214        pentabase=this->GetBasalElement();
     5215       
    52025216        GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
    52035217        GetInputListOnVertices(&pressure[0],PressureEnum);
     5218        pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum);
     5219
    52045220        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    52055221        latentheat=matpar->GetLatentHeat();
    … …  
    52075223        for(int iv=0;iv<NUMVERTICES;iv++){
    52085224                matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
    5209 
     5225                dw=DrainageFunctionWaterfraction(waterfraction[iv], dt);
    52105226                /*drain water fraction & update enthalpy*/
    5211                 waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt);
     5227                waterfraction[iv]-=dw;
    52125228                matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
     5229                /*add drained water to watercolumn at base*/
     5230                dwc=dw*this->IceVolume();
     5231                watercolumnbase[iv%NUMVERTICES2D]+=dwc;
    52135232        }
    52145233        /*feed updated results back into model*/
    52155234        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    52165235        this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    5217         // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    // temperature should not change during drainage...
     5236        pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum));
     5237
     5238        delete pentabase;
    52185239}
    52195240/*}}}*/

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r16322 r16359  
    385385}
    386386/*}}}*/
     387/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
     388IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1){
     389        /*returns kappa depending on distribution of enthalpy over edge of element
     390                lambda is the barycentric coordinate that solves H0+(H1-H0)*lambda=H_pureice.
     391                it represents the fraction of the ice column which is temperate/cold like H0.
     392                if lambda<0 or lambda>1, then the whole ice column is cold or temperate, respectively.
     393        */
     394        IssmDouble kappa, kappa0, kappa1;
     395        if (enthalpy0!=enthalpy1){
     396                IssmDouble lambda=(PureIceEnthalpy(pressure0)-enthalpy0)/(enthalpy1-enthalpy0);
     397                if ((lambda>=0.) && (lambda<=1.)){
     398                        kappa0=GetEnthalpyDiffusionParameter(enthalpy0,pressure0);
     399                        kappa1=GetEnthalpyDiffusionParameter(enthalpy1,pressure1);
     400                        kappa=lambda*kappa0+(1.-lambda)*kappa1;
     401                }
     402        }
     403        else
     404                kappa=GetEnthalpyDiffusionParameter(enthalpy0, pressure0);
     405        return kappa;
     406}
     407/*}}}*/
    387408/*FUNCTION Matpar::EnthalpyToThermal {{{*/
    388409void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r16350 r16359  
    126126                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    127127                IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
     128                IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1);
    128129                IssmDouble GetLithosphereShearModulus();
    129130                IssmDouble GetLithosphereDensity();

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

    r16144 r16359  
    1515    IssmDouble w0=0.01, w1=0.02, w2=0.03;
    1616    IssmDouble Dret, D0=0, D1=0.005, D2=0.05;
     17                IssmDouble yts=365*24*60*60;
     18                dt/=yts;
    1719
    1820    /*get drainage function value*/
Note: See TracChangeset for help on using the changeset viewer.