Joughin sliding law being added, not finalized

BasalDrag.m

@@ -1,7 +1,15 @@
+
+
+
+
+
+
+
+
 function [taubx,tauby,dtaubxdu,dtaubxdv,dtaubydu,dtaubydv,dtaubxdh,dtaubydh,taubxo,taubyo,taubxa,taubya] = ...
-        BasalDrag(CtrlVar,MUA,He,delta,h,B,H,rho,rhow,ub,vb,C,m,uo,vo,Co,mo,ua,va,Ca,ma,q,g,muk)
+        BasalDrag(CtrlVar,MUA,He,delta,h,B,H,rho,rhow,ub,vb,C,m,uo,vo,Co,mo,ua,va,Ca,ma,q,g,muk,V0)
 
-    narginchk(24,24)
+    narginchk(25,25)
 
     %%
     % Returns basal drag and the directional derivatives of basal drag with respect to u,
@@ -119,13 +127,13 @@
 
                 U=speed;
                 N=N0(CtrlVar,h,H,rho,rhow,g) ;
-                dFuvdC=(U.^(1.0./m-1.0).*muk.^(m+1.0).*N.^(m+1.0).*He.*(muk.^m.*N.^m+U.*(He.*(C+C0).^(-1.0./m)).^m).^(-1.0./m-1.0).*(C+C0).^(-1.0./m-1.0))./m;
+                dFuvdC=(U.^(1./m-1).*muk.^(m+1).*N.^(m+1).*He.*(muk.^m.*N.^m+U.*(He.*(C+C0).^(-1./m)).^m).^(-1./m-1).*(C+C0).^(-1./m-1))./m;
 
             case {"rCW-N0","Umbi"} % reciprocal Coulumb-Weertman with zeroth-order hydrology
 
                 U=speed;
                 N=N0(CtrlVar,h,H,rho,rhow,g) ;
-                dFuvdC=(U.^(1.0./m).*muk.^2.*N.^2.*He.*1.0./(U.^(1.0./m).*He+muk.*N.*(C+C0).^(1.0./m)).^2.*(C+C0).^(1.0./m-1.0))./(U.*m) ;
+                dFuvdC=(U.^(1./m).*muk.^2.*N.^2.*He.*1./(U.^(1./m).*He+muk.*N.*(C+C0).^(1./m)).^2.*(C+C0).^(1./m-1))./(U.*m) ;
 
             case {"Tsai","minCW-N0"}
 
@@ -147,6 +155,12 @@
 
                 fprintf("Inversion using Coulomb sliding law not implemented. \n")
                 error("BasalDrag:InvalidCase","Inversion using Tsai sliding law not implemented.")
+
+            case {"Joughan","rCW-v0"}
+
+                 U=speed;
+                 dFuvdC=   (U.^(1./m-1).*He.*(U+V0).^(-1./m).*(C+C0).^(-1./m-1))./m ;
+
 
             otherwise
 
@@ -210,21 +224,28 @@
             dtaubydui=dtaubxdvi;  % just symmetry, always true, both Weertman and Coulom
             dtaubxdhi(isCoulomb)=dtaubxdhiC(isCoulomb);
             dtaubydhi(isCoulomb)=dtaubydhiC(isCoulomb);
-            
-            
+
+
         case {"rpCW-N0","Cornford"}
-            
+
             [N,dNdh]=N0(CtrlVar,h,H,rho,rhow,g) ;
             [taubxi,taubyi,dtaubxdui,dtaubydvi,dtaubxdvi,dtaubydui,dtaubxdhi,dtaubydhi] =  rpCWN0(C,CtrlVar.Czero,N,dNdh,He,delta,m,muk,ub,vb,CtrlVar.SpeedZero) ;
-            
+
         case {"rCW-N0","Umbi"} % reciprocal Coulumb-Weertman with zeroth-order hydrology
-            
+
             [N,dNdh]=N0(CtrlVar,h,H,rho,rhow,g) ;
             [taubxi,taubyi,dtaubxdui,dtaubydvi,dtaubxdvi,dtaubydui,dtaubxdhi,dtaubydhi] =  rCWN0(C,CtrlVar.Czero,N,dNdh,He,delta,m,muk,ub,vb,CtrlVar.SpeedZero) ;
-
-
+
+        case {"Joughin","rCW-v0"}
+
+            C0=CtrlVar.Czero;
+            u0=CtrlVar.SpeedZero;
+
+            [taubxi,taubyi,dtaubxdui,dtaubydvi,dtaubxdvi,dtaubydui,dtaubxdhi,dtaubydhi] = rCWV0(C,C0,V0,He,delta,m,ub,vb,u0)   ;
+
+
         otherwise
-            
+
             error("BasalDrag:CaseNotFound","what sliding law?")
     end
 

CalcBasalTraction.m

@@ -47,7 +47,7 @@
 delta = DiracDelta(CtrlVar.kH,F.h-hf,CtrlVar.Hh0) ;
 
 [tbx,tby] = ...
-    BasalDrag(CtrlVar,MUA,He,delta,F.h,F.B,F.S-F.B,F.rho,F.rhow,F.ub,F.vb,F.C,F.m,F.uo,F.vo,F.Co,F.mo,F.ua,F.va,F.Ca,F.ma,F.q,F.g,F.muk);
+    BasalDrag(CtrlVar,MUA,He,delta,F.h,F.B,F.S-F.B,F.rho,F.rhow,F.ub,F.vb,F.C,F.m,F.uo,F.vo,F.Co,F.mo,F.ua,F.va,F.Ca,F.ma,F.q,F.g,F.muk,F.V0);
 
 tb=sqrt(tbx.^2+tby.^2);
 

GetSlipperyDistribution.m

@@ -42,6 +42,7 @@
         else
             [UserVar,F.C,F.m,F.q]=DefineSlipperyDistribution(UserVar,CtrlVar,MUA,F);
         end
+
     case 5
 
         if NargInputFile>4
@@ -53,7 +54,19 @@
             [UserVar,F.C,F.m,F.q,F.muk]=DefineSlipperyDistribution(UserVar,CtrlVar,MUA,F);
 
         end
-
+
+    case 6
+
+        if NargInputFile>4
+
+            [UserVar,F.C,F.m,F.q,F.muk,F.V0]=DefineSlipperyDistribution(UserVar,CtrlVar,MUA,CtrlVar.time,F.s,F.b,F.h,F.S,F.B,F.rho,F.rhow,F.GF);
+
+        else
+
+            [UserVar,F.C,F.m,F.q,F.muk.F.V0]=DefineSlipperyDistribution(UserVar,CtrlVar,MUA,F);
+
+        end
+
     otherwise
 
         error('Ua:GetSlipperyDistribution','DefineSlipperyDistribution must return between 3 and 5 arguments')
@@ -63,7 +76,7 @@
 F.Cmax=CtrlVar.Cmax;
 F.Cmin=CtrlVar.Cmin;
 
-[F.C,F.m,F.q,F.muk]=TestSlipperinessInputValues(CtrlVar,MUA,F.C,F.m,F.q,F.muk);
+[F.C,F.m,F.q,F.muk,F.V0]=TestSlipperinessInputValues(CtrlVar,MUA,F.C,F.m,F.q,F.muk,F.V0);
 
 
 

TestSlipperinessInputValues.m

@@ -1,8 +1,8 @@
 
-function [C,m,q,muk]=TestSlipperinessInputValues(CtrlVar,MUA,C,m,q,muk)
+function [C,m,q,muk,V0]=TestSlipperinessInputValues(CtrlVar,MUA,C,m,q,muk,V0)
 
-narginchk(4,6)
-nargoutchk(2,4)
+narginchk(4,7)
+nargoutchk(2,5)
 
 
 
@@ -18,6 +18,10 @@
     muk=[];
 end
 
+if nargin<7
+    V0=[];
+end
+
 
 if numel(C)==1
     %fprintf(' C given by user is a scalar. Assuming that C is same everywhere. \n')
@@ -30,7 +34,7 @@
 
 
 if numel(m)==1
-    
+
     if  CtrlVar.CisElementBased
         m=m+zeros(MUA.Nele,1);
     else
@@ -41,90 +45,119 @@
 
 
 if isempty(q)
-    
+
     pattern=["Budd","W-N0"];
     if contains(CtrlVar.SlidingLaw,pattern)
         fprintf("Input Error:  \t For sliding law: %s \n \t \t \t \t q must be defined in DefineSlipperiness.m \n",CtrlVar.SlidingLaw)
         fprintf("\t \t \t \t and in an inverse run in DefineInputsForInverseRun.m as well. \n")
         error("Incorrect inputs")
-        
+
     end
-    
+
 else
-    
+
     if numel(q)==1
-        
+
         if  CtrlVar.CisElementBased
             q=q+zeros(MUA.Nele,1);
         else
             q=q+zeros(MUA.Nnodes,1);
         end
     end
-    
+
 end
 
 
 
 if isempty(muk)
-    
+
     pattern=["Tsai","Coulomb","Cornford","Umbi","minCW-N0","rpCW-N0","rCW-N0"]  ;
     if contains(CtrlVar.SlidingLaw,pattern)
         fprintf("Input Error:  \t For sliding law: %s \n \t \t \t \t muk must be defined in DefineSlipperiness.m \n",CtrlVar.SlidingLaw)
         fprintf("\t \t \t \t and in an inverse run in DefineInputsForInverseRun.m as well. \n")
         error("Incorrect inputs")
-        
+
     end
-    
+
 else
-    
+
     if numel(muk)==1
-        
+
         if  CtrlVar.CisElementBased
             muk=muk+zeros(MUA.Nele,1);
         else
             muk=muk+zeros(MUA.Nnodes,1);
         end
     end
-
+
+end
+
+pattern=["Joughin","rCW-V0"];
+if ~contains(CtrlVar.SlidingLaw,pattern)
+    V0=[];
 end
 
+if isempty(V0)
+
+    pattern=["Joughin","rCW-V0"];
+    if contains(CtrlVar.SlidingLaw,pattern)
+        fprintf("Input Error:  \t For sliding law: %s \n \t \t \t \t V0 must be defined in DefineSlipperiness.m \n",CtrlVar.SlidingLaw)
+        fprintf("\t \t \t \t and in an inverse run in DefineInputsForInverseRun.m as well. \n")
+        error("Incorrect inputs")
+
+    end
+
+else
+
+    if numel(V0)==1
+
+        if  CtrlVar.CisElementBased
+            V0=V0+zeros(MUA.Nele,1);
+        else
+            V0=V0+zeros(MUA.Nnodes,1);
+        end
+    end
+
+end
+
+
 
 
 if CtrlVar.AutomaticallyMapAGlenBetweenNodesAndEleIfEnteredIncorrectly
-    
+
     if CtrlVar.CisElementBased
-        
+
         if nC==MUA.Nnodes
             fprintf('\n Note:  C is element based, but entered on input as a nodal variable.\n')
             fprintf('        C will be mapped from nodes to elements by averaging over neighbouring nodes. \n')
             C=Nodes2EleMean(MUA.connectivity,C);
         end
-        
+
         if nm==MUA.Nnodes
             fprintf('\n Note:  m is element based, but entered on input as a nodal variable.\n')
             fprintf('        m will be mapped from nodes to elements by averaging over neighbouring nodes. \n')
             m=Nodes2EleMean(MUA.connectivity,m);
         end
-        
+
     else
-        
+
         if nC==MUA.Nele || nm==MUA.Nele
-            
+
             M=Ele2Nodes(MUA.connectivity,MUA.Nnodes);
-            
+
             if nC==MUA.Nele
-                
+
                 fprintf('\n Note:  C is nodal based, but entered on input as an element variable.\n')
                 fprintf('        C will be mapped from elements to nodes by averaging over neighbouring elements. \n')
-                
+
                 C=M*C;
             end
-            
+
             if nm==MUA.Nele
-                
+
                 fprintf('\n Note:  m is nodal based, but entered on input as an element variable.\n')
                 fprintf('        m will be mapped from elements to nodes by averaging over neighbouring elements. \n')
-                
+
                 m=M*m;
             end
         end
@@ -140,7 +173,7 @@
 elseif ~CtrlVar.CisElementBased && ~(length(MUA.coordinates) == length(C))
     save TestSave ;
     error(' C is node-based but input does not have same number of elements as there are nodes in mesh. All variables saved in TestSave.mat ')
-    
+
 end
 
 
@@ -162,7 +195,7 @@
     if niU>0
         fprintf('TestCInputValues: On input %i C values are larger than largest allowed value (%g).  \n',niU);
         fprintf(' These values have been reset to the maximum allowed value of %g .\n',niU,CtrlVar.Cmax);
-        
+
     end
     if niL>0
         fprintf('TestCInputValues: On input %i C values are smaller than smallest allowed value (%g). \n',niL,CtrlVar.Cmin);

Ua2D_DefaultParameters.m

@@ -85,7 +85,7 @@
 % model, where N=rho g (h-h_f) where h_f is the flotation thickness. 
 %
 CtrlVar.SlidingLaw="Weertman" ;
-CtrlVar.MustBe.SlidingLaw=["Weertman","Budd","Tsai","Coulomb","Cornford","Umbi","W","W-N0","minCW-N0","C","rpCW-N0","rCW-N0"]  ;
+CtrlVar.MustBe.SlidingLaw=["Weertman","Budd","Tsai","Coulomb","Cornford","Umbi","Joughin","W","W-N0","minCW-N0","C","rpCW-N0","rCW-N0","rCW-V0"]  ;
 %% Boundary conditions
 CtrlVar.UpdateBoundaryConditionsAtEachTimeStep=0;  % if true, `DefineBoundaryConditions.m' is called at the beginning of each time step to update the boundary conditions.
                                                    % otherwise boundary conditions are only updated at the beginning of the run (also at the beginning or a restart run).

UaFields.m

@@ -53,6 +53,7 @@
 
         q=[]; 
         muk=[];
+        V0=[] ; % This is a parameter in Joughin's sliding law, rCW-V0
 
         Co=[];
         mo=[]

dIdCq.m

@@ -48,6 +48,12 @@
     muknod=mnod*0 ;
 end
 
+if ~isempty(F.V0)
+    V0nod=reshape(F.V0(MUA.connectivity,1),MUA.Nele,MUA.nod);
+else
+    V0nod=mnod*0 ;
+end
+
 
 Bnod=reshape(F.B(MUA.connectivity,1),MUA.Nele,MUA.nod);
 Snod=reshape(F.S(MUA.connectivity,1),MUA.Nele,MUA.nod);
@@ -74,6 +80,7 @@
     mint=mnod*fun;
     qint=qnod*fun;
     mukint=muknod*fun;
+    V0int=V0nod*fun;
     Bint=Bnod*fun;
     Sint=Snod*fun;
     Hint=Sint-Bint;
@@ -94,7 +101,7 @@
     % setting this CtrlVar field to true ensures that BasalDrag.m returns the (point) derivative
     CtrlVar.Inverse.dFuvdClambda=true;
     Ctemp= ...
-        BasalDrag(CtrlVar,MUA,Heint,[],hint,Bint,Hint,rhoint,F.rhow,uint,vint,Cint,mint,[],[],[],[],[],[],[],[],qint,F.g,mukint);
+        BasalDrag(CtrlVar,MUA,Heint,[],hint,Bint,Hint,rhoint,F.rhow,uint,vint,Cint,mint,[],[],[],[],[],[],[],[],qint,F.g,mukint,V0int);
     CtrlVar.Inverse.dFuvdClambda=false;
 
     if ~CtrlVar.DevelopmentVersion

rCWN0.m

@@ -26,7 +26,7 @@
 t30 = t11.*t13.*(-1.0./2.0);
 t15=delta ; % t15 = dirac(t14);
 t16=He ; % t16 = heaviside(t14);
-% t19=0 ; % t19 = (t14 == 0.0);  ?! for some reason the simplificatoin gives a kronecka delta function
+% t19=0 ; % t19 = (t14 == 0.0);  ?! for some reason the simplification gives a kronecka delta function
 t22 = mu.*t2.*t17;
 %t23 = t11.*t18;
 %t24 = t12.*t18;