Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 17026) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 17027) @@ -246,5 +246,5 @@ /*Enthalpy diffusion parameter*/ - IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.); + IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.); /* Start looping on the number of gaussian points: */ @@ -443,5 +443,5 @@ if(stabilization==2){ diameter=element->MinEdgeLength(xyz_list); - kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.); + kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>=0.); } @@ -829,10 +829,4 @@ } - /*drain excess water fraction: */ - //for(int i=0;ielements->Size();i++){ - // element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - // element->DrainWaterfraction(); - //} - /*Update basal dirichlet BCs for enthalpy: */ for(int i=0;ielements->Size();i++){ @@ -841,9 +835,4 @@ } }/*}}}*/ - - - - - void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/ @@ -878,5 +867,5 @@ Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); - IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>0.); + IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>=0.); element->NormalBase(&normal_base[0],xyz_list_base); element->VerticalSegmentIndices(&pairindices,&numsegments); @@ -990,5 +979,5 @@ basalmeltingrate[vertexdown]=meltingrate_enthalpy[is]; watercolumn[vertexdown]+=meltingrate_enthalpy[is]; - } + } _assert_(watercolumn[vertexdown]>=0.); } @@ -1002,12 +991,13 @@ delete gauss; delete friction; - xDelete(enthalpy); - xDelete(pressure); - xDelete(watercolumn); - xDelete(basalmeltingrate); + xDelete(enthalpy); + xDelete(pressure); + xDelete(watercolumn); + xDelete(basalmeltingrate); xDelete(meltingrate_enthalpy); xDelete(heating); xDelete(drainrate_column); xDelete(drainrate_element); + xDelete(xyz_list); }/*}}}*/ @@ -1070,82 +1060,72 @@ }/*}}}*/ - - - - void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/ - /* /\*Intermediary*\/ */ - /* bool isdynamicbasalspc,setspc; */ - /* int numindices, numindicesup; */ - /* IssmDouble pressure, pressureup; */ - /* IssmDouble h_pmp, enthalpy, enthalpyup; */ - /* IssmDouble watercolumn; */ - /* int *indices = NULL, *indicesup = NULL; */ + /*Intermediary*/ + bool isdynamicbasalspc,setspc; + int numindices, numindicesup; + IssmDouble pressure, pressureup; + IssmDouble h_pmp, enthalpy, enthalpyup; + IssmDouble watercolumn; + int *indices = NULL, *indicesup = NULL; + Node* node = NULL; - /* /\* Only update Constraints at the base of grounded ice*\/ */ - /* if(!IsOnBed() || IsFloating()) return; */ - - /* /\*Check wether dynamic basal boundary conditions are activated *\/ */ - /* element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); */ - /* if(!isdynamicbasalspc) return; */ - - /* /\*Fetch indices of basal & surface nodes for this finite element*\/ */ - /* // BasalNodeIndices(&numindices,&indices,this->element_type); */ - /* // SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); */ - /* // _assert_(numindices==numindicesup); */ - - /* /\*Get parameters and inputs: *\/ */ - /* Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); */ - /* Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); */ - /* Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); */ - - /* /\*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*\/ */ - /* GaussPenta* gauss=new GaussPenta(); */ - /* GaussPenta* gaussup=new GaussPenta(); */ - /* for(int i=0;iGaussNode(this->element_type,indices[i]); */ - /* gaussup->GaussNode(this->element_type,indicesup[i]); */ - - /* /\*Check wether there is a temperate layer at the base or not *\/ */ - /* /\*check if node is temperate, if not, continue*\/ */ - /* enthalpy_input->GetInputValue(&enthalpy, gauss); */ - /* pressure_input->GetInputValue(&pressure, gauss); */ - /* watercolumn_input->GetInputValue(&watercolumn,gauss); */ - /* setspc = false; */ - /* // TODO: add case H0.; */ - /* if (enthalpy>=element->PureIceEnthalpy(pressure)){ */ - /* /\*check if upper node is temperate, too. */ - /* if yes, then we have a temperate layer of positive thickness and reset the spc. */ - /* if not, apply dirichlet BC.*\/ */ - /* enthalpy_input->GetInputValue(&enthalpyup, gaussup); */ - /* pressure_input->GetInputValue(&pressureup, gaussup); */ - /* setspc=((enthalpyupPureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; */ - /* } */ - - /* if (setspc) { */ - /* /\*Calculate enthalpy at pressure melting point *\/ */ - /* h_pmp=matpar->PureIceEnthalpy(pressure); */ - /* /\*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*\/ */ - /* //element->ApplyConstraint(indices[i],1,h_pmp); */ - /* //nodes[indices[i]]->ApplyConstraint(1,h_pmp); */ - /* } */ - /* else { */ - /* /\*remove spc*\/ */ - /* //element->RemoveConstraint(indices[i],0); */ - /* //nodes[indices[i]]->DofInFSet(0); */ - /* } */ - /* } */ - - /* /\*Free ressources:*\/ */ - /* xDelete(indices); */ - /* xDelete(indicesup); */ - /* delete gauss; */ - /* delete gaussup; */ -}/*}}}*/ - - - - + /* Only update Constraints at the base of grounded ice*/ + if(!(element->IsOnBed()) || element->IsFloating()) return; + + /*Check wether dynamic basal boundary conditions are activated */ + element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); + if(!isdynamicbasalspc) return; + + /*Fetch indices of basal & surface nodes for this finite element*/ + Penta *penta = (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element* + penta->BasalNodeIndices(&numindices,&indices,element->GetElementType()); + penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType()); + _assert_(numindices==numindicesup); + + /*Get parameters and inputs: */ + Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); + Input* enthalpy_input=element->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* watercolumn_input=element->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + + /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/ + GaussPenta* gauss=new GaussPenta(); + GaussPenta* gaussup=new GaussPenta(); + for(int i=0;iGaussNode(element->GetElementType(),indices[i]); + gaussup->GaussNode(element->GetElementType(),indicesup[i]); + + /*Check wether there is a temperate layer at the base or not */ + /*check if node is temperate, else continue*/ + enthalpy_input->GetInputValue(&enthalpy, gauss); + pressure_input->GetInputValue(&pressure, gauss); + watercolumn_input->GetInputValue(&watercolumn,gauss); + h_pmp=PureIceEnthalpy(element,pressure); + if (enthalpy>=h_pmp){ + /*check if upper node is temperate, too. + if yes, then we have a temperate layer of positive thickness and reset the spc. + if not, apply dirichlet BC.*/ + enthalpy_input->GetInputValue(&enthalpyup, gaussup); + pressure_input->GetInputValue(&pressureup, gaussup); + setspc=((enthalpyup=0.))?true:false; + } + else if(watercolumn>0.) // case H0. + setspc=true; + else + setspc = false; + + node=element->GetNode(indices[i]); + if (setspc) + node->ApplyConstraint(1,h_pmp); /*apply spc*/ //nodes[indices[i]]->ApplyConstraint(1,h_pmp); + else + node->DofInFSet(0); /*remove spc*/ //nodes[indices[i]]->DofInFSet(0); + } + + /*Free ressources:*/ + xDelete(indices); + xDelete(indicesup); + delete gauss; + delete gaussup; +}/*}}}*/ /*Intermediaries*/ @@ -1167,7 +1147,6 @@ int iv; IssmDouble lambda; /* fraction of cold ice */ - IssmDouble kappa ,kappa_c,kappa_t; /* enthalpy conductivities */ + IssmDouble kappa,kappa_c,kappa_t; /* enthalpy conductivities */ IssmDouble Hc,Ht; - /*Get pressures and enthalpies on vertices*/ @@ -1209,6 +1188,6 @@ _assert_((Hc+Ht)>0.); lambda = Hc/(Hc+Ht); - kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t); - } + kappa = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1 + } /*Clean up and return*/ @@ -1218,6 +1197,5 @@ xDelete(enthalpies); return kappa; -} -/*}}}*/ +}/*}}}*/ IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17026) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17027) @@ -168,4 +168,5 @@ virtual IssmDouble GetYcoord(Gauss* gauss)=0; virtual void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; + virtual int GetElementType(void)=0; virtual IssmDouble SurfaceArea(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17026) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17027) @@ -141,4 +141,5 @@ IssmDouble GetYcoord(Gauss* gauss){_error_("Not implemented");}; IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented yet");}; + int GetElementType(void){_error_("not implemented yet");}; Gauss* NewGauss(void){_error_("not implemented yet");}; Gauss* NewGauss(int order); Index: /issm/trunk-jpl/src/c/cores/thermal_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/thermal_core.cpp (revision 17026) +++ /issm/trunk-jpl/src/c/cores/thermal_core.cpp (revision 17027) @@ -6,4 +6,5 @@ #include "../toolkits/toolkits.h" #include "../classes/classes.h" +#include "../analyses/analyses.h" #include "../shared/shared.h" #include "../modules/modules.h" @@ -17,4 +18,5 @@ int solution_type,numoutputs; char** requested_outputs = NULL; + EnthalpyAnalysis * enthalpy_analysis = NULL; /*first recover parameters common to all solutions*/ @@ -47,5 +49,7 @@ /*Post process*/ - if(solution_type!=SteadystateSolutionEnum) PostprocessingEnthalpyx(femmodel); + enthalpy_analysis = new EnthalpyAnalysis(); + enthalpy_analysis->PostProcessing(femmodel); + delete enthalpy_analysis; } else{ @@ -66,3 +70,4 @@ /*Free ressources:*/ if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} + } Index: /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 17026) +++ /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 17027) @@ -53,5 +53,5 @@ if (enthalpy < H_sp){ Tstar = referencetemperature + enthalpy/heatcapacity - betaCC*pressure; - waterfraction = 0; + waterfraction = 0.; } else{ @@ -63,14 +63,13 @@ /*Get A*/ - if(Tstar<263.15){ - A=3.61*pow(10.,-13.) * exp( -6.*pow(10.,4.)/(R*Tstar)); + if(Tstar<=263.15){ + A=3.61e-13 * exp( -6.e+4/(R*Tstar)); } else{ - A=1.73*pow(10., 3.) * exp(-13.9*pow(10.,4.)/(R*Tstar)); + A=1.73e3 * exp(-13.9e+4/(R*Tstar)); } - A*=(1 + 181.25*waterfraction); + A*=(1. + 181.25*waterfraction); /*Convert to B*/ - _assert_(n>0); B=pow(A,-1./n); Index: /issm/trunk-jpl/src/m/materials/arrhenius.m =================================================================== --- /issm/trunk-jpl/src/m/materials/arrhenius.m (revision 17026) +++ /issm/trunk-jpl/src/m/materials/arrhenius.m (revision 17027) @@ -31,8 +31,12 @@ end -wf_max=.1; -if(waterfraction>wf_max) +wf_max=1.; +if(any(waterfraction>wf_max)) error(['waterfraction exceeds permitted maximum of ' num2str(wf_max) '.']); end + +%limit waterfraction to 1% +pos1p=find(waterfraction>0.01); +waterfraction(pos1p)=0.01; pos=find((temperature0)); % cold, wet ice @@ -71,4 +75,7 @@ Qa=GetQa(T); A0=GetA0(T); + if(w>0.01) + w=0.01; + end A=A0.*exp(-Qa./(R*T)).*(1+181.25*w); end