Changeset 17650

Timestamp:

04/04/14 15:40:43 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: done with X-TH LA

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)
• shared/Numerics/NewtonSolveDnorm.cpp (modified) (3 diffs)
• solutionsequences/solutionsequence_la_theta.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -7 +7 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 10
+//#define FSANALYTICAL 2
 
 /*Model processing*/
@@ -3209 +3209 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
-        int         i,meshtype,dim;
+        int         i,meshtype,dim,fe_FS;
         IssmDouble  x_coord,y_coord,z_coord;
         IssmDouble  Jdet,forcex,forcey,forcez;
         IssmDouble *xyz_list = NULL;
 
+        element->FindParam(&fe_FS,FlowequationFeFSEnum);
         element->FindParam(&meshtype,MeshTypeEnum);
         switch(meshtype){
@@ -3271 +3272 @@
         xDelete<IssmDouble>(vbasis);
         xDelete<IssmDouble>(xyz_list);
+        if(fe_FS==XTaylorHoodEnum){
+                ElementVector* pe2=CreatePVectorFSViscousXTH(element);
+                ElementVector* pe3 = new ElementVector(pe,pe2);
+                delete pe;
+                delete pe2;
+                return pe3;
+        }
         return pe;
 }/*}}}*/
@@ -3552 +3560 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(Dstar);
-xDelete<IssmDouble>(d);
+        xDelete<IssmDouble>(d);
         xDelete<IssmDouble>(D);
         xDelete<IssmDouble>(tau);
@@ -4382 +4390 @@
                                                         +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
                         }
+                        /*Initial guess cannot be 0 otherwise log(0)  - inf*/
+                        if(dnorm==0.) dnorm=1.;
                         NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n,dnorm);
+                        _assert_(dnorm>=0.);
+                        _assert_(!xIsNan<IssmDouble>(dnorm));
 
                         /*Create Ke*/
@@ -4414 +4426 @@
                         Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
                         Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
+                        for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
+                        for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_yy[i]));
+                        for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
                         element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
                         element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
@@ -4420 +4435 @@
                 else{
                         _assert_(tnumnodes==4);
-                        Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
-                        Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
-                        Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
-                        Matrix3x3Solve(&d_zz[0],Ke,pe_zz);
-                        Matrix3x3Solve(&d_xz[0],Ke,pe_xz);
-                        Matrix3x3Solve(&d_yz[0],Ke,pe_yz);
+                        Matrix4x4Solve(&d_xx[0],Ke,pe_xx);
+                        Matrix4x4Solve(&d_yy[0],Ke,pe_yy);
+                        Matrix4x4Solve(&d_xy[0],Ke,pe_xy);
+                        Matrix4x4Solve(&d_zz[0],Ke,pe_zz);
+                        Matrix4x4Solve(&d_xz[0],Ke,pe_xz);
+                        Matrix4x4Solve(&d_yz[0],Ke,pe_yz);
                         element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
                         element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);

issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp

@@ -15 +15 @@
         /*trivial solution*/
         if(c3==0.){
-                *pdnorm =0.;
+                *pdnorm = 0.;
                 return 0;
         }
@@ -26 +26 @@
 
         /*Initial guess*/
+        _assert_(dnorm>0.); 
         IssmDouble y1 = log10(dnorm);
 
@@ -41 +42 @@
                 }
 
-                if(counter>50) _error_("Could not converge");
+                if(counter>50) break;
         }
+
+        /*Avoid extremely large values that indicate non convergence*/
+        if(y2>50.) y2 = 50;
 
         /*Assign output pointer*/

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

@@ -19 +19 @@
         Vector<IssmDouble>*  df     = NULL;
         Vector<IssmDouble>*  ys     = NULL;
+        IssmDouble           eps_rel;
         int  configuration_type,max_nonlinear_iterations;
 
@@ -26 +27 @@
         /*Recover parameters: */
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
         femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
+        femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
 
         /*Update constraints and initialize d and tau if necessary*/
@@ -33 +36 @@
 
         /*Convergence criterion*/
-        int        count   = 0;
-        IssmDouble eps_rel = .001;
-        femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
+        int  count = 0;
         GetSolutionFromInputsx(&ug,femmodel);
 
         while(true){
+                count++;
 
                 /*save pointer to old velocity*/
@@ -47 +49 @@
                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
                 Reduceloadx(pf, Kfs, ys); delete Kfs;
-
-        //      pf->Echo();
-        //      IssmDouble* temp=Kff->ToSerial();
-        //      int m,n;
-        //      Kff->GetSize(&m,&n);
-        //      printarray(temp,m,n);
-        //      xDelete<IssmDouble>(temp);
-
                 Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
-                //uf->Echo();
                 delete Kff; delete pf; delete df;
-
                 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
 
                 /*Update solution*/
-                _printf_("ug norm = " << ug->Norm(NORM_TWO) << "\n");
                 InputUpdateFromSolutionx(femmodel,ug); 
 
@@ -81 +72 @@
                 }
 
-                count++;
-                if(count>=max_nonlinear_iterations) break;
+                if(count>=max_nonlinear_iterations){
+                        _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+                        break;
+                }
         }
+
+        if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
 
         delete ug;