Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 11346) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 11347) @@ -161,4 +161,5 @@ ThermalSpctemperatureEnum, ThermalStabilizationEnum, + ThermalIsenthalpyEnum, ThicknessEnum, TimesteppingCflCoefficientEnum, Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11347) @@ -165,4 +165,5 @@ case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; case ThermalStabilizationEnum : return "ThermalStabilization"; + case ThermalIsenthalpyEnum : return "ThermalIsenthalpy"; case ThicknessEnum : return "Thickness"; case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 11347) @@ -79,4 +79,5 @@ parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum)); parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum)); + parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum)); parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 11347) @@ -20,5 +20,5 @@ int i,analysis_type,dim,verbose; - bool isthermal,isprognostic,isdiagnostic,isgroundingline; + bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy; /*output: */ @@ -39,4 +39,5 @@ iomodel->Constant(&verbose,VerboseEnum); iomodel->Constant(&isthermal,TransientIsthermalEnum); + iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum); iomodel->Constant(&isprognostic,TransientIsprognosticEnum); iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum); @@ -52,6 +53,11 @@ if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && dim==2) continue; if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && dim==2) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && dim==2) continue; if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isthermal==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isthermal==false) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isthermal==false) continue; + if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isenthalpy==true) continue; + if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isenthalpy==true) continue; + if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isenthalpy==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==PrognosticAnalysisEnum && isprognostic==false && isgroundingline==false) continue; if(solution_type==TransientSolutionEnum && analysis_type==DiagnosticHorizAnalysisEnum && isdiagnostic==false) continue; Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11347) @@ -163,4 +163,5 @@ else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum; + else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum; else if (strcmp(name,"Thickness")==0) return ThicknessEnum; else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11347) @@ -4085,8 +4085,10 @@ double rho_ice,heatcapacity,geothermalflux_value; double basalfriction,alpha2,vx,vy; - double scalar; + double scalar,enthalpy,enthalpyup; + double pressure,pressureup; double basis[NUMVERTICES]; Friction* friction=NULL; GaussPenta* gauss=NULL; + GaussPenta* gaussup=NULL; /* Geothermal flux on ice sheet base and basal friction */ @@ -4106,4 +4108,6 @@ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); @@ -4113,25 +4117,43 @@ /* Start looping on the number of gauss 2d (nodes on the bedrock) */ gauss=new GaussPenta(0,1,2,2); + gaussup=new GaussPenta(3,4,5,2); for(ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); + gaussup->GaussPoint(ig); GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); GetNodalFunctionsP1(&basis[0], gauss); - 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)/(heatcapacity*rho_ice); - if(dt) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; + enthalpy_input->GetInputValue(&enthalpy,gauss); + pressure_input->GetInputValue(&pressure,gauss); +// if(enthalpy>matpar->PureIceEnthalpy(pressure)){ +// enthalpy_input->GetInputValue(&enthalpyup,gaussup); +// pressure_input->GetInputValue(&pressureup,gaussup); +// if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){ +// //do nothing, don't add heatflux +// } +// else{ +// //need to change spcenthalpy according to Aschwanden +// //NEED TO UPDATE +// } +// } +// else{ + 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)/(heatcapacity*rho_ice); + if(dt) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*basis[i]; +// } } /*Clean up and return*/ delete gauss; + delete gaussup; delete friction; return pe; @@ -4309,5 +4331,5 @@ matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]); if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector"); - if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector"); + // if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector"); } Index: /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp (revision 11347) @@ -420,5 +420,5 @@ } else{ - return 0.1*thermalconductivity/(rho_ice*heatcapacity); + return 1*thermalconductivity/(rho_ice*heatcapacity); } } Index: /issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp (revision 11347) @@ -95,5 +95,5 @@ case TransientSolutionEnum: - numanalyses=8; + numanalyses=9; analyses=(int*)xmalloc(numanalyses*sizeof(int)); analyses[0]=DiagnosticHorizAnalysisEnum; @@ -105,4 +105,5 @@ analyses[6]=MeltingAnalysisEnum; analyses[7]=PrognosticAnalysisEnum; + analyses[8]=EnthalpyAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 11347) @@ -24,5 +24,5 @@ /*parameters: */ double finaltime,dt,yts; - bool control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline; + bool control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy; bool dakota_analysis=false; bool time_adapt=false; @@ -51,4 +51,5 @@ femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); @@ -91,5 +92,10 @@ _printf_(VerboseSolution()," computing temperatures:\n"); #ifdef _HAVE_THERMAL_ - thermal_core_step(femmodel,step,time); + if(isenthalpy==0){ + thermal_core_step(femmodel,step,time); + } + else{ + enthalpy_core_step(femmodel,step,time); + } #else _error_("ISSM was not compiled with thermal capabilities. Exiting"); @@ -135,5 +141,7 @@ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum,step,time); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,step,time); + if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time); InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum,step,time); Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 11346) +++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 11347) @@ -85,4 +85,7 @@ if(count>=max_nonlinear_iterations){ _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + converged=true; + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); break; } Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 11346) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 11347) @@ -62,5 +62,5 @@ checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]); end - if ismember(EnthalpyAnalysisEnum,analyses), + if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy), checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]); end Index: /issm/trunk-jpl/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.m (revision 11346) +++ /issm/trunk-jpl/src/m/classes/thermal.m (revision 11347) @@ -12,4 +12,5 @@ penalty_lock = 0; penalty_factor = 0; + isenthalpy = 0; end methods @@ -42,4 +43,7 @@ %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor obj.penalty_factor=3; + + %Should we use cold ice (default) or enthalpy formulation + obj.isenthalpy=0; end % }}} function checkconsistency(obj,md,solution,analyses) % {{{ @@ -52,4 +56,5 @@ if ismember(EnthalpyAnalysisEnum,analyses), checkfield(md,'thermal.spctemperature','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*md.geometry.thickness,'message','spctemperature should be below the adjusted melting point'); + checkfield(md,'thermal.isenthalpy','numel',1,'values',[0 1]); end end % }}} @@ -62,4 +67,5 @@ fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)'); + fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); end % }}} @@ -71,4 +77,5 @@ WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean'); end % }}} end Index: /issm/trunk-jpl/src/m/enum/ThermalIsenthalpyEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/ThermalIsenthalpyEnum.m (revision 11347) +++ /issm/trunk-jpl/src/m/enum/ThermalIsenthalpyEnum.m (revision 11347) @@ -0,0 +1,11 @@ +function macro=ThermalIsenthalpyEnum() +%THERMALISENTHALPYENUM - Enum of ThermalIsenthalpy +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=ThermalIsenthalpyEnum() + +macro=StringToEnum('ThermalIsenthalpy'); Index: /issm/trunk-jpl/src/m/solutions/AnalysisConfiguration.m =================================================================== --- /issm/trunk-jpl/src/m/solutions/AnalysisConfiguration.m (revision 11346) +++ /issm/trunk-jpl/src/m/solutions/AnalysisConfiguration.m (revision 11347) @@ -42,6 +42,6 @@ case TransientSolutionEnum, - numanalyses=8; - analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum]; + numanalyses=9; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum;EnthalpyAnalysisEnum]; case FlaimSolutionEnum, Index: /issm/trunk-jpl/src/m/solutions/transient_core.m =================================================================== --- /issm/trunk-jpl/src/m/solutions/transient_core.m (revision 11346) +++ /issm/trunk-jpl/src/m/solutions/transient_core.m (revision 11347) @@ -19,4 +19,5 @@ isthermal=femmodel.parameters.TransientIsthermal; isgroundingline=femmodel.parameters.TransientIsgroundingline; + isenthalpy=femmodel.parameters.ThermalIsenthalpy; groundinglinemigration=femmodel.parameters.GroundinglineMigration; @@ -57,5 +58,9 @@ if (isthermal & dim==3) issmprintf(VerboseSolution,'\n%s',[' computing temperature']); - femmodel=thermal_core_step(femmodel); + if (isenthalpy==0), + femmodel=thermal_core_step(femmodel); + else + femmodel=enthalpy_core_step(femmodel); + end end @@ -90,4 +95,6 @@ femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time); if (dim==3 & isthermal), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,step,time);end + if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,WaterfractionEnum,step,time);end + if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,EnthalpyEnum,step,time);end femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time); femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceforcingsMassBalanceEnum,step,time);