Changeset 1669

Timestamp:

08/13/09 14:50:04 (15 years ago)

Author:

Mathieu Morlighem

Message:

separated convergence from diagnstic_core_nonlinear

Location:

issm/trunk/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/parallel/convergence.cpp (added)
• c/parallel/diagnostic_core_nonlinear.cpp (modified) (6 diffs)
• c/parallel/parallel.h (modified) (1 diff)
• m/solutions/cielo/convergence.m (modified) (4 diffs)

issm/trunk/src/c/Makefile.am

@@ -580 +580 @@
                                         ./parallel/diagnostic_core_linear.cpp\
                                         ./parallel/diagnostic_core_nonlinear.cpp\
+                                        ./parallel/convergence.cpp\
                                         ./parallel/thermal_core.cpp\
                                         ./parallel/thermal_core_nonlinear.cpp\

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

@@ -10 +10 @@
 #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 +21 @@
         Vec old_uf=NULL; 
         DataSet* loads=NULL;
-        Vec KU=NULL;
-        Vec KUF=NULL;
-        Vec KUold=NULL;
-        Vec KUoldF=NULL;
 
         /*intermediary: */
@@ -37 +34 @@
         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 +102 @@
                 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 +113 @@
 
                 /*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){
-                        if (debug) _printf_("%-50s%g%s%g%s\n","   mechanical equilibrium convergence criterion",res*100," < ",eps_res*100," \%");
-                        converged=1;
-                }
-                else{ 
-                        if (debug) _printf_("%-50s%g%s%g%s\n","   mechanical equilibrium convergence criterion",res*100," > ",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){
-                        if (debug) _printf_("%-50s%g%s%g%s\n","   Convergence criterion: norm(du)/norm(u)",ndu/nu*100," < ",eps_rel*100," \%");
-                        converged=1;
-                }
-                else{ 
-                        if (debug) _printf_("%-50s%g%s%g%s\n","   Convergence criterion: norm(du)/norm(u)",ndu/nu*100," > ",eps_rel*100," \%");
-                        converged=0;
-                }
-
-                /*Absolute criterion (Optional) = max(du)*/
-                if (!isnan(eps_abs)){
-                        if ((nduinf*yts)<eps_abs){
-                                if (debug) _printf_("%-50s%g%s%g%s\n","   Convergence criterion: max(du)",nduinf*yts," < ",eps_abs," m/yr");
-                        }
-                        else{
-                                if (debug) _printf_("%-50s%g%s%g%s\n","   Convergence criterion: max(du)",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 +135 @@
         
                 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);
 

issm/trunk/src/c/parallel/parallel.h

@@ -20 +20 @@
 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);

issm/trunk/src/m/solutions/cielo/convergence.m

@@ -16 +16 @@
         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 +33 @@
 
 %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 +59 @@
 
 %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 +68 @@
         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