Changeset 25366

Timestamp:

08/10/20 15:46:03 (5 years ago)

Author:

Cheng Gong

Message:

NEW: add inversion for Weertman friction law. In order to achieve better convergence, Weertman friction law is modified to sigma_b=C^{2 |u_b|}(1/m-1)u_b

Location:

issm/trunk-jpl

Files:

• src/c/analyses/AdjointHorizAnalysis.cpp (modified) (2 diffs)
• src/c/classes/Loads/Friction.cpp (modified) (5 diffs)
• src/c/classes/Loads/Friction.h (modified) (1 diff)
• src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (modified) (1 diff)
• src/m/classes/frictionweertman.m (modified) (2 diffs)
• src/m/classes/frictionweertmantemp.m (modified) (2 diffs)
• test/Archives/Archive357.arch (added)
• test/NightlyRun/test357.m (added)

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

@@ -1713 +1713 @@
         basalelement->GradientIndexing(&vertexpidlist[0],control_index);
 
-        /* get the friction law: if 11-Schoof, use a special name for the coefficient*/
+        /* get the friction law: if 2-Weertman, 11-Schoof, use a special name for the coefficient*/
         element->FindParam(&frictionlaw, FrictionLawEnum);
         Input2* dragcoefficient_input;
         switch(frictionlaw) {
+                case 2:
                 case 11:
                         dragcoefficient_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoefficient_input);
@@ -1997 +1998 @@
                         dragcoeff_input = basalelement->GetInput2(FrictionCoefficientEnum); _assert_(dragcoeff_input);
                         break;
+                case 2:
                 case 11:
                         dragcoeff_input = basalelement->GetInput2(FrictionCEnum); _assert_(dragcoeff_input);

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

@@ -46 +46 @@
                         GetAlphaViscousComplement(palpha_complement,gauss);
                         break;
+                case 2:
+                        GetAlphaWeertmanComplement(palpha_complement, gauss);
+                        break;
                 case 3:
                         GetAlphaHydroComplement(palpha_complement,gauss);
@@ -134 +137 @@
 void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
 
-        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
+        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
          * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
-         * alpha_complement= Neff ^r * vel ^s*/
+         * alpha_complement= Neff ^r * vel ^(s-1)*/
 
         /*diverse: */
@@ -194 +197 @@
                 alpha_complement= pow(vmag, m-1.)*pow((1 + pow(C/(Cmax*Neff),1./m)*vmag), -m-1.);
         }
+        /*Assign output pointers:*/
+        *palpha_complement=alpha_complement;
+}/*}}}*/
+void Friction::GetAlphaWeertmanComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
+
+        /* Compute the complement of Weertman's law for inversion
+         * alpha2 = C^2 * vel^(1/m-1)
+         * alpha_complement = vel^(1/m-1)
+        */
+        /*diverse: */
+        IssmDouble  m;
+        IssmDouble  alpha_complement;
+
+        /*Recover parameters: */
+        element->GetInputValue(&m,gauss,FrictionMEnum);
+
+        /*Get effective pressure*/
+        IssmDouble vmag = VelMag(gauss);
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(vmag==0. && (1./m-1.)<0.) alpha_complement=0.;
+        else alpha_complement= pow(vmag, 1.0/m-1.);
+        
         /*Assign output pointers:*/
         *palpha_complement=alpha_complement;
@@ -530 +556 @@
 void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
-        /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
+        /*This routine calculates the basal friction coefficient alpha2= C^2 |v|^(1/m-1) */
 
         /*diverse: */
@@ -545 +571 @@
         /*Check to prevent dividing by zero if vmag==0*/
         if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
-        else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
+        else alpha2=C*C*pow(vmag,(1./m-1.));
 
         /*Assign output pointers:*/

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

@@ -29 +29 @@
                 void  GetAlphaViscousComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlphaSchoofComplement(IssmDouble* alpha_complement,Gauss* gauss);
+                void  GetAlphaWeertmanComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);

issm/trunk-jpl/src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp

@@ -67 +67 @@
         Input2* drag_input;
         switch(frictionlaw) {
+                case 2:
                 case 11:
                         drag_input = basalelement->GetInput2(FrictionCEnum); _assert_(drag_input);

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

@@ -9 +9 @@
                 m = NaN;
         end
+   methods (Static)
+      function self = loadobj(self) % {{{
+                        disp('Warning: the Weertman friciton law is updated to Sigma_b = C^2*|u_b|^(1/m-1)*u_b, since 2020-08-10');
+                end
+        end %}}}
         methods
                 function self = frictionweertman(varargin) % {{{
@@ -32 +37 @@
                         md = checkfield(md,'fieldname','friction.m','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
                 end % }}}
-                function disp(self) % {{{
+                function disp(self) % {{{z
                         disp('Weertman sliding law parameters:');
                         disp('   Weertman''s sliding law reads:');
-                        disp('      v_b = C * Sigma_b^m');
+                        disp('      v_b = C_w * Sigma_b^m');
                         disp('   In ISSM, this law is rewritten as:');
-                        disp('      Sigma_b = C^(-1/m) * |u_b|^(1/m-1)  u_b');
+                        disp('      Sigma_b = C^2 * |u_b|^(1/m-1) * u_b');
+                        disp('   where C_w=C^(-2m)');
                         disp(' ');
                         fielddisplay(self,'C','friction coefficient [SI]');

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

@@ -10 +10 @@
                 m = NaN;
         end
+   methods (Static)
+      function self = loadobj(self) % {{{
+                        disp('Warning: the Weertman friciton law is updated to Sigma_b = C^2*|u_b|^(1/m-1)*u_b, since 2020-08-10');
+                        disp(' and friction coefficient in the temprature dependent version is updated accordingly.');
+                end
+        end %}}}
         methods
                 function self = frictionweertmantemp(varargin) % {{{
@@ -31 +37 @@
                 function disp(self) % {{{
                         disp('Weertman sliding law parameters:');
-                        disp('      Sigma_b = C^(-1/m) * |u_b|^(1/m-1)  u_b * 1/f(T)');
+                        disp('      Sigma_b = C^2 * |u_b|^(1/m-1)  u_b * 1/f(T)');
                         disp(' ');
                         fielddisplay(self,'gamma','submelt sliding parameter f(T) = exp((T-Tpmp)/gamma)');