Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 1665) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 1666) @@ -70,4 +70,10 @@ parameters->AddObject(param); + /*eps_res: */ + count++; + param= new Param(count,"eps_res",DOUBLE); + param->SetDouble(model->eps_res); + parameters->AddObject(param); + /*eps_rel: */ count++; Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 1665) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 1666) @@ -128,4 +128,5 @@ model->debug=0; model->plot=0; + model->eps_res=0; model->eps_rel=0; model->eps_abs=0; @@ -346,4 +347,5 @@ ModelFetchData((void**)&model->mincontrolconstraint,NULL,NULL,model_handle,"mincontrolconstraint","Scalar",NULL); ModelFetchData((void**)&model->maxcontrolconstraint,NULL,NULL,model_handle,"maxcontrolconstraint","Scalar",NULL); + ModelFetchData((void**)&model->eps_res,NULL,NULL,model_handle,"eps_res","Scalar",NULL); ModelFetchData((void**)&model->eps_rel,NULL,NULL,model_handle,"eps_rel","Scalar",NULL); ModelFetchData((void**)&model->eps_abs,NULL,NULL,model_handle,"eps_abs","Scalar",NULL); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1665) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1666) @@ -126,4 +126,5 @@ int debug; int plot; + double eps_res; double eps_rel; double eps_abs; Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1665) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1666) @@ -42,5 +42,5 @@ double nKUoldF; double nF; - double solver_residue,mechanical_residue; + double solver_residue,res; /*parameters:*/ @@ -48,5 +48,5 @@ char* solver_string=NULL; int debug=0; - double eps_rel,eps_abs,yts; + double eps_res,eps_rel,eps_abs,yts; int dofs[4]={1,1,0,0}; //recover vx,vy by default. vz and pressure may be. @@ -58,4 +58,5 @@ fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); fem->parameters->FindParam((void*)&solver_string,"solverstring"); + fem->parameters->FindParam((void*)&eps_res,"eps_res"); fem->parameters->FindParam((void*)&eps_rel,"eps_rel"); fem->parameters->FindParam((void*)&eps_abs,"eps_abs"); @@ -117,4 +118,19 @@ if (debug) _printf_(" solving\n"); Solverx(&uf, Kff, pf, old_uf, solver_string); + + /*print solver residue = norm(KU-F)/norm(F)*/ + if (debug>1){ + //compute KUF = KU - F = K*U - F + VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU); + VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf); + //compute norm(KUF), norm(F) and residue + VecNorm(KUF,NORM_2,&nKUF); + VecNorm(pf,NORM_2,&nF); + solver_residue=nKUF/nF; + //clean up + VecFree(&KU); + VecFree(&KUF); + _printf_("%s%g\n"," solver residue: norm(K[n]U[n]-F)/norm(F)",solver_residue); + } if (debug) _printf_(" merging solution from f to g set\n"); @@ -130,28 +146,23 @@ /*Figure out if convergence is reached.*/ - /*solver residue = norm(KU-F)/norm(F)*/ - //compute KUF = KU - F = K*U - F - VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU); - VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf); - //compute norm(KUF), norm(F) and residue - VecNorm(KUF,NORM_2,&nKUF); - VecNorm(pf,NORM_2,&nF); - solver_residue=nKUF/nF; - //clean up - VecFree(&KU); - VecFree(&KUF); - if (debug) _printf_("%-50s%g\n"," Convergence criterion: norm(KU-F)/norm(F)",solver_residue); - - /*mechanical residue = norm(KUold-F)/norm(F)*/ - //compute KUoldF = KUold - F = K*Uold - F + /*1 mechanical equilibrium residue = norm(K[n]U[n-1]-F)/norm(F)*/ + //compute K[n]U[n-1]F = K[n]U[n-1] - F VecDuplicate(uf,&KUold); MatMultPatch(Kff,old_uf,KUold); VecDuplicate(KUold,&KUoldF);VecCopy(KUold,KUoldF); VecAYPX(KUoldF,-1.0,pf); //compute norm(KUF), norm(F) and residue VecNorm(KUoldF,NORM_2,&nKUoldF); - mechanical_residue=nKUoldF/nF; + res=nKUoldF/nF; //clean up VecFree(&KUold); VecFree(&KUoldF); - if (debug) _printf_("%-50s%g%s\n"," Convergence criterion: norm(KUold-F)/norm(F)",mechanical_residue*100," \%"); + if (isnan(res)) throw ErrorException(__FUNCT__,exprintf("mechanical equilibrium convergence criterion is NaN! ")); + if((res) ",eps_res*100," \%"); + converged=1;//TTTTTTTTTTTTTTTTTTFTYICGWFYUOGCYEUOGCYEUWGCHLWGWI + } /*relative criterion = norm(du)/norm(u)*/ Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 1665) +++ /issm/trunk/src/m/classes/@model/model.m (revision 1666) @@ -167,4 +167,5 @@ %Statics parameters + md.eps_res=0; md.eps_rel=0; md.eps_abs=0; Index: /issm/trunk/src/m/classes/@model/setdefaultparameters.m =================================================================== --- /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 1665) +++ /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 1666) @@ -56,4 +56,7 @@ %Solver parameters +%mechanical residue convergence criterion norm(K(uold)uold - F)/norm(F) +md.eps_res=0.001; + %relative convergence criterion ((vel(n-1)-vel(n))/vel(n)) md.eps_rel=0.01; Index: /issm/trunk/src/m/classes/public/display/displaydiagnostic.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 1665) +++ /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 1666) @@ -12,6 +12,7 @@ disp(sprintf('\n %s','Newton convergence criteria:')); -fielddisplay(md,'eps_rel','velocity relative convergence criterion'); -fielddisplay(md,'eps_abs','velocity absolute convergence criterion'); +fielddisplay(md,'eps_res','mechanical equilibrium residue convergence criterion'); +fielddisplay(md,'eps_rel','velocity relative convergence criterion, NaN -> not applied'); +fielddisplay(md,'eps_abs','velocity absolute convergence criterion, NaN -> not applied'); fielddisplay(md,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'); Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 1665) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 1666) @@ -65,5 +65,5 @@ fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',... - 'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; + 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; for i=1:length(fields), if ~isempty(md.(fields{i})), @@ -75,7 +75,7 @@ end -%FIELDS > 0 +%FIELDS >= 0 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',... - 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',... + 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',... 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'}; for i=1:length(fields), @@ -94,5 +94,5 @@ %FIELDS ~=0 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',... - 'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',... + 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',... 'sparsity','deltaH','DeltaH','timeacc','timedec'}; for i=1:length(fields), Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 1665) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 1666) @@ -119,4 +119,5 @@ WriteData(fid,md.mincontrolconstraint,'Scalar','mincontrolconstraint'); WriteData(fid,md.maxcontrolconstraint,'Scalar','maxcontrolconstraint'); +WriteData(fid,md.eps_res,'Scalar','eps_res'); WriteData(fid,md.eps_rel,'Scalar','eps_rel'); WriteData(fid,md.eps_abs,'Scalar','eps_abs'); Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1665) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1666) @@ -60,11 +60,17 @@ %Figure out if convergence have been reached - %Residue criterion + %Solver solver_residue=norm(K_ff*soln(count).u_f-p_f,2)/norm(p_f,2); - displaystring(m.parameters.debug,'%-53s%g',' convergence criterion: norm(KU-F)/norm(F)',solver_residue); + displaystring(m.parameters.debug>1,'%-60s%g',' convergence criterion: norm(K[n]U[n]-F)/norm(F)',solver_residue); %Force equilibrium - mechanical_residue=norm(K_ff*soln(count-1).u_f-p_f,2)/norm(p_f,2); - displaystring(m.parameters.debug,'%-53s%g%s',' convergence criterion: norm(KUold-F)/norm(F)',mechanical_residue*100,' %'); + res=norm(K_ff*soln(count-1).u_f-p_f,2)/norm(p_f,2); + if (res<=m.parameters.eps_res), + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' mechanical equilibrium convergence criterion',res*100,' < ',m.parameters.eps_res*100,' %'); + converged=1; + else + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' mechanical equilibrium convergence criterion',res*100,' > ',m.parameters.eps_res*100,' %'); + converged=1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ATTTENETIONHiohfiow + end %Relative criterion @@ -73,8 +79,8 @@ nu=norm(soln(count-1).u_g,2); if (ndu/nu<=m.parameters.eps_rel), - displaystring(m.parameters.debug,'%-53s%g%s%g%s',' convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',m.parameters.eps_rel*100,' %'); + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',m.parameters.eps_rel*100,' %'); converged=1; else - displaystring(m.parameters.debug,'%-53s%g%s%g%s',' convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',m.parameters.eps_rel*100,' %'); + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',m.parameters.eps_rel*100,' %'); converged=0; end @@ -84,10 +90,11 @@ if ~isnan(m.parameters.eps_abs), if (nduinf<=m.parameters.eps_abs), - displaystring(m.parameters.debug,'%-53s%g%s%g%s',' convergence criterion: max(du)',nduinf,' < ',m.parameters.eps_abs,' m/yr'); + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' < ',m.parameters.eps_abs,' m/yr'); else - displaystring(m.parameters.debug,'%-53s%g%s%g%s',' convergence criterion: max(du)',nduinf,' > ',m.parameters.eps_abs,' m/yr'); + displaystring(m.parameters.debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' > ',m.parameters.eps_abs,' m/yr'); converged=0; end end + disp(' ') %rift convergence criterion