solutionsequence_thermal_nonlinear.cpp File Reference

: core of the thermal solution More...

#include "./solutionsequences.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

Go to the source code of this file.

Functions
void	solutionsequence_thermal_nonlinear (FemModel *femmodel)

Detailed Description

: core of the thermal solution

Definition in file solutionsequence_thermal_nonlinear.cpp.

Function Documentation

solutionsequence_thermal_nonlinear()

void solutionsequence_thermal_nonlinear

(

FemModel *

femmodel

)

Definition at line 11 of file solutionsequence_thermal_nonlinear.cpp.

                                                           {
 
   /*solution : */
   Vector<IssmDouble>* tg=NULL; 
   Vector<IssmDouble>* tf=NULL; 
   Vector<IssmDouble>* tf_old=NULL; 
   Vector<IssmDouble>* ys=NULL; 
   IssmDouble melting_offset;
 
   /*intermediary: */
   Matrix<IssmDouble>* Kff=NULL;
   Matrix<IssmDouble>* Kfs=NULL;
   Vector<IssmDouble>* pf=NULL;
   Vector<IssmDouble>* df=NULL;
 
   bool converged;
   bool isenthalpy, isdynamicbasalspc;
   int constraints_converged;
   int num_unstable_constraints;
   int count;
   int thermal_penalty_threshold;
   int thermal_maxiter;
 
   /*parameters:*/
   IssmDouble eps_rel;
 
   /*Recover parameters: */
   femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
   femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
 
   converged=false;
   InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
 
   if(isenthalpy){
      femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
      femmodel->parameters->FindParam(&eps_rel,ThermalReltolEnum);
      femmodel->UpdateConstraintsx();
 
      //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);
      InputUpdateFromSolutionx(femmodel,tg);
   }
   else{
      femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
      ResetZigzagCounterx(femmodel);
      femmodel->UpdateConstraintsx();
   }
 
   count=1;
 
   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){
            delete tf_old;
            Reducevectorgtofx(&tf_old, tg, femmodel->nodes,femmodel->parameters);
         }
      }
      else SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
      delete tg;
      CreateNodalConstraintsx(&ys,femmodel->nodes);
      Reduceloadx(pf, Kfs, ys); delete Kfs;
      femmodel->profiler->Start(SOLVER);
      Solverx(&tf, Kff, pf, tf_old, df, femmodel->parameters);
      femmodel->profiler->Stop(SOLVER);
      Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
      if(isenthalpy){ 
         convergence(&converged,Kff,pf,tf,tf_old,0.05,eps_rel,NAN); 
         InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
      }
      delete Kff; delete pf; delete df;
      InputUpdateFromSolutionx(femmodel,tg);
      ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
      if(VerboseConvergence()) _printf0_("   number of unstable constraints: " << num_unstable_constraints << "\n");
 
      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){
            break;
         }
         else if(dt==0.){
            EnthalpyAnalysis::ComputeBasalMeltingrate(femmodel);
            EnthalpyAnalysis::UpdateBasalConstraints(femmodel);
         }
      }
      else{ // dry ice method
         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;
      }
   }
 
   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: */
   delete tg;
   delete tf;
   delete tf_old;
}