Changeset 18589

Timestamp:

10/08/14 02:06:27 (10 years ago)

Author:

jbondzio

Message:

CHG: use solutionsequence_thermal_nonlinear for enthalpy

Location:

issm/trunk-jpl/src/c

Files:

• analyses/EnthalpyAnalysis.cpp (modified) (1 diff)
• cores/thermal_core.cpp (modified) (1 diff)
• solutionsequences/convergence.cpp (modified) (2 diffs)
• solutionsequences/solutionsequence_thermal_nonlinear.cpp (modified) (4 diffs)

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

@@ -15 +15 @@
 
         parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalMaxiterEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(ThermalReltolEnum));
         parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
         parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));

issm/trunk-jpl/src/c/cores/thermal_core.cpp ¶

@@ -31 +31 @@
                 if(VerboseSolution()) _printf0_("   computing enthalpy\n");
                 femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
-                solutionsequence_nonlinear(femmodel,true);
+                solutionsequence_thermal_nonlinear(femmodel);
 
                 /*transfer enthalpy to enthalpy picard for the next step: */

issm/trunk-jpl/src/c/solutionsequences/convergence.cpp ¶

@@ -23 +23 @@
         IssmDouble nF;
         IssmDouble solver_residue,res;
+        int analysis_type;
 
         /*convergence options*/
@@ -41 +42 @@
         /*get convergence options*/
         parameters->FindParam(&eps_res,StressbalanceRestolEnum);
-        parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
+        // get analysis type
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        if(analysis_type==StressbalanceAnalysisEnum)
+                parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
+        else if(analysis_type==EnthalpyAnalysisEnum)
+                parameters->FindParam(&eps_rel,ThermalReltolEnum);
         parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
         parameters->FindParam(&yts,ConstantsYtsEnum);

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

@@ -3 +3 @@
  */ 
 
+#include "./solutionsequences.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
@@ -24 +25 @@
 
         bool converged;
+        bool isenthalpy, isdynamicbasalspc;
         int constraints_converged;
         int num_unstable_constraints;
@@ -29 +31 @@
         int thermal_penalty_threshold;
         int thermal_maxiter;
+        IssmDouble thermal_reltol;
 
         /*parameters:*/
@@ -34 +37 @@
 
         /*Recover parameters: */
-        femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
+        femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
         femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
 
+        converged=false;
+        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+
+        if(isenthalpy){
+                femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+                femmodel->parameters->FindParam(&thermal_reltol,ThermalReltolEnum);
+                femmodel->UpdateConstraintsx();
+
+                //Update the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+                GetSolutionFromInputsx(&tg,femmodel);
+                Reducevectorgtofx(&tf, tg, femmodel->nodes,femmodel->parameters);
+                InputUpdateFromSolutionx(femmodel,tg);
+        }
+        else{
+                femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
+                InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
+                femmodel->UpdateConstraintsx();
+        }
+
         count=1;
-        converged=false;
+        
+        for(;;){
+                delete tf_old;tf_old=tf;
+                delete tg;
 
-        InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
-        InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
-        femmodel->UpdateConstraintsx();
-
-        for(;;){
-
-                delete tf_old; tf_old=tf;
-                SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
+                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;
-                Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
-                delete Kff;delete pf;delete tg; delete df;
+                Solverx(&tf, Kff, pf, tf_old, df, femmodel->parameters);
                 Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
+                if(isenthalpy){ 
+                        convergence(&converged,Kff,pf,tf,tf_old,femmodel->parameters); delete Kff; delete pf; delete df;
+                        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+                }
                 InputUpdateFromSolutionx(femmodel,tg);
+                ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
+                if(VerboseConvergence()) _printf0_("   number of unstable constraints: " << num_unstable_constraints << "\n");
 
-                ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
-
-                if (!converged){
-                        if(VerboseConvergence()) _printf0_("   #unstable constraints = " << num_unstable_constraints << "\n");
-                        if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
-                        if (count>=thermal_maxiter){
+                if(isenthalpy){
+                        /*Increase count: */
+                        count++;
+                        if(converged==true){
+                                bool max_iteration_state=false;
+                                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;
+                        }
+                        if(count>=thermal_maxiter){
+                                _printf0_("   maximum number of nonlinear iterations (" << thermal_maxiter << ") exceeded\n"); 
                                 converged=true;
-                                _printf0_("   maximum number of iterations (" << thermal_maxiter << ") exceeded\n"); 
+                                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;
                         }
                 }
-                count++;
-
-                InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
-
-                if(converged)break;
+                else{
+                        if(!converged){
+                                if(num_unstable_constraints<=thermal_penalty_threshold) converged=true;
+                                if(count>=thermal_maxiter){
+                                        converged=true;
+                                        _printf0_("   maximum number of iterations (" << thermal_maxiter << ") exceeded\n"); 
+                                }
+                        }
+                        count++;
+                        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+                        if(converged)break;
+                }
         }
 
-        InputUpdateFromSolutionx(femmodel,tg);
-        femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);
+        if(isenthalpy){
+                if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
+        }
+        else{
+                InputUpdateFromSolutionx(femmodel,tg);
+                femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);
+        }
 
         /*Free ressources: */