Changeset 17629

Timestamp:

04/02/14 16:46:22 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing Newton algo with first guess

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)
• shared/Numerics/NewtonSolveDnorm.cpp (modified) (2 diffs)
• shared/Numerics/numerics.h (modified) (1 diff)

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

@@ -3552 +3552 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(Dstar);
-        xDelete<IssmDouble>(d);
+xDelete<IssmDouble>(d);
         xDelete<IssmDouble>(D);
         xDelete<IssmDouble>(tau);
@@ -4244 +4244 @@
         int         dim,tausize,meshtype;
         IssmDouble  epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
+        IssmDouble  epsxx_old,epsyy_old,epszz_old,epsxy_old,epsxz_old,epsyz_old;
         IssmDouble  sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
         IssmDouble  dvx[3],dvy[3],dvz[3],B,n;
@@ -4287 +4288 @@
                 }
 
+                /*Get previous d*/
+                Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
+                Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
+                Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
+                Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+                if(dim==3){
+                        epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
+                        epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
+                        epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+                }
+
                 /*Get tau*/
                 Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
@@ -4312 +4324 @@
                                 sigmapxz_input->GetInputValue(&sigmapxz,gauss);
                                 sigmapyz_input->GetInputValue(&sigmapyz,gauss);
+                        }
+
+                        /*Get previous d*/
+                        epsxx_input->GetInputValue(&epsxx_old,gauss);
+                        epsyy_input->GetInputValue(&epsyy_old,gauss);
+                        epsxy_input->GetInputValue(&epsxy_old,gauss);
+                        if(dim==3){
+                                epszz_input->GetInputValue(&epszz_old,gauss);
+                                epsxz_input->GetInputValue(&epsxz_old,gauss);
+                                epsyz_input->GetInputValue(&epsyz_old,gauss);
                         }
 
@@ -4333 +4355 @@
                         B_input->GetInputValue(&B,gauss);
                         n_input->GetInputValue(&n,gauss);
-                        coef1 = (B/2.)*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0
+                        coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2)  )^(n-1)/n ) 
                         coef2 = r;
                         if(dim==2){
@@ -4353 +4375 @@
                         }
                         IssmDouble dnorm;
-                        NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n);
+                        if(dim==2){
+                                dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old );
+                        }
+                        else{
+                                dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old 
+                                                        +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
+                        }
+                        NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n,dnorm);
 
                         /*Create Ke*/
@@ -4412 +4441 @@
                 IssmDouble* tau_xz = NULL;
                 IssmDouble* tau_yz = NULL;
-                Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
-                Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
-                Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
-                Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
-                if(dim==3){
-                        epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
-                        epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
-                        epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
-                }
                 delete gauss;
                 gauss = element->NewGauss();

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

@@ -3 +3 @@
 #include "../Exceptions/exceptions.h"
 
-int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n){
+int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm){
         /* solve the following equation using Newton-Raphson
          *
@@ -26 +26 @@
 
         /*Initial guess*/
-        IssmDouble y1 = 0.;
+        IssmDouble y1 = log10(dnorm);
 
         while(true){

issm/trunk-jpl/src/c/shared/Numerics/numerics.h

@@ -38 +38 @@
 int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
 
-int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n);
+int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm);
 
 #endif //ifndef _NUMERICS_H_