Changeset 18620

Timestamp:

10/13/14 08:55:00 (10 years ago)

Author:

jbondzio

Message:

CHG: dynamic basal constraints for enthalpy steadystate

Location:

issm/trunk-jpl

Files:

• src/c/analyses/EnthalpyAnalysis.cpp (modified) (2 diffs)
• src/c/analyses/EnthalpyAnalysis.h (modified) (1 diff)
• src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (modified) (4 diffs)
• src/m/classes/thermal.m (modified) (1 diff)
• src/m/classes/thermal.py (modified) (1 diff)
• test/Archives/Archive431.nc (modified) ( previous)
• test/Archives/Archive432.nc (modified) ( previous)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -1185 +1185 @@
                 UpdateBasalConstraints(element);
         }
+        femmodel->UpdateConstraintsx();
 }/*}}}*/
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
@@ -1394 +1395 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
 
-        if(dt==0.){ // steadystate case
-                state=0; //TODO: add consistent steadystate basal condition scheme
-//              if(enthalpy<PureIceEnthalpy(element,pressure)){ /*is base cold?*/
-//                      if(watercolumn<=0.) state=0; // cold, dry base
-//                      else state=1; // cold, wet base (refreezing)
-//              }
-//              else{ /*base is temperate, check if upper node is temperate, too.*/
-//                      if(enthalpyup<PureIceEnthalpy(element,pressureup))
-//                              if(meltingrate<0.)      state=2; // refreezing temperate base (non-physical, only for steadystate solver)
-//                              else                                            state=3; // melting temperate base with no temperate layer
-//                      else
-//                              state=4; // melting temperate base with temperate layer of positive thickness
-//              }
-        }
-        else{ // transient case
-                if(enthalpy<PureIceEnthalpy(element,pressure)){
-                        if(watercolumn<=0.) state=0; // cold, dry base
-                        else state=1; // cold, wet base (refreezing)
-                }
-                else{
-                        if(enthalpyup<PureIceEnthalpy(element,pressureup))      state=3; // temperate base, but no temperate layer
-                        else state=4; // temperate layer with positive thickness
-                }
+        if(enthalpy<PureIceEnthalpy(element,pressure)){
+                if(watercolumn<=0.) state=0; // cold, dry base
+                else state=1; // cold, wet base (refreezing)
+        }
+        else{
+                if(enthalpyup<PureIceEnthalpy(element,pressureup)){
+                        if((dt==0.) && (meltingrate<0.)) state=2;       // refreezing temperate base (non-physical, only for steadystate solver)
+                        else    state=3; // temperate base, but no temperate layer
+                }
+                else state=4; // temperate layer with positive thickness
         }
 

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

@@ -51 +51 @@
                 static void UpdateBasalConstraintsTransient(Element* element);
                 static void UpdateBasalConstraintsSteadystate(Element* element);
-//              static int GetThermalBasalCondition(Element* element, Gauss* gauss, Gauss* gaussup, int enthalpy_enum);
-//              static int GetThermalBasalCondition(Element* element, GaussPenta* gauss, GaussPenta* gaussup, int enthalpy_enum);
                 static int GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
                 static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp

@@ -50 +50 @@
                 femmodel->UpdateConstraintsx();
 
-                //Update the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+                //Update the solution to make sure that tf and tf_old are similar (for next step in transient or steadystate)
                 GetSolutionFromInputsx(&tg,femmodel);
                 Reducevectorgtofx(&tf, tg, femmodel->nodes,femmodel->parameters);
@@ -65 +65 @@
         for(;;){
                 delete tf_old;tf_old=tf;
+
+                if(isenthalpy){ 
+                        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+                        /*Update old solution, such that sizes of tf_old and tf are comparable*/
+                        if(isdynamicbasalspc)   Reducevectorgtofx(&tf_old, tg, femmodel->nodes,femmodel->parameters);
+                }
+                else SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
                 delete tg;
-
-                if(isenthalpy) SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
-                else SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
                 Reduceloadx(pf, Kfs, ys); delete Kfs;
@@ -82 +86 @@
                 if(VerboseConvergence()) _printf0_("   number of unstable constraints: " << num_unstable_constraints << "\n");
 
-                if(isenthalpy){
-                        /*Increase count: */
+                if(isenthalpy){ // enthalpy method
+                        IssmDouble dt;
+                        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
                         count++;
+                        bool max_iteration_state=false;
+                        if(count>=thermal_maxiter){
+                                _printf0_("   maximum number of nonlinear iterations (" << thermal_maxiter << ") exceeded\n"); 
+                                converged=true;
+                                InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+                                InputUpdateFromSolutionx(femmodel,tg);          
+                                max_iteration_state=true;
+                        }
                         if(converged==true){
-                                bool max_iteration_state=false;
                                 int step; IssmDouble time;
                                 femmodel->parameters->FindParam(&time,TimeEnum);
@@ -93 +106 @@
                                 break;
                         }
-                        if(count>=thermal_maxiter){
-                                _printf0_("   maximum number of nonlinear iterations (" << thermal_maxiter << ") exceeded\n"); 
-                                converged=true;
-                                InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
-                                InputUpdateFromSolutionx(femmodel,tg);          
-                                bool max_iteration_state=true;
-                                int step; IssmDouble time;
-                                femmodel->parameters->FindParam(&time,TimeEnum);
-                                femmodel->parameters->FindParam(&step,StepEnum);
-                                femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, step, time));
-                                break;
+                        else if(dt==0.){
+                                EnthalpyAnalysis::ComputeBasalMeltingrate(femmodel);
+                                EnthalpyAnalysis::UpdateBasalConstraints(femmodel);
                         }
                 }
-                else{
+                else{ // dry ice method
                         if(!converged){
                                 if(num_unstable_constraints<=thermal_penalty_threshold) converged=true;

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

@@ -100 +100 @@
                                 if(md.thermal.isenthalpy)
                                         if isnan(md.stressbalance.reltol),
-                                                md = checkmessage(md,['for a steadystate computation, stressbalance.reltol (relative convergence criterion) must be defined!']);
-                                        end
+                                                md = checkmessage(md,['for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!']);
+                                        end 
                                         md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message','reltol must be larger than zero');
-                                        md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[1], 'message',['for enthalpy run thermal.isdynamicbasalspc should be 1']);
                                 end
             end

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

@@ -95 +95 @@
                         if(md.thermal.isenthalpy):
                                 if numpy.isnan(md.stressbalance.reltol):
-                                        md.checkmessage("for a steadystate computation, stressbalance.reltol (relative convergence criterion) must be defined!")
+                                        md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!")
                                 md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero");
-                                md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel', [1],'values',[1], 'message',"for enthalpy run thermal.isdynamicbasalspc should be true")
                 md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1)