Changeset 23840

Timestamp:

04/11/19 14:40:05 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed Schoof

Location:

issm/trunk-jpl/src

Files:

• c/analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• c/classes/Loads/Friction.cpp (modified) (14 diffs)
• m/classes/frictionschoof.m (modified) (4 diffs)

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

@@ -871 +871 @@
                         break;
                 case 11:
+                        iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
                         iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
                         iomodel->FetchDataToInput(elements,"md.friction.Cmax",FrictionCmaxEnum);
+                        break;
                 default:
                         _error_("friction law "<< frictionlaw <<" not supported");
@@ -968 +970 @@
                         break;
                 case 11:
-                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.m",FrictionMEnum));
+                        parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
                         break;
                 default: _error_("Friction law "<<frictionlaw<<" not implemented yet");

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -56 +56 @@
         }
 
+
+        /*Checks*/
+        _assert_(!xIsNan<IssmDouble>(*palpha_complement));
+        _assert_(!xIsInf<IssmDouble>(*palpha_complement));
+
 }/*}}}*/
 void Friction::GetAlphaHydroComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
@@ -93 +98 @@
                                         (Gamma/As - (alpha*q_exp*pow(Chi,q_exp-1.)* Gamma * Gamma/As));
 
-        _assert_(!xIsNan<IssmDouble>(alpha_complement));
         /*Assign output pointers:*/
         *palpha_complement=alpha_complement;
@@ -154 +158 @@
         else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
 
-        _assert_(!xIsNan<IssmDouble>(alpha_complement));
-        _assert_(!xIsInf<IssmDouble>(alpha_complement));
-
         /*Assign output pointers:*/
         *palpha_complement=alpha_complement;
@@ -200 +201 @@
                         _error_("Friction law "<< this->law <<" not supported");
         }
+
+        /*Checks*/
+        _assert_(!xIsNan<IssmDouble>(*palpha2));
+        _assert_(!xIsInf<IssmDouble>(*palpha2));
+        _assert_(*palpha2>=0);
 
 }/*}}}*/
@@ -245 +251 @@
         if(alpha2_coulomb<alpha2) alpha2=alpha2_coulomb;
 
-        _assert_(!xIsNan<IssmDouble>(alpha2));
-
         /*Assign output pointers:*/
         *palpha2=alpha2;
@@ -292 +296 @@
         else    alpha2=Neff * C_param * pow(Gamma,1./n) * pow(vmag,-1);
 
-        _assert_(!xIsNan<IssmDouble>(alpha2));
         /*Assign output pointers:*/
         *palpha2=alpha2;
@@ -324 +327 @@
 
         alpha2=drag_coefficient*drag_coefficient*Neff;
-        _assert_(!xIsNan<IssmDouble>(alpha2));
 
         /*Assign output pointers:*/
@@ -396 +398 @@
         alpha2 = alpha2 * exp((deltaTref - deltaT)/(2*gamma));
 
-        /*Final checks in debuging mode*/
-        _assert_(!xIsNan<IssmDouble>(alpha2));
-        _assert_(!xIsInf<IssmDouble>(alpha2));
-
         /*Assign output pointers:*/
         *palpha2=alpha2;
@@ -430 +428 @@
         if(vmag==0. && (s-1.)<0.) alpha2=0.;
         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
-        _assert_(!xIsNan<IssmDouble>(alpha2));
 
         /*Assign output pointers:*/
@@ -477 +474 @@
         if(vmag==0. && (s-1.)<0.) alpha2=0.;
         else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
-        _assert_(!xIsNan<IssmDouble>(alpha2));
 
         /*Assign output pointers:*/
@@ -500 +496 @@
         if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
         else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
-        _assert_(!xIsNan<IssmDouble>(alpha2));
 
         /*Assign output pointers:*/
@@ -595 +590 @@
         IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q));
 
-        /*Final checks in debuging mode*/
-        _assert_(!xIsNan<IssmDouble>(alpha2));
-        _assert_(!xIsInf<IssmDouble>(alpha2));
-
         /*Assign output pointers:*/
         *palpha2=alpha2;
@@ -613 +604 @@
 
         /*diverse: */
-        IssmDouble  C,Cmax,m;
+        IssmDouble  C,Cmax,m,alpha2;
 
         /*Recover parameters: */
@@ -625 +616 @@
 
         /*Compute alpha^2*/
-        IssmDouble alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
-
-        /*Checks*/
-        _assert_(!xIsNan<IssmDouble>(alpha2));
-        _assert_(alpha2>=0);
+        if(ub<1e-10){
+                alpha2 = 0.;
+        }
+        else{
+                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
+        }
 
         /*Assign output pointers:*/

issm/trunk-jpl/src/m/classes/frictionschoof.m

@@ -8 +8 @@
                 C    = NaN;
                 Cmax = NaN;
-                m    = 0.;
+                m    = NaN;
         end
         methods
@@ -25 +25 @@
                 function self = setdefaultparameters(self) % {{{
 
-                        %default m is 1/n = 1/3
-                        self.m = 1./3.;
-
                 end % }}}
                 function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -35 +32 @@
                         md = checkfield(md,'fieldname','friction.C','timeseries',1,'NaN',1,'Inf',1,'>',0.);
                         md = checkfield(md,'fieldname','friction.Cmax','timeseries',1,'NaN',1,'Inf',1,'>',0.);
-                        md = checkfield(md,'fieldname','friction.m','NaN',1,'Inf',1,'numel',1,'>',0.);
+                        md = checkfield(md,'fieldname','friction.m','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofelements,1]);
                 end % }}}
                 function disp(self) % {{{
@@ -55 +52 @@
                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','C','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','Cmax','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
-                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','m','format','Double');
+                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','m','format','DoubleMat','mattype',2);