Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1668) +++ /issm/trunk/src/c/Makefile.am (revision 1669) @@ -580,4 +580,5 @@ ./parallel/diagnostic_core_linear.cpp\ ./parallel/diagnostic_core_nonlinear.cpp\ + ./parallel/convergence.cpp\ ./parallel/thermal_core.cpp\ ./parallel/thermal_core_nonlinear.cpp\ Index: /issm/trunk/src/c/parallel/convergence.cpp =================================================================== --- /issm/trunk/src/c/parallel/convergence.cpp (revision 1669) +++ /issm/trunk/src/c/parallel/convergence.cpp (revision 1669) @@ -0,0 +1,138 @@ +/*!\file: convergence.cpp + * \brief: figure out if convergence has been reached + */ + +#undef __FUNCT__ +#define __FUNCT__ "convergence" + +#include "../objects/objects.h" +#include "../issm.h" + +int convergence(int* pconverged, Mat Kff,Vec pf,Vec uf,Vec old_uf,DataSet* parameters){ + + /*output*/ + int converged=0; + + /*intermediary*/ + Vec KU=NULL; + Vec KUF=NULL; + Vec KUold=NULL; + Vec KUoldF=NULL; + Vec duf=NULL; + double ndu,nduinf,nu; + double nKUF; + double nKUoldF; + double nF; + double solver_residue,res; + + /*convergence options*/ + double eps_res; + double eps_rel; + double eps_abs; + int debug; + double yts; + + /*get convergence options*/ + parameters->FindParam((void*)&eps_res,"eps_res"); + parameters->FindParam((void*)&eps_rel,"eps_rel"); + parameters->FindParam((void*)&eps_abs,"eps_abs"); + parameters->FindParam((void*)&yts,"yts"); + parameters->FindParam((void*)&debug,"debug"); + + /*Display solver caracteristics*/ + if (debug) _printf_("\n"); + 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; + _printf_("\n%s%g\n"," solver residue: norm(K[n]U[n]-F)/norm(F)",solver_residue); + + //clean up + VecFree(&KU); + VecFree(&KUF); + } + + /*Force equilibrium (Mandatory)*/ + + //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); + VecNorm(KUoldF,NORM_2,&nKUoldF); + VecNorm(pf,NORM_2,&nF); + res=nKUoldF/nF; + if (isnan(res)) throw ErrorException(__FUNCT__,exprintf("mechanical equilibrium convergence criterion is NaN! ")); + + //clean up + VecFree(&KUold); + VecFree(&KUoldF); + + //print + if(res ",eps_res*100," \%"); + converged=1;//TTTTTTTTTTTTTTTTTTFTYICGWFYUOGCYEUOGCYEUWGCHLWGWI + } + + /*Relative criterion (optional)*/ + if (!isnan(eps_rel) || (debug>1)){ + + //compute norm(du)/norm(u) + VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); + VecNorm(duf,NORM_2,&ndu); VecNorm(old_uf,NORM_2,&nu); + if (isnan(ndu) || isnan(nu)) throw ErrorException(__FUNCT__,exprintf("convergence criterion is NaN! ")); + + //clean up + VecFree(&duf); + + //print + if (!isnan(eps_rel)){ + if((ndu/nu) ",eps_rel*100," \%"); + converged=0; + } + } + else _printf_("%-50s%g%s\n"," Convergence criterion: norm(du)/norm(u)",ndu/nu*100," \%"); + + } + + /*Absolute criterion (Optional) = max(du)*/ + if (!isnan(eps_abs) || (debug>1)){ + + //compute max(du) + VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); + VecNorm(duf,NORM_2,&ndu); VecNorm(duf,NORM_INFINITY,&nduinf); + if (isnan(ndu) || isnan(nu)) throw ErrorException(__FUNCT__,exprintf("convergence criterion is NaN! ")); + + //clean up + VecFree(&duf); + + //print + if (!isnan(eps_rel)){ + if ((nduinf*yts) ",eps_abs," m/yr"); + converged=0; + } + } + else _printf_("%-50s%g%s\n"," Convergence criterion: max(du)",nduinf*yts," m/yr"); + + } + + /*assign output*/ + *pconverged=converged; +} Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1668) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1669) @@ -10,4 +10,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" #include "../issm.h" +#include "./parallel.h" void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ @@ -20,8 +21,4 @@ Vec old_uf=NULL; DataSet* loads=NULL; - Vec KU=NULL; - Vec KUF=NULL; - Vec KUold=NULL; - Vec KUoldF=NULL; /*intermediary: */ @@ -37,30 +34,16 @@ int numberofnodes; - Vec dug=NULL; - double ndu,nduinf,nu; - double nKUF; - double nKUoldF; - double nF; - double solver_residue,res; - /*parameters:*/ int kflag,pflag,connectivity,numberofdofspernode; char* solver_string=NULL; int debug=0; - 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. /*Recover parameters: */ kflag=1; pflag=1; - fem->parameters->FindParam((void*)&connectivity,"connectivity"); fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); 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"); - fem->parameters->FindParam((void*)&yts,"yts"); - fem->parameters->FindParam((void*)&debug,"debug"); /*Copy loads for backup: */ @@ -119,21 +102,6 @@ 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); - } - + //Merge back to g set if (debug) _printf_(" merging solution from f to g set\n"); - //Merge back to g set Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); @@ -145,59 +113,14 @@ /*Figure out if convergence is reached.*/ - - /*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); - res=nKUoldF/nF; - //clean up - VecFree(&KUold); - VecFree(&KUoldF); - 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)*/ - VecDuplicate(old_ug,&dug);VecCopy(old_ug,dug); VecAYPX(dug,-1.0,ug); - VecNorm(dug,NORM_2,&ndu);VecNorm(old_ug,NORM_2,&nu);VecNorm(dug,NORM_INFINITY,&nduinf); VecFree(&dug); - if (isnan(ndu) || isnan(nu)) throw ErrorException(__FUNCT__,exprintf("convergence criterion is NaN! ")); - if((ndu/nu) ",eps_rel*100," \%"); - converged=0; - } - - /*Absolute criterion (Optional) = max(du)*/ - if (!isnan(eps_abs)){ - if ((nduinf*yts) ",eps_abs," m/yr"); - converged=0; - } - } + convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); + MatFree(&Kff);VecFree(&pf); //rift convergence if (!constraints_converged) converged=0; - //no need for Kff and pf anymore - MatFree(&Kff);VecFree(&pf); - /*Increase count: */ count++; if(converged==1)break; - + } @@ -212,7 +135,5 @@ SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); - MatFree(&Kgg);VecFree(&pg); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 1668) +++ /issm/trunk/src/c/parallel/parallel.h (revision 1669) @@ -20,4 +20,5 @@ void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,DataSet* parameters); void transient_core(DataSet* results,FemModel* fems, ParameterInputs* inputs); Index: /issm/trunk/src/m/solutions/cielo/convergence.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/convergence.m (revision 1668) +++ /issm/trunk/src/m/solutions/cielo/convergence.m (revision 1669) @@ -16,5 +16,5 @@ solver_res=norm(K_ff*u_f-p_f,2)/norm(p_f,2); displaystring(debug>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); - displaystring(debug>1,'%s%g',' solver residue: norm(K[n]U[n]-F)/norm(F)',solver_res); + displaystring(debug>1,'%s%g',' solver residue norm(K[n]U[n]-F)/norm(F): ',solver_res); end @@ -33,19 +33,25 @@ %Relative criterion (optional) -if ~isnan(eps_rel), +if ((~isnan(eps_rel)) | (debug>1)), + %compute ndu/nu duf=u_f-u_f_old; ndu=norm(duf,2); - nu=norm(u_f,2); + nu=norm(u_f_old,2); if isnan(ndu), error('convergence error message: relative convergence criterion is NaN!'); end - if (ndu/nu<=eps_rel), - displaystring(debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',eps_rel*100,' %'); - converged=1; + %print criterion + if ~isnan(eps_rel), + if (ndu/nu<=eps_rel), + displaystring(debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',eps_rel*100,' %'); + converged=1; + else + displaystring(debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',eps_rel*100,' %'); + converged=0; + end else - displaystring(debug,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',eps_rel*100,' %'); - converged=0; + displaystring(debug,'%-60s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' %'); end @@ -53,6 +59,8 @@ %Absolute criterion (optional) -if ~isnan(eps_abs), +if ((~isnan(eps_abs)) | (debug>1)), + %compute max(du) + duf=u_f-u_f_old; nduinf=norm(duf,inf)*yts; if isnan(nduinf), @@ -60,9 +68,14 @@ end - if (nduinf<=eps_abs), - displaystring(debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' < ',eps_abs,' m/yr'); + %print criterion + if ~isnan(eps_abs), + if (nduinf<=eps_abs), + displaystring(debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' < ',eps_abs,' m/yr'); + else + displaystring(debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' > ',eps_abs,' m/yr'); + converged=0; + end else - displaystring(debug,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' > ',eps_abs,' m/yr'); - converged=0; + displaystring(debug,'%-60s%g%s',' absolute convergence criterion: max(du)',nduinf,' m/yr'); end