Changeset 1666


Ignore:
Timestamp:
08/13/09 12:51:46 (16 years ago)
Author:
Mathieu Morlighem
Message:

Added mechanical equilibrium residue convergence criterion

Location:
issm/trunk/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

    r1648 r1666  
    7070        parameters->AddObject(param);
    7171
     72        /*eps_res: */
     73        count++;
     74        param= new Param(count,"eps_res",DOUBLE);
     75        param->SetDouble(model->eps_res);
     76        parameters->AddObject(param);
     77
    7278        /*eps_rel: */
    7379        count++;

  • issm/trunk/src/c/ModelProcessorx/Model.cpp

    r1648 r1666  
    128128        model->debug=0;
    129129        model->plot=0;
     130        model->eps_res=0;
    130131        model->eps_rel=0;
    131132        model->eps_abs=0;
    … …  
    346347        ModelFetchData((void**)&model->mincontrolconstraint,NULL,NULL,model_handle,"mincontrolconstraint","Scalar",NULL);
    347348        ModelFetchData((void**)&model->maxcontrolconstraint,NULL,NULL,model_handle,"maxcontrolconstraint","Scalar",NULL);
     349        ModelFetchData((void**)&model->eps_res,NULL,NULL,model_handle,"eps_res","Scalar",NULL);
    348350        ModelFetchData((void**)&model->eps_rel,NULL,NULL,model_handle,"eps_rel","Scalar",NULL);
    349351        ModelFetchData((void**)&model->eps_abs,NULL,NULL,model_handle,"eps_abs","Scalar",NULL);

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r1648 r1666  
    126126        int     debug;
    127127        int     plot;
     128        double  eps_res;
    128129        double  eps_rel;
    129130        double  eps_abs;

  • issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp

    r1665 r1666  
    4242        double nKUoldF;
    4343        double nF;
    44         double solver_residue,mechanical_residue;
     44        double solver_residue,res;
    4545
    4646        /*parameters:*/
    … …  
    4848        char* solver_string=NULL;
    4949        int debug=0;
    50         double eps_rel,eps_abs,yts;
     50        double eps_res,eps_rel,eps_abs,yts;
    5151        int dofs[4]={1,1,0,0}; //recover vx,vy by default. vz and pressure may be.
    5252
    … …  
    5858        fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    5959        fem->parameters->FindParam((void*)&solver_string,"solverstring");
     60        fem->parameters->FindParam((void*)&eps_res,"eps_res");
    6061        fem->parameters->FindParam((void*)&eps_rel,"eps_rel");
    6162        fem->parameters->FindParam((void*)&eps_abs,"eps_abs");
    … …  
    117118                if (debug) _printf_("   solving\n");
    118119                Solverx(&uf, Kff, pf, old_uf, solver_string);
     120
     121                /*print solver residue = norm(KU-F)/norm(F)*/
     122                if (debug>1){
     123                        //compute KUF = KU - F = K*U - F
     124                        VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU);
     125                        VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf);
     126                        //compute norm(KUF), norm(F) and residue
     127                        VecNorm(KUF,NORM_2,&nKUF);
     128                        VecNorm(pf,NORM_2,&nF);
     129                        solver_residue=nKUF/nF;
     130                        //clean up
     131                        VecFree(&KU);
     132                        VecFree(&KUF);
     133                        _printf_("%s%g\n","   solver residue: norm(K[n]U[n]-F)/norm(F)",solver_residue);
     134                }
    119135               
    120136                if (debug) _printf_("   merging solution from f to g set\n");
    … …  
    130146                /*Figure out if convergence is reached.*/
    131147
    132                 /*solver residue = norm(KU-F)/norm(F)*/
    133                 //compute KUF = KU - F = K*U - F
    134                 VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU);
    135                 VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf);
    136                 //compute norm(KUF), norm(F) and residue
    137                 VecNorm(KUF,NORM_2,&nKUF);
    138                 VecNorm(pf,NORM_2,&nF);
    139                 solver_residue=nKUF/nF;
    140                 //clean up
    141                 VecFree(&KU);
    142                 VecFree(&KUF);
    143                 if (debug) _printf_("%-50s%g\n","   Convergence criterion: norm(KU-F)/norm(F)",solver_residue);
    144 
    145                 /*mechanical residue = norm(KUold-F)/norm(F)*/
    146                 //compute KUoldF = KUold - F = K*Uold - F
     148                /*1 mechanical equilibrium residue = norm(K[n]U[n-1]-F)/norm(F)*/
     149                //compute K[n]U[n-1]F = K[n]U[n-1] - F
    147150                VecDuplicate(uf,&KUold); MatMultPatch(Kff,old_uf,KUold);
    148151                VecDuplicate(KUold,&KUoldF);VecCopy(KUold,KUoldF); VecAYPX(KUoldF,-1.0,pf);
    149152                //compute norm(KUF), norm(F) and residue
    150153                VecNorm(KUoldF,NORM_2,&nKUoldF);
    151                 mechanical_residue=nKUoldF/nF;
     154                res=nKUoldF/nF;
    152155                //clean up
    153156                VecFree(&KUold);
    154157                VecFree(&KUoldF);
    155                 if (debug) _printf_("%-50s%g%s\n","   Convergence criterion: norm(KUold-F)/norm(F)",mechanical_residue*100," \%");
     158                if (isnan(res)) throw ErrorException(__FUNCT__,exprintf("mechanical equilibrium convergence criterion is NaN! "));
     159                if((res)<eps_res){
     160                        if (debug) _printf_("%-50s%g%s%g%s\n","   mechanical equilibrium convergence criterion",res*100," < ",eps_res*100," \%");
     161                        converged=1;
     162                }
     163                else{
     164                        if (debug) _printf_("%-50s%g%s%g%s\n","   mechanical equilibrium convergence criterion",res*100," > ",eps_res*100," \%");
     165                        converged=1;//TTTTTTTTTTTTTTTTTTFTYICGWFYUOGCYEUOGCYEUWGCHLWGWI
     166                }
    156167
    157168                /*relative criterion = norm(du)/norm(u)*/

  • issm/trunk/src/m/classes/@model/model.m

    r1629 r1666  
    167167
    168168        %Statics parameters
     169        md.eps_res=0;
    169170        md.eps_rel=0;
    170171        md.eps_abs=0;

  • issm/trunk/src/m/classes/@model/setdefaultparameters.m

    r1319 r1666  
    5656%Solver parameters
    5757
     58%mechanical residue convergence criterion norm(K(uold)uold - F)/norm(F)
     59md.eps_res=0.001;
     60
    5861%relative convergence criterion ((vel(n-1)-vel(n))/vel(n))
    5962md.eps_rel=0.01;

  • issm/trunk/src/m/classes/public/display/displaydiagnostic.m

    r1277 r1666  
    1212
    1313disp(sprintf('\n      %s','Newton convergence criteria:'));
    14 fielddisplay(md,'eps_rel','velocity relative convergence criterion');
    15 fielddisplay(md,'eps_abs','velocity absolute convergence criterion');
     14fielddisplay(md,'eps_res','mechanical equilibrium residue convergence criterion');
     15fielddisplay(md,'eps_rel','velocity relative convergence criterion, NaN -> not applied');
     16fielddisplay(md,'eps_abs','velocity absolute convergence criterion, NaN -> not applied');
    1617fielddisplay(md,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');
    1718

  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r1651 r1666  
    6565fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
    6666        'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
    67         'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
     67        'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    6868for i=1:length(fields),
    6969        if ~isempty(md.(fields{i})),
    … …  
    7575end
    7676
    77 %FIELDS > 0
     77%FIELDS >= 0
    7878fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
    79         'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
     79        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
    8080        'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    8181for i=1:length(fields),
    … …  
    9494%FIELDS ~=0
    9595fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
    96         'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',...
     96        'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
    9797        'sparsity','deltaH','DeltaH','timeacc','timedec'};
    9898for i=1:length(fields),

  • issm/trunk/src/m/classes/public/marshall.m

    r1650 r1666  
    119119WriteData(fid,md.mincontrolconstraint,'Scalar','mincontrolconstraint');
    120120WriteData(fid,md.maxcontrolconstraint,'Scalar','maxcontrolconstraint');
     121WriteData(fid,md.eps_res,'Scalar','eps_res');
    121122WriteData(fid,md.eps_rel,'Scalar','eps_rel');
    122123WriteData(fid,md.eps_abs,'Scalar','eps_abs');

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

    r1665 r1666  
    6060                %Figure out if convergence have been reached
    6161
    62                 %Residue criterion
     62                %Solver
    6363                solver_residue=norm(K_ff*soln(count).u_f-p_f,2)/norm(p_f,2);
    64                 displaystring(m.parameters.debug,'%-53s%g','      convergence criterion: norm(KU-F)/norm(F)',solver_residue);
     64                displaystring(m.parameters.debug>1,'%-60s%g','      convergence criterion: norm(K[n]U[n]-F)/norm(F)',solver_residue);
    6565
    6666                %Force equilibrium
    67                 mechanical_residue=norm(K_ff*soln(count-1).u_f-p_f,2)/norm(p_f,2);
    68                 displaystring(m.parameters.debug,'%-53s%g%s','      convergence criterion: norm(KUold-F)/norm(F)',mechanical_residue*100,' %');
     67                res=norm(K_ff*soln(count-1).u_f-p_f,2)/norm(p_f,2);
     68                if (res<=m.parameters.eps_res),
     69                        displaystring(m.parameters.debug,'%-60s%g%s%g%s','      mechanical equilibrium convergence criterion',res*100,' < ',m.parameters.eps_res*100,' %');
     70                        converged=1;
     71                else
     72                        displaystring(m.parameters.debug,'%-60s%g%s%g%s','      mechanical equilibrium convergence criterion',res*100,' > ',m.parameters.eps_res*100,' %');
     73                        converged=1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ATTTENETIONHiohfiow
     74                end
    6975
    7076                %Relative criterion
    … …  
    7379                nu=norm(soln(count-1).u_g,2);
    7480                if (ndu/nu<=m.parameters.eps_rel),
    75                         displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',m.parameters.eps_rel*100,' %');
     81                        displaystring(m.parameters.debug,'%-60s%g%s%g%s','      relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',m.parameters.eps_rel*100,' %');
    7682                        converged=1;
    7783                else
    78                         displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',m.parameters.eps_rel*100,' %');
     84                        displaystring(m.parameters.debug,'%-60s%g%s%g%s','      relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',m.parameters.eps_rel*100,' %');
    7985                        converged=0;
    8086                end
    … …  
    8490                if ~isnan(m.parameters.eps_abs),
    8591                        if (nduinf<=m.parameters.eps_abs),
    86                                 displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: max(du)',nduinf,' < ',m.parameters.eps_abs,' m/yr');
     92                                displaystring(m.parameters.debug,'%-60s%g%s%g%s','      absolute convergence criterion: max(du)',nduinf,' < ',m.parameters.eps_abs,' m/yr');
    8793                        else
    88                                 displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: max(du)',nduinf,' > ',m.parameters.eps_abs,' m/yr');
     94                                displaystring(m.parameters.debug,'%-60s%g%s%g%s','      absolute convergence criterion: max(du)',nduinf,' > ',m.parameters.eps_abs,' m/yr');
    8995                                converged=0;
    9096                        end
    9197                end
     98                disp(' ')
    9299
    93100                %rift convergence criterion
Note: See TracChangeset for help on using the changeset viewer.