Changeset 17553

Timestamp:

03/26/14 16:40:35 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with first version of XTH algorithm

Location:

issm/trunk-jpl/src

Files:

• c/analyses/StressbalanceAnalysis.cpp (modified) (6 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (2 diffs)
• c/shared/Numerics/NewtonSolveDnorm.cpp (modified) (3 diffs)
• c/solutionsequences/solutionsequence_la_theta.cpp (modified) (1 diff)
• m/enum/AugmentedLagrangianREnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)

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

@@ -2900 +2900 @@
         IssmDouble *xyz_list = NULL;
 
-        /*FIXME this should change*/
-        r = 1.;
-
         /*Get problem dimension*/
         element->FindParam(&meshtype,MeshTypeEnum);
+        element->FindParam(&r,AugmentedLagrangianREnum);
         switch(meshtype){
                 case Mesh2DverticalEnum: dim = 2; epssize = 3; break;
@@ -3289 +3287 @@
                 ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
                 ElementVector* pe4=CreatePVectorFSViscousXTH(element);
-                ElementVector* pe = new ElementVector(petemp,pe4);
+                pe = new ElementVector(petemp,pe4);
                 delete pe1;
                 delete pe2;
@@ -3403 +3401 @@
         }
 
-        /*FIXME*/
-        r = 1.;
-
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->NumberofNodesVelocity();
@@ -3427 +3422 @@
 
         /*Retrieve all inputs and parameters*/
+        element->FindParam(&r,AugmentedLagrangianREnum);
         element->GetVerticesCoordinates(&xyz_list);
 
@@ -4260 +4256 @@
         Gauss*      gauss = NULL;
 
+        parameters->FindParam(&r,AugmentedLagrangianREnum);
         parameters->FindParam(&meshtype,MeshTypeEnum);
         switch(meshtype){
@@ -4267 +4264 @@
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
-
-        /*FIXME*/
-        r = 1.;
 
         for(int i=0;i<elements->Size();i++){

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -624 +624 @@
         AbsoluteEnum,
         IncrementalEnum,
+        AugmentedLagrangianREnum,
         /*}}}*/
         /*Grounding Line{{{*/

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -603 +603 @@
                 case AbsoluteEnum : return "Absolute";
                 case IncrementalEnum : return "Incremental";
+                case AugmentedLagrangianREnum : return "AugmentedLagrangianR";
                 case AgressiveMigrationEnum : return "AgressiveMigration";
                 case NoneEnum : return "None";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -615 +615 @@
               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
+              else if (strcmp(name,"AugmentedLagrangianR")==0) return AugmentedLagrangianREnum;
               else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
               else if (strcmp(name,"None")==0) return NoneEnum;
@@ -628 +629 @@
               else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
               else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
-              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
+              if (strcmp(name,"Colinear")==0) return ColinearEnum;
+              else if (strcmp(name,"ControlSteady")==0) return ControlSteadyEnum;
               else if (strcmp(name,"Fset")==0) return FsetEnum;
               else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;

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

@@ -1 +1 @@
 #include <math.h>
 #include "./types.h"
+#include "../Exceptions/exceptions.h"
 
 int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n){
@@ -12 +13 @@
          */
 
+        /*trivial solution*/
+        if(c3==0.){
+                *pdnorm =0.;
+                return 0;
+        }
 
         /*Intermediaries*/
+        int        counter = 0;
         IssmDouble s = (1.+n)/n;
         IssmDouble y2;
@@ -24 +31 @@
 
                 /*Newton step*/
-                y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10)*pow(pow(10.,y1),s-1) + c2*log(10.)*pow(10.,y1));
+                y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10)*pow(pow(10.,y1),s-1.) + c2*log(10.)*pow(10.,y1));
 
-                if( fabs(y2-y1)/fabs(y1)<threshold ){
+                if( fabs(y2-y1)/(fabs(y2))<threshold ){
                         break;
                 }
                 else{
                         y1=y2;
+                        counter++;
                 }
+
+                if(counter>50) _error_("Could not converge");
         }
 
         /*Assign output pointer*/
         *pdnorm = pow(10.,y2);
+        return 0;
 }

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

@@ -30 +30 @@
         analysis->InitializeXTH(femmodel->elements,femmodel->parameters);
 
-        /*Solve KU=F*/
-        SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
-        CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-        Reduceloadx(pf, Kfs, ys); delete Kfs;
-        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
-        delete Kff; delete pf; delete df;
-        Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+        /*Convergence criterion*/
+        int count = 0;
+        femmodel->parameters->SetParam(10.,AugmentedLagrangianREnum);
 
-        /*Update solution*/
-        InputUpdateFromSolutionx(femmodel,ug); 
+        while(true){
 
-        /*Update d and tau accordingly*/
-        analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
+                /*Solve KU=F*/
+                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+                Reduceloadx(pf, Kfs, ys); delete Kfs;
+                Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
+                delete Kff; delete pf; delete df;
+                Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+
+                /*Update solution*/
+                printf("ug norm = %g\n",ug->Norm(NORM_TWO));
+                InputUpdateFromSolutionx(femmodel,ug); 
+
+                /*Update d and tau accordingly*/
+                analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
+
+                count++;
+                if(count>5) break;
+        }
 
         /*Check for convergence*/
-        _error_("STOP");
+        //_error_("STOP");
 
         delete ug;  

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -595 +595 @@
 def AbsoluteEnum(): return StringToEnum("Absolute")[0]
 def IncrementalEnum(): return StringToEnum("Incremental")[0]
+def AugmentedLagrangianREnum(): return StringToEnum("AugmentedLagrangianR")[0]
 def AgressiveMigrationEnum(): return StringToEnum("AgressiveMigration")[0]
 def NoneEnum(): return StringToEnum("None")[0]