Changeset 24313 for issm/trunk/src/c/modules/Solverx/Solverx.cpp

Timestamp:

11/01/19 12:01:57 (5 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 24310

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/modules/Solverx/Solverx.cpp (modified) (1 diff)

issm/trunk/src/c

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota

issm/trunk/src/c/modules/Solverx/Solverx.cpp

@@ -11 +11 @@
 #include "./Solverx.h"
 #include "../../shared/shared.h"
+#include "../../classes/Params/Parameters.h"
 
 void    Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){
 
-        /*intermediary: */
-        Solver<IssmDouble> *solver=NULL;
-
-        /*output: */
-        Vector<IssmDouble> *uf=NULL;
-
-        if(VerboseModule()) _printf0_("   Solving matrix system\n");
-
-        /*Initialize solver: */
-        solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
+        /*Create Solver Object*/
+        Solver<IssmDouble>* solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
 
         /*Solve:*/
-        uf=solver->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;
 }
+bool checkconvergence(Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Parameters* parameters){
+
+        /*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 << " => Trying recovery solver\n");
+                return false;
+        }
+        else{
+                return true;
+        }
+
+}

@@ -5 +5 @@
 *.obj
 *.exe
+*.mod
 appscan.*
 ouncemake*
@@ -14 +15 @@
 probe.results
 stXXXX*
-*.deps
-*.dirstamp
+.deps
+.dirstamp
 .libs
 issm
@@ -21 +22 @@
 issm_slr
 issm_ocean
-lnb_param.mod
-lovenb_sub.mod
-model.mod
-util.mod
+issm_dakota