Changeset 1414

Timestamp:

07/29/09 15:29:34 (16 years ago)

Author:

Mathieu Morlighem

Message:

Added mechanical residue criterion

File:

• issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (modified) (5 diffs)

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

@@ -22 +22 @@
         Vec KU=NULL;
         Vec KUF=NULL;
+        Vec KUold=NULL;
+        Vec KUoldF=NULL;
 
         /*intermediary: */
@@ -38 +40 @@
         double ndu,nduinf,nu;
         double nKUF;
+        double nKUoldF;
         double nF;
-        double residue;
+        double solver_residue,mechanical_residue;
 
         /*parameters:*/
@@ -113 +116 @@
                 if (debug) _printf_("   solving\n");
                 Solverx(&uf, Kff, pf, old_uf, solver_string);
-
-                /*Compute residue*/
-                //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);
-                //comput norm(KUF), norm(F) and residue
-                VecNorm(KUF,NORM_2,&nKUF); VecFree(&KUF);
-                VecNorm(pf,NORM_2,&nF);
-                residue=nKUF/nF;
-        
-                //no need for Kff and pf anymore
-                MatFree(&Kff);VecFree(&pf);
                 
                 if (debug) _printf_("   merging solution from f to g set\n");
@@ -134 +125 @@
                 if (old_ug) inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes);
                 inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
-                
                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type); 
 
-                //Figure out if convergence is reached.
+                /*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); VecFree(&KUF);
+                VecNorm(pf,NORM_2,&nF);
+                solver_residue=nKUF/nF;
+                if (debug) _printf_("%-50s%g\n","   Convergence criterion: norm(KU-F)/norm(F)",solver_residue);
+
+                /*mechanical residue = norm(KUold-F)/norm(F)*/
+                if (count>1){
+                        //compute KUoldF = KUold - F = K*Uold - 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); VecFree(&KUoldF);
+                        mechanical_residue=nKUoldF/nF;
+                }
+                else{
+                        mechanical_residue=nKUoldF/nF;
+                }
+                if (debug) _printf_("%-50s%g%s\n","   Convergence criterion: norm(KUold-F)/norm(F)",mechanical_residue*100," \%");
+
+                /*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! "));
-
-                //residue
-                if (debug) _printf_("%s%g\n","   Convergence criterion: norm(KU-F)/norm(F)=",residue);
-
-                //relative criterion
                 if((ndu/nu)<eps_rel){
-                        if (debug) _printf_("%s%g%s%g\n","   Convergence criterion: norm(du)/norm(u)=",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_("%s%g%s%g\n","   Convergence criterion: norm(du)/norm(u)=",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=0;
                 }
 
-                //Absolute criterion (Optional)
+                /*Absolute criterion (Optional) = max(du)*/
                 if (!isnan(eps_abs)){
                         if ((nduinf*yts)<eps_abs){
-                                if (debug) _printf_("%s%g%s%g\n","   Convergence criterion: max(du)=",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_("%s%g%s%g\n","   Convergence criterion: max(du)=",nduinf*yts," > ",eps_abs);
+                                if (debug) _printf_("%-50s%g%s%g%s\n","   Convergence criterion: max(du)",nduinf*yts," > ",eps_abs," m/yr");
                                 converged=0;
                         }
@@ -172 +183 @@
                 count++;
                 if(converged==1)break;
-
+        
+                //no need for Kff and pf anymore
+                MatFree(&Kff);VecFree(&pf);
         }