Changeset 18619


Ignore:
Timestamp:
10/13/14 08:46:11 (10 years ago)
Author:
Mathieu Morlighem
Message:

CHG: convergence criteria must now be provided to convergence

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

Legend:

Unmodified
Added
Removed

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

    r18601 r18619  
    77#include "../shared/shared.h"
    88
    9 void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,Parameters* parameters){
     9void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){
    1010
    1111        /*output*/
    … …  
    2525        int analysis_type;
    2626
    27         /*convergence options*/
    28         IssmDouble eps_res;
    29         IssmDouble eps_rel;
    30         IssmDouble eps_abs;
    31         IssmDouble yts;
    32 
    3327        if(VerboseModule()) _printf0_("   checking convergence\n");
    3428
    … …  
    3933                return;
    4034        }
    41 
    42         /*get convergence options*/
    43         parameters->FindParam(&eps_res,StressbalanceRestolEnum);
    44         // get analysis type
    45         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    46         if(analysis_type==StressbalanceAnalysisEnum)
    47                 parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
    48         else if(analysis_type==EnthalpyAnalysisEnum)
    49                 parameters->FindParam(&eps_rel,ThermalReltolEnum);
    50         else
    51                 eps_rel = NAN;
    52 
    53         parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
    54         parameters->FindParam(&yts,ConstantsYtsEnum);
    5535
    5636        /*Display solver caracteristics*/
    … …  
    139119                //print
    140120                if (!xIsNan<IssmDouble>(eps_abs)){
    141                         if ((nduinf*yts)<eps_abs){
    142                                 if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " < " << eps_abs << " m/yr\n");
     121                        if ((nduinf)<eps_abs){
     122                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf << " < " << eps_abs << "\n");
    143123                        }
    144124                        else{
    145                                 if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " > " << eps_abs << " m/yr\n");
     125                                if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf << " > " << eps_abs << "\n");
    146126                                converged=false;
    147127                        }
    148128                }
    149                 else  _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " m/yr\n");
     129                else  _printf0_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf << "\n");
    150130
    151131        }

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

    r17816 r18619  
    3030        /*parameters:*/
    3131        int max_nonlinear_iterations;
    32         int  configuration_type;
     32        int configuration_type;
     33        IssmDouble eps_res,eps_rel,eps_abs;
    3334
    3435        /*Recover parameters: */
    … …  
    3637        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    3738        femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
     39        femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
     40        femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
     41        femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
    3842        femmodel->UpdateConstraintsx();
    3943
    … …  
    8387
    8488                /*Check convergence*/
    85                 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
     89                convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
    8690                delete Kff; delete pf;
    8791                if(converged==true){   

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

    r17055 r18619  
    3131        int max_nonlinear_iterations;
    3232        int configuration_type;
     33        IssmDouble eps_res,eps_rel,eps_abs;
     34
    3335
    3436        /*Recover parameters: */
    3537        femmodel->parameters->FindParam(&min_mechanical_constraints,StressbalanceRiftPenaltyThresholdEnum);
    3638        femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
     39        femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
     40        femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
     41        femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
    3742        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    3843        femmodel->UpdateConstraintsx();
    … …  
    6671                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
    6772
    68                 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); delete Kff; delete pf; delete df;
     73                convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
    6974                InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    7075                InputUpdateFromSolutionx(femmodel,ug);

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

    r16126 r18619  
    3333        int  max_nonlinear_iterations;
    3434        int  configuration_type;
     35        IssmDouble eps_res,eps_rel,eps_abs;
    3536
    3637        /*Recover parameters: */
    3738        femmodel->parameters->FindParam(&min_mechanical_constraints,StressbalanceRiftPenaltyThresholdEnum);
    3839        femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
     40        femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
     41        femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
     42        femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
    3943        femmodel->UpdateConstraintsx();
    4044
    … …  
    6872                InputUpdateFromSolutionx(femmodel,ug_horiz);
    6973
    70                 convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); delete Kff_horiz; delete pf_horiz; delete df_horiz;
     74                convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,eps_res,eps_rel,eps_abs); delete Kff_horiz; delete pf_horiz; delete df_horiz;
    7175
    7276                /*Second compute vertical velocity: */

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

    r18590 r18619  
    3131        int thermal_penalty_threshold;
    3232        int thermal_maxiter;
    33         IssmDouble thermal_reltol;
    3433
    3534        /*parameters:*/
    3635        int  configuration_type;
     36        IssmDouble eps_rel;
    3737
    3838        /*Recover parameters: */
    … …  
    4040        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    4141        femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
     42        femmodel->parameters->FindParam(&eps_rel,ThermalReltolEnum);
    4243
    4344        converged=false;
    … …  
    4647        if(isenthalpy){
    4748                femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
    48                 femmodel->parameters->FindParam(&thermal_reltol,ThermalReltolEnum);
     49                femmodel->parameters->FindParam(&eps_rel,ThermalReltolEnum);
    4950                femmodel->UpdateConstraintsx();
    5051
    … …  
    7374                Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
    7475                if(isenthalpy){
    75                         convergence(&converged,Kff,pf,tf,tf_old,femmodel->parameters);
     76                        convergence(&converged,Kff,pf,tf,tf_old,0.05,eps_rel,NAN);
    7677                        InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
    7778                }

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h

    r18344 r18619  
    2424
    2525/*convergence*/
    26 void convergence(bool* pconverged, Matrix<IssmDouble>* K_ff,Vector<IssmDouble>* p_f,Vector<IssmDouble>* u_f,Vector<IssmDouble>* u_f_old,Parameters* parameters);
     26void convergence(bool* pconverged, Matrix<IssmDouble>* K_ff,Vector<IssmDouble>* p_f,Vector<IssmDouble>* u_f,Vector<IssmDouble>* u_f_old,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs);
    2727
    2828#endif

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

    r18579 r18619  
    223223                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','restol','format','Double');
    224224                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','reltol','format','Double');
    225                         WriteData(fid,'object',obj,'class','stressbalance','fieldname','abstol','format','Double');
     225                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','abstol','format','Double','scale',1./yts);
    226226                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','isnewton','format','Integer');
    227227                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','FSreconditioning','format','Double');

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

    r17686 r18619  
    173173                WriteData(fid,'object',self,'class','stressbalance','fieldname','restol','format','Double')
    174174                WriteData(fid,'object',self,'class','stressbalance','fieldname','reltol','format','Double')
    175                 WriteData(fid,'object',self,'class','stressbalance','fieldname','abstol','format','Double')
     175                WriteData(fid,'object',self,'class','stressbalance','fieldname','abstol','format','Double','scale',1./yts)
    176176                WriteData(fid,'object',self,'class','stressbalance','fieldname','isnewton','format','Integer')
    177177                WriteData(fid,'object',self,'class','stressbalance','fieldname','FSreconditioning','format','Double')
Note: See TracChangeset for help on using the changeset viewer.