Index: /issm/trunk-jpl/src/c/modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp (revision 24082) @@ -22,20 +22,18 @@ char word1[1000]; char word2[1000]; - int my_rank; - int i; /*intermediary: */ - int *analyses = NULL; - char **strings = NULL; - char *string = NULL; - char **toolkits = NULL; - char *toolkit = NULL; - char *newstring = NULL; - char *catstring = NULL; - int numanalyses; - int stringlength,toolkitlength; + int *analyses = NULL; + char **strings = NULL; + char *string = NULL; + char **toolkits = NULL; + char *toolkit = NULL; + char *newstring = NULL; + char *catstring = NULL; + int numanalyses; + int stringlength,toolkitlength; /*Get my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); if(my_rank==0){ @@ -59,6 +57,6 @@ strings = xNew(numanalyses); toolkits = xNew(numanalyses); - for(i=0;i(strings[i]); - for(i=0;i(toolkits[i]); + for(int i=0;i(strings[i]); + for(int i=0;i(toolkits[i]); xDelete(strings); xDelete(toolkits); Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 24082) @@ -11,23 +11,65 @@ #include "./Solverx.h" #include "../../shared/shared.h" +#include "../../classes/Params/Parameters.h" void Solverx(Vector** puf, Matrix* Kff, Vector* pf, Vector* uf0,Vector* df, Parameters* parameters){ - /*intermediary: */ - Solver *solver=NULL; - - /*output: */ - Vector *uf=NULL; - - if(VerboseModule()) _printf0_(" Solving matrix system\n"); - - /*Initialize solver: */ - solver=new Solver(Kff,pf,uf0,df,parameters); + /*Create Solver Object*/ + Solver* solver=new Solver(Kff,pf,uf0,df,parameters); /*Solve:*/ - uf=solver->Solve(); + if(VerboseModule()) _printf0_(" Solving matrix system\n"); + Vector* 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* Kff,Vector* pf,Vector* uf,Parameters* parameters){ + + /*Recover parameters: */ + IssmDouble solver_residue_threshold; + parameters->FindParam(&solver_residue_threshold,SettingsSolverResidueThresholdEnum); + + /*don't check convergence if NaN*/ + if(xIsNan(solver_residue_threshold)) return true; + + /*compute KUF = KU - F = K*U - F*/ + Vector* KU = uf->Duplicate(); Kff->MatMult(uf,KU); + Vector* 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=nKUF/nF; + if(VerboseConvergence()) _printf0_("\n solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n"); + if(xIsNan(solver_residue)) _error_("Solver residue is NaN"); + + /*clean up*/ + + /*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; + } + +} Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 24081) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 24082) @@ -16,4 +16,5 @@ /* local prototypes: */ void Solverx(Vector** puf, Matrix* Kff, Vector* pf, Vector* uf0,Vector* df, Parameters* parameters); +bool checkconvergence(Matrix* Kff,Vector* pf,Vector* uf,Parameters* parameters); #endif /* _SOLVERX_H */ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24081) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24082) @@ -1159,4 +1159,5 @@ syn keyword cConstant RegionaloutputEnum syn keyword cConstant RegularEnum +syn keyword cConstant RecoveryAnalysisEnum syn keyword cConstant RiftfrontEnum syn keyword cConstant SIAApproximationEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24081) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24082) @@ -1157,4 +1157,5 @@ RegionaloutputEnum, RegularEnum, + RecoveryAnalysisEnum, RiftfrontEnum, SIAApproximationEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24082) @@ -1161,4 +1161,5 @@ case RegionaloutputEnum : return "Regionaloutput"; case RegularEnum : return "Regular"; + case RecoveryAnalysisEnum : return "RecoveryAnalysis"; case RiftfrontEnum : return "Riftfront"; case SIAApproximationEnum : return "SIAApproximation"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24082) @@ -1188,4 +1188,5 @@ else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else if (strcmp(name,"Regular")==0) return RegularEnum; + else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum; else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum; - else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; + if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; + else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 24082) @@ -13,14 +13,11 @@ /*intermediary*/ - Vector* KU=NULL; - Vector* KUF=NULL; Vector* KUold=NULL; Vector* KUoldF=NULL; Vector* duf=NULL; IssmDouble ndu,nduinf,nu; - IssmDouble nKUF; IssmDouble nKUoldF; IssmDouble nF; - IssmDouble solver_residue,res; + IssmDouble res; int analysis_type; @@ -32,28 +29,4 @@ *pconverged=true; return; - } - - - /*Display solver caracteristics*/ - if (VerboseConvergence()){ - - - - /*compute KUF = KU - F = K*U - F*/ - KU=uf->Duplicate(); Kff->MatMult(uf,KU); - KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0); - - /*compute norm(KUF), norm(F) and residue*/ - nKUF=KUF->Norm(NORM_TWO); - nF=pf->Norm(NORM_TWO); - solver_residue=nKUF/nF; - _printf0_("\n" << " solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n"); - if(xIsNan(solver_residue)){ - //Kff->Echo(); - } - - /*clean up*/ - delete KU; - delete KUF; } Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 24081) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 24082) @@ -18,15 +18,6 @@ Vector* df = NULL; Vector* ys = NULL; - IssmDouble solver_residue_threshold; - - /*solver convergence test: */ - IssmDouble nKUF; - IssmDouble nF; - IssmDouble solver_residue; - Vector* KU=NULL; - Vector* KUF=NULL; /*Recover parameters: */ - femmodel->parameters->FindParam(&solver_residue_threshold,SettingsSolverResidueThresholdEnum); femmodel->UpdateConstraintsx(); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); @@ -38,19 +29,5 @@ femmodel->profiler->Stop(SOLVER); - /*Check that the solver converged nicely: */ - - //compute KUF = KU - F = K*U - F - KU=uf->Duplicate(); Kff->MatMult(uf,KU); - KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0); - - //compute norm(KUF), norm(F) and residue - nKUF=KUF->Norm(NORM_TWO); - nF=pf->Norm(NORM_TWO); - solver_residue=nKUF/nF; - - if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_(" solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n"); - - //clean up - delete KU; delete KUF; + /*Clean up*/ delete Kff; delete pf; delete df; Index: /issm/trunk-jpl/src/m/classes/toolkits.js =================================================================== --- /issm/trunk-jpl/src/m/classes/toolkits.js (revision 24081) +++ /issm/trunk-jpl/src/m/classes/toolkits.js (revision 24082) @@ -29,4 +29,6 @@ } } + + this.RecoveryAnalysis = this.DefaultAnalysis; }// }}} this.disp = function(){// {{{ @@ -94,5 +96,6 @@ //properties // {{{ - this.DefaultAnalysis = []; + this.DefaultAnalysis = []; + this.RecoveryAnalysis = []; //The other properties are dynamic this.setdefaultparameters(); Index: /issm/trunk-jpl/src/m/classes/toolkits.m =================================================================== --- /issm/trunk-jpl/src/m/classes/toolkits.m (revision 24081) +++ /issm/trunk-jpl/src/m/classes/toolkits.m (revision 24082) @@ -6,5 +6,6 @@ classdef toolkits < dynamicprops properties (SetAccess=public) - DefaultAnalysis = struct(); + DefaultAnalysis = struct(); + RecoveryAnalysis = struct(); %The other properties are dynamic end @@ -54,4 +55,8 @@ end + %Use same solver for Recovery mode + self.RecoveryAnalysis = self.DefaultAnalysis; + + end % }}} function disp(self) % {{{ @@ -73,10 +78,10 @@ end % }}} function ToolkitsFile(toolkits,filename) % {{{ - %TOOLKITSFILE - build toolkits file - % - % Build a Petsc compatible options file, from the toolkits model field + return options string. - % This file will also be used when the toolkit used is 'issm' instead of 'petsc' - % - % Usage: ToolkitsFile(toolkits,filename); + %TOOLKITSFILE - build toolkits file + % + % Build a Petsc compatible options file, from the toolkits model field + return options string. + % This file will also be used when the toolkit used is 'issm' instead of 'petsc' + % + % Usage: ToolkitsFile(toolkits,filename); %open file for writing @@ -122,7 +127,7 @@ fclose(fid); end %}}} - function savemodeljs(self,fid,modelname) % {{{ - - analyses=properties(self); + function savemodeljs(self,fid,modelname) % {{{ + + analyses=properties(self); for i=1:numel(analyses), if isempty(fieldnames(self.(analyses{i}))) @@ -132,6 +137,6 @@ end end - - end % }}} + + end % }}} end end Index: /issm/trunk-jpl/src/m/classes/toolkits.py =================================================================== --- /issm/trunk-jpl/src/m/classes/toolkits.py (revision 24081) +++ /issm/trunk-jpl/src/m/classes/toolkits.py (revision 24082) @@ -31,5 +31,8 @@ raise IOError("ToolkitsFile error: need at least Mumps or Gsl to define issm solver type") - #The other properties are dynamic + #Use same solver for Recovery mode + self.RecoveryAnalysis = self.DefaultAnalysis + + #The other properties are dynamic # }}} def __repr__(self): # {{{