Changeset 11347

Timestamp:

02/07/12 14:30:43 (13 years ago)

Author:

seroussi

Message:

added possibility to used enthalpy in transient solutions

Location:

issm/trunk-jpl/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/ModelProcessorx.cpp (modified) (3 diffs)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (4 diffs)
• c/objects/Materials/Matpar.cpp (modified) (1 diff)
• c/solutions/AnalysisConfiguration.cpp (modified) (2 diffs)
• c/solutions/transient_core.cpp (modified) (4 diffs)
• c/solvers/solver_nonlinear.cpp (modified) (1 diff)
• m/classes/initialization.m (modified) (1 diff)
• m/classes/thermal.m (modified) (5 diffs)
• m/enum/ThermalIsenthalpyEnum.m (added)
• m/solutions/AnalysisConfiguration.m (modified) (1 diff)
• m/solutions/transient_core.m (modified) (3 diffs)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h ¶

@@ -161 +161 @@
         ThermalSpctemperatureEnum,
         ThermalStabilizationEnum,
+        ThermalIsenthalpyEnum,
         ThicknessEnum,
         TimesteppingCflCoefficientEnum,

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp ¶

@@ -165 +165 @@
                 case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
                 case ThermalStabilizationEnum : return "ThermalStabilization";
+                case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
                 case ThicknessEnum : return "Thickness";
                 case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient";

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -79 +79 @@
         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));

issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp ¶

@@ -20 +20 @@
 
         int   i,analysis_type,dim,verbose;
-        bool  isthermal,isprognostic,isdiagnostic,isgroundingline;
+        bool  isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy;
         
         /*output: */
@@ -39 +39 @@
         iomodel->Constant(&verbose,VerboseEnum);
         iomodel->Constant(&isthermal,TransientIsthermalEnum);
+        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
         iomodel->Constant(&isprognostic,TransientIsprognosticEnum);
         iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum);
@@ -52 +53 @@
                 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;

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp ¶

@@ -163 +163 @@
         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;

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp ¶

@@ -4085 +4085 @@
         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 +4108 @@
         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 +4117 @@
         /* 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();ig<gauss->end();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;i<numdof;i++) pe->values[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;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+//              }
         }
 
         /*Clean up and return*/
         delete gauss;
+        delete gaussup;
         delete friction;
         return pe;
@@ -4309 +4331 @@
                         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");
                 }
                         

issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp ¶

@@ -420 +420 @@
         }
         else{
-                return 0.1*thermalconductivity/(rho_ice*heatcapacity);
+                return 1*thermalconductivity/(rho_ice*heatcapacity);
         }
 }

issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp ¶

@@ -95 +95 @@
 
                 case TransientSolutionEnum:
-                        numanalyses=8;
+                        numanalyses=9;
                         analyses=(int*)xmalloc(numanalyses*sizeof(int));
                         analyses[0]=DiagnosticHorizAnalysisEnum;
@@ -105 +105 @@
                         analyses[6]=MeltingAnalysisEnum;
                         analyses[7]=PrognosticAnalysisEnum;
+                        analyses[8]=EnthalpyAnalysisEnum;
                         break;
                 

issm/trunk-jpl/src/c/solutions/transient_core.cpp ¶

@@ -24 +24 @@
         /*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 +51 @@
         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 +92 @@
                         _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 +141 @@
                         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);

issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp ¶

@@ -85 +85 @@
                 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;
                 }

issm/trunk-jpl/src/m/classes/initialization.m ¶

@@ -62 +62 @@
                                 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

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

@@ -12 +12 @@
                 penalty_lock      = 0;
                 penalty_factor    = 0;
+                isenthalpy        = 0;
         end
         methods
@@ -42 +43 @@
                         %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 +56 @@
                         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 +67 @@
                         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 +77 @@
                         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

issm/trunk-jpl/src/m/solutions/AnalysisConfiguration.m ¶

@@ -42 +42 @@
 
         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,

issm/trunk-jpl/src/m/solutions/transient_core.m ¶

@@ -19 +19 @@
         isthermal=femmodel.parameters.TransientIsthermal;
         isgroundingline=femmodel.parameters.TransientIsgroundingline;
+        isenthalpy=femmodel.parameters.ThermalIsenthalpy;
         groundinglinemigration=femmodel.parameters.GroundinglineMigration;
 
@@ -57 +58 @@
                 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 +95 @@
                         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);