Changeset 24082


Ignore:
Timestamp:
07/11/19 11:37:13 (6 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added md.toolkits.RecoveryAnalysis, which is a recovery mode when default solver fails

Location:
issm/trunk-jpl/src
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/ParseToolkitsOptionsx/ParseToolkitsOptionsx.cpp

    r22844 r24082  
    2222        char word1[1000];
    2323        char word2[1000];
    24         int my_rank;
    25         int i;
    2624
    2725        /*intermediary: */
    28         int         *analyses     = NULL;
    29         char       **strings      = NULL;
    30         char        *string       = NULL;
    31         char       **toolkits     = NULL;
    32         char        *toolkit      = NULL;
    33         char        *newstring    = NULL;
    34         char        *catstring    = NULL;
    35         int          numanalyses;
    36         int          stringlength,toolkitlength;
     26        int   *analyses     = NULL;
     27        char **strings      = NULL;
     28        char  *string       = NULL;
     29        char **toolkits     = NULL;
     30        char  *toolkit      = NULL;
     31        char  *newstring    = NULL;
     32        char  *catstring    = NULL;
     33        int    numanalyses;
     34        int    stringlength,toolkitlength;
    3735
    3836        /*Get my_rank:*/
    39         my_rank=IssmComm::GetRank();
     37        int my_rank=IssmComm::GetRank();
    4038
    4139        if(my_rank==0){
    … …  
    5957                strings  = xNew<char*>(numanalyses);
    6058                toolkits = xNew<char*>(numanalyses);
    61                 for(i=0;i<numanalyses;i++) strings[i]  = NULL;
    62                 for(i=0;i<numanalyses;i++) toolkits[i] = NULL;
     59                for(int i=0;i<numanalyses;i++) strings[i]  = NULL;
     60                for(int i=0;i<numanalyses;i++) toolkits[i] = NULL;
    6361
    6462                /*Go back to beginning of file:*/
    … …  
    119117        }
    120118        ISSM_MPI_Bcast(analyses,numanalyses,ISSM_MPI_INT,0,IssmComm::GetComm());
    121         for(i=0;i<numanalyses;i++){
     119        for(int i=0;i<numanalyses;i++){
    122120                char* toolkit=toolkits[i];
    123121                if(my_rank==0){
    … …  
    147145
    148146        /*Clean up and return*/
    149         for(i=0;i<numanalyses;i++) xDelete<char>(strings[i]);
    150         for(i=0;i<numanalyses;i++) xDelete<char>(toolkits[i]);
     147        for(int i=0;i<numanalyses;i++) xDelete<char>(strings[i]);
     148        for(int i=0;i<numanalyses;i++) xDelete<char>(toolkits[i]);
    151149        xDelete<char*>(strings);
    152150        xDelete<char*>(toolkits);

  • issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp

    r15104 r24082  
    1111#include "./Solverx.h"
    1212#include "../../shared/shared.h"
     13#include "../../classes/Params/Parameters.h"
    1314
    1415void    Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){
    1516
    16         /*intermediary: */
    17         Solver<IssmDouble> *solver=NULL;
    18 
    19         /*output: */
    20         Vector<IssmDouble> *uf=NULL;
    21 
    22         if(VerboseModule()) _printf0_("   Solving matrix system\n");
    23 
    24         /*Initialize solver: */
    25         solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
     17        /*Create Solver Object*/
     18        Solver<IssmDouble>* solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
    2619
    2720        /*Solve:*/
    28         uf=solver->Solve();
     21        if(VerboseModule()) _printf0_("   Solving matrix system\n");
     22        Vector<IssmDouble>* uf=solver->Solve();
     23
     24        /*Check convergence, if failed, try recovery model*/
     25        if(!checkconvergence(Kff,pf,uf,parameters)){
     26
     27                _printf0_("   WARNING: Solver failed, Trying Recovery Mode\n");
     28                ToolkitsOptionsFromAnalysis(parameters,RecoveryAnalysisEnum);
     29                delete uf;
     30                uf=solver->Solve();
     31
     32                if(!checkconvergence(Kff,pf,uf,parameters)) _error_("Recovery solver failed...");
     33        }
    2934
    3035        /*clean up and assign output pointers:*/
     36        _assert_(puf);
    3137        delete solver;
    3238        *puf=uf;
    3339}
     40bool checkconvergence(Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Parameters* parameters){
     41
     42        /*Recover parameters: */
     43        IssmDouble solver_residue_threshold;
     44        parameters->FindParam(&solver_residue_threshold,SettingsSolverResidueThresholdEnum);
     45
     46        /*don't check convergence if NaN*/
     47        if(xIsNan<IssmPDouble>(solver_residue_threshold)) return true;
     48
     49        /*compute KUF = KU - F = K*U - F*/
     50        Vector<IssmDouble>* KU  = uf->Duplicate(); Kff->MatMult(uf,KU);
     51        Vector<IssmDouble>* KUF = KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.);
     52        delete KU;
     53
     54        /*compute norm(KUF), norm(F)*/
     55        IssmDouble nKUF=KUF->Norm(NORM_TWO);
     56        IssmDouble nF=pf->Norm(NORM_TWO);
     57        delete KUF;
     58
     59        /*Check solver residue*/
     60        IssmDouble solver_residue=nKUF/nF;
     61        if(VerboseConvergence()) _printf0_("\n   solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n");
     62        if(xIsNan<IssmDouble>(solver_residue)) _error_("Solver residue is NaN");
     63
     64        /*clean up*/
     65
     66        /*Check convergence*/
     67        if(solver_residue>solver_residue_threshold){
     68                _printf0_("   solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << " => Trying recovery solver\n");
     69                return false;
     70        }
     71        else{
     72                return true;
     73        }
     74
     75}

  • issm/trunk-jpl/src/c/modules/Solverx/Solverx.h

    r15051 r24082  
    1616/* local prototypes: */
    1717void    Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters);
     18bool checkconvergence(Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Parameters* parameters);
    1819
    1920#endif  /* _SOLVERX_H */

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24080 r24082  
    11591159syn keyword cConstant RegionaloutputEnum
    11601160syn keyword cConstant RegularEnum
     1161syn keyword cConstant RecoveryAnalysisEnum
    11611162syn keyword cConstant RiftfrontEnum
    11621163syn keyword cConstant SIAApproximationEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r24080 r24082  
    11571157        RegionaloutputEnum,
    11581158        RegularEnum,
     1159        RecoveryAnalysisEnum,
    11591160        RiftfrontEnum,
    11601161        SIAApproximationEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r24080 r24082  
    11611161                case RegionaloutputEnum : return "Regionaloutput";
    11621162                case RegularEnum : return "Regular";
     1163                case RecoveryAnalysisEnum : return "RecoveryAnalysis";
    11631164                case RiftfrontEnum : return "Riftfront";
    11641165                case SIAApproximationEnum : return "SIAApproximation";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r24080 r24082  
    11881188              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
    11891189              else if (strcmp(name,"Regular")==0) return RegularEnum;
     1190              else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    11901191              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    11911192              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    … …  
    12431244              else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
    12441245              else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
    1245               else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
     1249              if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
     1250              else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
    12501251              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    12511252              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;

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

    r24080 r24082  
    1313
    1414        /*intermediary*/
    15         Vector<IssmDouble>* KU=NULL;
    16         Vector<IssmDouble>* KUF=NULL;
    1715        Vector<IssmDouble>* KUold=NULL;
    1816        Vector<IssmDouble>* KUoldF=NULL;
    1917        Vector<IssmDouble>* duf=NULL;
    2018        IssmDouble ndu,nduinf,nu;
    21         IssmDouble nKUF;
    2219        IssmDouble nKUoldF;
    2320        IssmDouble nF;
    24         IssmDouble solver_residue,res;
     21        IssmDouble res;
    2522        int analysis_type;
    2623
    … …  
    3229                *pconverged=true;
    3330                return;
    34         }
    35 
    36 
    37         /*Display solver caracteristics*/
    38         if (VerboseConvergence()){
    39 
    40        
    41 
    42                 /*compute KUF = KU - F = K*U - F*/
    43                 KU=uf->Duplicate(); Kff->MatMult(uf,KU);
    44                 KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
    45 
    46                 /*compute norm(KUF), norm(F) and residue*/
    47                 nKUF=KUF->Norm(NORM_TWO);
    48                 nF=pf->Norm(NORM_TWO);
    49                 solver_residue=nKUF/nF;
    50                 _printf0_("\n" << "   solver residue: norm(KU-F)/norm(F)=" << solver_residue << "\n");
    51                 if(xIsNan<IssmDouble>(solver_residue)){
    52                         //Kff->Echo();
    53                 }
    54 
    55                 /*clean up*/
    56                 delete KU;
    57                 delete KUF;
    5831        }
    5932

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp

    r23587 r24082  
    1818        Vector<IssmDouble>*  df  = NULL;
    1919        Vector<IssmDouble>*  ys  = NULL;
    20         IssmDouble solver_residue_threshold;
    21 
    22         /*solver convergence test: */
    23         IssmDouble nKUF;
    24         IssmDouble nF;
    25         IssmDouble solver_residue;
    26         Vector<IssmDouble>* KU=NULL;
    27         Vector<IssmDouble>* KUF=NULL;
    2820
    2921        /*Recover parameters: */
    30         femmodel->parameters->FindParam(&solver_residue_threshold,SettingsSolverResidueThresholdEnum);
    3122        femmodel->UpdateConstraintsx();
    3223        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    … …  
    3829        femmodel->profiler->Stop(SOLVER);
    3930
    40         /*Check that the solver converged nicely: */
    41 
    42         //compute KUF = KU - F = K*U - F
    43         KU=uf->Duplicate(); Kff->MatMult(uf,KU);
    44         KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
    45 
    46         //compute norm(KUF), norm(F) and residue
    47         nKUF=KUF->Norm(NORM_TWO);
    48         nF=pf->Norm(NORM_TWO);
    49         solver_residue=nKUF/nF;
    50 
    51         if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_("   solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n");
    52 
    53         //clean up
    54         delete KU; delete KUF;
     31        /*Clean up*/
    5532        delete Kff; delete pf; delete df;
    5633

  • issm/trunk-jpl/src/m/classes/toolkits.js

    r21065 r24082  
    2929                        }
    3030                }
     31
     32                this.RecoveryAnalysis = this.DefaultAnalysis;
    3133        }// }}}
    3234        this.disp = function(){// {{{
    … …  
    9496        //properties
    9597        // {{{
    96         this.DefaultAnalysis           = [];
     98        this.DefaultAnalysis  = [];
     99        this.RecoveryAnalysis = [];
    97100        //The other properties are dynamic
    98101        this.setdefaultparameters();

  • issm/trunk-jpl/src/m/classes/toolkits.m

    r21049 r24082  
    66classdef toolkits < dynamicprops
    77    properties (SetAccess=public)
    8                  DefaultAnalysis           = struct();
     8                 DefaultAnalysis  = struct();
     9                 RecoveryAnalysis = struct();
    910                 %The other properties are dynamic
    1011         end
    … …  
    5455                         end
    5556
     57                         %Use same solver for Recovery mode
     58                         self.RecoveryAnalysis = self.DefaultAnalysis;
     59
     60
    5661                 end % }}}
    5762                 function disp(self) % {{{
    … …  
    7378                 end % }}}
    7479                 function ToolkitsFile(toolkits,filename) % {{{
    75                 %TOOLKITSFILE - build toolkits file
    76                 %
    77                 %   Build a Petsc compatible options file, from the toolkits model field  + return options string.
    78                 %   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
    79                 %
    80                 %   Usage:     ToolkitsFile(toolkits,filename);
     80                        %TOOLKITSFILE - build toolkits file
     81                        %
     82                        %   Build a Petsc compatible options file, from the toolkits model field  + return options string.
     83                        %   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
     84                        %
     85                        %   Usage:     ToolkitsFile(toolkits,filename);
    8186
    8287                         %open file for writing
    … …  
    122127                         fclose(fid);
    123128                 end %}}}
    124                 function savemodeljs(self,fid,modelname) % {{{
    125        
    126                         analyses=properties(self);
     129                 function savemodeljs(self,fid,modelname) % {{{
     130
     131                         analyses=properties(self);
    127132                         for i=1:numel(analyses),
    128133                                 if isempty(fieldnames(self.(analyses{i})))
    … …  
    132137                                 end
    133138                         end
    134        
    135                 end % }}}
     139
     140                 end % }}}
    136141         end
    137142 end

  • issm/trunk-jpl/src/m/classes/toolkits.py

    r23716 r24082  
    3131                                raise IOError("ToolkitsFile error: need at least Mumps or Gsl to define issm solver type")
    3232
    33                 #The other properties are dynamic
     33                #Use same solver for Recovery mode
     34                self.RecoveryAnalysis = self.DefaultAnalysis
     35
     36                #The other properties are dynamic
    3437        # }}}
    3538        def __repr__(self):    # {{{
Note: See TracChangeset for help on using the changeset viewer.