Changeset 17603

Timestamp:

03/31/14 13:33:46 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added convergence test for LA Theta algorithm

File:

• issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp

@@ -11 +11 @@
 
         /*intermediary: */
-        Matrix<IssmDouble>*  Kff = NULL;
-        Matrix<IssmDouble>*  Kfs = NULL;
-        Vector<IssmDouble>*  ug  = NULL;
-        Vector<IssmDouble>*  uf  = NULL;
-        Vector<IssmDouble>*  pf  = NULL;
-        Vector<IssmDouble>*  df  = NULL;
-        Vector<IssmDouble>*  ys  = NULL;
+        Matrix<IssmDouble>*  Kff    = NULL;
+        Matrix<IssmDouble>*  Kfs    = NULL;
+        Vector<IssmDouble>*  ug_old = NULL;
+        Vector<IssmDouble>*  ug     = NULL;
+        Vector<IssmDouble>*  uf     = NULL;
+        Vector<IssmDouble>*  pf     = NULL;
+        Vector<IssmDouble>*  df     = NULL;
+        Vector<IssmDouble>*  ys     = NULL;
         int  configuration_type;
 
@@ -31 +32 @@
 
         /*Convergence criterion*/
-        int count = 0;
-        femmodel->parameters->SetParam(reCast<IssmDouble>(10.),AugmentedLagrangianREnum);
+        int        count   = 0;
+        IssmDouble eps_rel = .001;
+        femmodel->parameters->SetParam(0.6,AugmentedLagrangianREnum);
+        GetSolutionFromInputsx(&ug,femmodel);
 
         while(true){
+
+                /*save pointer to old velocity*/
+                delete ug_old;ug_old=ug;
 
                 /*Solve KU=F*/
@@ -51 +57 @@
                 analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
 
+                /*Check for convergence*/
+                Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
+                IssmDouble ndu=dug->Norm(NORM_TWO);   delete dug;
+                IssmDouble nu =ug_old->Norm(NORM_TWO);
+                if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+                if((ndu/nu)<eps_rel){
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
+                        break;
+                }
+                else{ 
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+                }
+
                 count++;
-                if(count>5) break;
+                if(count>20) break;
         }
 
-        /*Check for convergence*/
-        //_error_("STOP");
-
         delete ug;  
+        delete ug_old;  
         delete analysis;
 }