diff --git a/BasalDrag.m b/BasalDrag.m index dea9cdb0..8149021f 100755 --- a/BasalDrag.m +++ b/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 diff --git a/CalcBasalTraction.m b/CalcBasalTraction.m index 282c9c7a..0db87378 100755 --- a/CalcBasalTraction.m +++ b/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); diff --git a/GetSlipperyDistribution.m b/GetSlipperyDistribution.m index 92ed0fcd..0a85cf40 100755 --- a/GetSlipperyDistribution.m +++ b/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); diff --git a/TestSlipperinessInputValues.m b/TestSlipperinessInputValues.m index 0a1b2cdb..c8c2409e 100755 --- a/TestSlipperinessInputValues.m +++ b/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); diff --git a/Ua2D_DefaultParameters.m b/Ua2D_DefaultParameters.m index af2d3661..32c435b6 100755 --- a/Ua2D_DefaultParameters.m +++ b/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). diff --git a/UaFields.m b/UaFields.m index c9320435..1df8f996 100755 --- a/UaFields.m +++ b/UaFields.m @@ -53,6 +53,7 @@ q=[]; muk=[]; + V0=[] ; % This is a parameter in Joughin's sliding law, rCW-V0 Co=[]; mo=[] diff --git a/dIdCq.m b/dIdCq.m index 49a33bc5..ed00afab 100755 --- a/dIdCq.m +++ b/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 diff --git a/rCWN0.m b/rCWN0.m index 67cc77f6..62ee7c3f 100644 --- a/rCWN0.m +++ b/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; diff --git a/rCWV0.m b/rCWV0.m new file mode 100755 index 00000000..75b30a20 --- /dev/null +++ b/rCWV0.m @@ -0,0 +1,62 @@ + + + + +%function [Tauu,Tauv,dTauudu,dTauvdv,dTauudv,dTauvdu,dTauudh,dTauvdh] = rCWV0(C,C0,V0,h,hf,m,u,v,u0) +function [Tauu,Tauv,dTauudu,dTauvdv,dTauudv,dTauvdu,dTauudh,dTauvdh] = rCWV0(C,C0,V0,He,delta,m,u,v,u0) +%rCWV0 +% This function was generated by the Symbolic Math Toolbox version 23.2. +% 21-Nov-2023 16:30:23 + +t2 = V0.*m; +t3 = u.^2; +t4 = u0.^2; +t5 = v.^2; +t6 = C+C0; +t7 = -V0; +% t8 = -hf; +t9 = 1.0./m; +% t10 = h+t8; % h-hf +t11 = t2.*t4; +t12 = -t9; +t15 = t9./2.0; +t17 = t3+t4+t5; +t13 = delta; % t13 = dirac(t10); +t14 = He; % t14 = heaviside(t10); +t16 = t12-1.0; +t18 = t6.^t12; +t19 = t15-1.0./2.0; +t20 = t15-3.0./2.0; +t21 = sqrt(t17); +t22 = m.*t21; +t23 = V0+t21; +t27 = t17.^t19; +t28 = t17.^t20; +t24 = t4.*t22; +t25 = t2+t7+t22; +t26 = t23.^t12; +Tauu = t14.*t18.*t26.*t27.*u; +if nargout > 1 + Tauv = t14.*t18.*t26.*t27.*v; +end +if nargout > 2 + t29 = t23.^t16; + dTauudu = t9.*t14.*t18.*t28.*t29.*(t11+t24+V0.*t3+t2.*t5+t5.*t22); +end +if nargout > 3 + dTauvdv = t9.*t14.*t18.*t28.*t29.*(t11+t24+V0.*t5+t2.*t3+t3.*t22); +end +if nargout > 4 + t31 = t12.*t14.*t18.*t25.*t28.*t29.*u.*v; + dTauudv = t31; +end +if nargout > 5 + dTauvdu = t31; +end +if nargout > 6 + dTauudh = t13.*t18.*t26.*t27.*u; +end +if nargout > 7 + dTauvdh = t13.*t18.*t26.*t27.*v; +end +end diff --git a/uvMatrixAssemblySSTREAM.m b/uvMatrixAssemblySSTREAM.m index 9ef03c85..0a94c486 100644 --- a/uvMatrixAssemblySSTREAM.m +++ b/uvMatrixAssemblySSTREAM.m @@ -110,6 +110,11 @@ qnod=reshape(F.q(MUA.connectivity,1),MUA.Nele,MUA.nod); end +if ~isempty(F.V0) + V0nod=reshape(F.V0(MUA.connectivity,1),MUA.Nele,MUA.nod); +end + + if ~isempty(F.muk) muknod=reshape(F.muk(MUA.connectivity,1),MUA.Nele,MUA.nod); end @@ -195,7 +200,7 @@ end - + Cint=Cnod*fun; Cint(Cint