Solverx.h File Reference

solver More...

#include "../../toolkits/toolkits.h"

Go to the source code of this file.

Functions
void	Solverx (Vector< IssmDouble > **puf, Matrix< IssmDouble > *Kff, Vector< IssmDouble > *pf, Vector< IssmDouble > *uf0, Vector< IssmDouble > *df, Parameters *parameters)
bool	checkconvergence (Matrix< IssmDouble > *Kff, Vector< IssmDouble > *pf, Vector< IssmDouble > *uf, Parameters *parameters)

Detailed Description

solver

Definition in file Solverx.h.

Function Documentation

Solverx()

void Solverx	(	Vector< IssmDouble > **	puf,
		Matrix< IssmDouble > *	Kff,
		Vector< IssmDouble > *	pf,
		Vector< IssmDouble > *	uf0,
		Vector< IssmDouble > *	df,
		Parameters *	parameters
	)

Definition at line 15 of file Solverx.cpp.

                                                                                                                                                                {
 
   /*Create Solver Object*/
   Solver<IssmDouble>* solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
 
   /*Solve:*/
   if(VerboseModule()) _printf0_("   Solving matrix system\n");
   Vector<IssmDouble>* uf=solver->Solve();
 
   /*Check convergence, if failed, try recovery model*/
   if(!checkconvergence(Kff,pf,uf,parameters)){
 
      _printf0_("WARNING: Solver failed, Trying Recovery Mode\n");
      ToolkitsOptionsFromAnalysis(parameters,RecoveryAnalysisEnum);
      delete uf;
      uf=solver->Solve();
 
      if(!checkconvergence(Kff,pf,uf,parameters)) _error_("Recovery solver failed...");
   }
 
   /*clean up and assign output pointers:*/
   _assert_(puf);
   delete solver;
   *puf=uf;
}

checkconvergence()

bool checkconvergence	(	Matrix< IssmDouble > *	Kff,
		Vector< IssmDouble > *	pf,
		Vector< IssmDouble > *	uf,
		Parameters *	parameters
	)

Definition at line 40 of file Solverx.cpp.

                                                                                                                   {
 
   /*Recover parameters: */
   IssmDouble solver_residue_threshold;
   parameters->FindParam(&solver_residue_threshold,SettingsSolverResidueThresholdEnum);
 
   /*don't check convergence if NaN*/
   if(xIsNan<IssmDouble>(solver_residue_threshold)) return true;
 
   /*compute KUF = KU - F = K*U - F*/
   Vector<IssmDouble>* KU  = uf->Duplicate(); Kff->MatMult(uf,KU);
   Vector<IssmDouble>* KUF = KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.);
   delete KU;
 
   /*compute norm(KUF), norm(F)*/
   IssmDouble nKUF=KUF->Norm(NORM_TWO);
   IssmDouble nF=pf->Norm(NORM_TWO);
   delete KUF;
 
   /*Check solver residue*/
   IssmDouble solver_residue = 0.;
   if(nF>0.)  solver_residue = nKUF/(nF);
   if(VerboseConvergence()) _printf0_("\n   solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n");
   if(xIsNan<IssmDouble>(solver_residue)) _error_("Solver residue is NaN");
 
   /*Check convergence*/
   if(solver_residue>solver_residue_threshold){
      _printf0_("solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << " > "<<solver_residue_threshold<<" (md.settings.solver_residue_threshold)");
      return false;
   }
   else{
      return true;
   }
 
}