Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16030) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16031) @@ -4700,5 +4700,5 @@ IssmDouble scalar,enthalpy,enthalpyup; IssmDouble pressure,pressureup; - IssmDouble watercolumn; + IssmDouble watercolumn; IssmDouble basis[NUMVERTICES]; Friction* friction=NULL; @@ -4725,5 +4725,5 @@ Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); /*Build friction element, needed later: */ @@ -4746,5 +4746,5 @@ enthalpy_input->GetInputValue(&enthalpyup,gaussup); pressure_input->GetInputValue(&pressureup,gaussup); - if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ + if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ // temperate ice has positive thickness: grad enthalpy*n=0. } @@ -4754,19 +4754,19 @@ } else{ - watercolumn_input->GetInputValue(&watercolumn,gauss); - if(watercolumn==0.){ - /*add geothermal heat flux and basal friction*/ - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); - if(reCast(dt)) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; - } - else { /*do nothing (water layer acts as insulation)*/ } + watercolumn_input->GetInputValue(&watercolumn,gauss); + if(watercolumn==0.){ + /*add geothermal heat flux and basal friction*/ + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + if(reCast(dt)) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*basis[i]; + } + else /*do nothing (water layer acts as insulation)*/ } } @@ -5148,6 +5148,6 @@ const int numdof=NDOF1*NUMVERTICES; - bool converged=false; - int i,rheology_law; + bool converged=false; + int i,rheology_law; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble values[numdof]; @@ -5157,5 +5157,5 @@ IssmDouble B[numdof]; IssmDouble B_average,s_average; - int* doflist=NULL; + int* doflist=NULL; /*Get dof list: */ @@ -5208,12 +5208,13 @@ this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; - case LliboutryDuvalEnum: - 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()); - for(i=0;iGetN()); + for(i=0;imaterial->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet"); - } } @@ -5305,9 +5306,8 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble heatflux, geothermalflux_value; + IssmDouble heatflux, geothermalflux_value,kappa; IssmDouble vec_heatflux[3]; IssmDouble normal_base[3], d1enthalpy[3]; - IssmDouble kappa; - IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy; + IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy; IssmDouble enthalpy[NUMVERTICES], enthalpyup; IssmDouble pressure, pressureup; @@ -5315,7 +5315,5 @@ IssmDouble latentheat, rho_ice; IssmDouble basalfriction, alpha2; - IssmDouble vx,vy,vz; - IssmDouble connectivity; - IssmDouble dt; + IssmDouble vx,vy,vz,dt,connectivity; Friction *friction = NULL; @@ -5329,12 +5327,12 @@ /*Fetch parameters and inputs */ latentheat=matpar->GetLatentHeat(); - rho_ice=this->matpar->GetRhoIce(); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basalmeltingrate_input); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_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); + rho_ice=this->matpar->GetRhoIce(); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basalmeltingrate_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_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); Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); @@ -5351,76 +5349,74 @@ gauss->GaussVertex(iv); gaussup->GaussVertex(iv+3); - - checkpositivethickness=true; + + checkpositivethickness=true; watercolumn_input->GetInputValue(&watercolumn[iv], gauss); - basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss); + basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss); enthalpy_input->GetInputValue(&enthalpy[iv], gauss); pressure_input->GetInputValue(&pressure, gauss); /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/ - meltingrate_enthalpy=0.; + meltingrate_enthalpy=0.; if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure))){ - /*ensure that no ice is at TPureIceEnthalpy(pressure); - //meltingrate_enthalpy=0.; + /*ensure that no ice is at TPureIceEnthalpy(pressure); + //meltingrate_enthalpy=0.; } else if(enthalpy[iv]PureIceEnthalpy(pressure)){ - /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/ + /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/ meltingrate_enthalpy=0.; - checkpositivethickness=false; - } - else {/*do nothing, go to next check*/} - - if(checkpositivethickness){ - /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */ - istemperatelayer=false; - enthalpy_input->GetInputValue(&enthalpyup, gaussup); - pressure_input->GetInputValue(&pressureup, gaussup); - if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) - istemperatelayer=true; - if(istemperatelayer) - for(i=0;i<3;i++) vec_heatflux[i]=0.; - else{ - enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure); - for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; - } - - /*geothermal heatflux*/ - BedNormal(&normal_base[0],xyz_list_tria); - heatflux=0.; - for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - - /*basal friction*/ - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); - - matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure); - meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent - } + checkpositivethickness=false; + } + else {/*do nothing, go to next check*/} + + if(checkpositivethickness){ + /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */ + istemperatelayer=false; + enthalpy_input->GetInputValue(&enthalpyup, gaussup); + pressure_input->GetInputValue(&pressureup, gaussup); + if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) istemperatelayer=true; + if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; + else{ + enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure); + for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; + } + + /*geothermal heatflux*/ + BedNormal(&normal_base[0],xyz_list_tria); + heatflux=0.; + for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + + /*basal friction*/ + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); + + matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure); + meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent + } /*Update water column, basal meltingrate*/ - connectivity=IssmDouble(vertices[iv]->Connectivity()); - meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex - basalmeltingrate[iv]+=meltingrate_enthalpy; + connectivity=IssmDouble(vertices[iv]->Connectivity()); + meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex + basalmeltingrate[iv]+=meltingrate_enthalpy; this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); if(reCast(dt)) - watercolumn[iv]+=dt*meltingrate_enthalpy; + watercolumn[iv]+=dt*meltingrate_enthalpy; else - watercolumn[iv]+=meltingrate_enthalpy; + watercolumn[iv]+=meltingrate_enthalpy; } - /*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)); + /*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 gaussup; - delete friction; + delete friction; } /*}}}*/ @@ -5429,39 +5425,39 @@ /*Intermediaries*/ - bool isenthalpy; - IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES]; - IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; - IssmDouble latentheat, dt; - IssmDouble connectivity; - GaussPenta* gauss; - - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + bool isenthalpy; + IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES]; + IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; + IssmDouble latentheat, dt; + IssmDouble connectivity; + GaussPenta* gauss; + + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - - /*Check wether enthalpy is activated*/ + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + + /*Check wether enthalpy is activated*/ parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); if(!isenthalpy) return; - - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - latentheat=matpar->GetLatentHeat(); - - gauss=new GaussPenta(); - for(int iv=0;ivGaussVertex(iv); - enthalpy_input->GetInputValue(&enthalpy[iv], gauss); - pressure_input->GetInputValue(&pressure[iv], gauss); - matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); - - /*drain water fraction & update enthalpy*/ - // TODO: replace connectivity-wise draining by draining once per node per timestep - connectivity=IssmDouble(vertices[iv]->Connectivity()); - waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity; - matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); - } - /*feed updated results back into model*/ - this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); - this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); - // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum)); // temperature should not change during drainage... + + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + latentheat=matpar->GetLatentHeat(); + + gauss=new GaussPenta(); + for(int iv=0;ivGaussVertex(iv); + enthalpy_input->GetInputValue(&enthalpy[iv], gauss); + pressure_input->GetInputValue(&pressure[iv], gauss); + matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); + + /*drain water fraction & update enthalpy*/ + // TODO: replace connectivity-wise draining by draining once per node per timestep + connectivity=IssmDouble(vertices[iv]->Connectivity()); + waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity; + matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); + } + /*feed updated results back into model*/ + this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); + this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); + // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum)); // temperature should not change during drainage... } /*}}}*/