CHG: basalstress.m in matlab, python supports Schoof-type friction law.

src/m/mech/basalstress.m

@@ -9,14 +9,10 @@
 %   See also: plot_basaldrag
 
 %Check md.friction class
-if ~isa(md.friction,'friction')
-	error('Error: md.friction only supports "friction.m" class.');
+if ~(isa(md.friction,'friction') | isa(md.friction,'frictionschoof'))
+	error('Error: md.friction only supports "friction.m", "frictionschoof.m" classes.');
 end
 
-%compute exponents
-s=averaging(md,1./md.friction.p,0);
-r=averaging(md,md.friction.q./md.friction.p,0);
-
 %Compute effective pressure
 g        =md.constants.g;
 rho_ice  =md.materials.rho_ice;
@@ -48,7 +44,27 @@
 uby=md.initialization.vy/md.constants.yts;
 
 %compute basal drag (S.I.)
-alpha2 = (N.^r).*(md.friction.coefficient.^2).*(ub.^(s-1));
+switch(class(md.friction))
+	case 'friction'
+		%compute exponents
+		s=averaging(md,1./md.friction.p,0);
+		r=averaging(md,md.friction.q./md.friction.p,0);
+
+		alpha2 = (N.^r).*(md.friction.coefficient.^2).*(ub.^(s-1));
+	case 'frictionschoof'
+		if any(N < 0)
+			%NOTE: Sometimes, negative value of effective pressure N gives image number in alpha2. To prevent the image value in alpha2, we use small values.
+			warning('Find effective pressure value < 0. Treating minimum effective value with 0.1');
+			N = max(N, 0.1);
+		end
+		m=averaging(md,md.friction.m,0);
+		C=averaging(md,md.friction.C,0);
+		Cmax=averaging(md,md.friction.Cmax,0);
+
+		alpha2 = (C.^2 .* ub.^(m-1))./(1 + (C.^2./(Cmax.*N)).^(1./m).*ub).^(m);
+	otherwise
+		error('not supported yet');
+end
 b  =  alpha2.*ub;
 bx = -alpha2.*ubx;
 by = -alpha2.*uby;

src/m/mech/basalstress.py

@@ -1,4 +1,5 @@
 from friction import friction
+from frictionschoof import frictionschoof
 from averaging import averaging
 import numpy as np
 
@@ -18,10 +19,6 @@ def basalstress(md):
     if not isinstance(md.friction,friction):
         raise Exception('Error: md.friction only supports "friction.m" class.')
 
-    #compute exponents
-    s=averaging(md,1/md.friction.p,0)
-    r=averaging(md,md.friction.q/md.friction.p,0)
-
     #Compute effective pressure
     g        =md.constants.g
     rho_ice  =md.materials.rho_ice
@@ -51,18 +48,34 @@ def basalstress(md):
     ub=np.sqrt(md.initialization.vx**2+md.initialization.vy**2)/md.constants.yts
     ubx=md.initialization.vx/md.constants.yts
     uby=md.initialization.vy/md.constants.yts
+    ub         =np.ravel(ub)
+    ubx        =np.ravel(ubx)
+    uby        =np.ravel(uby)
 
     #compute basal drag (S.I.)
     #reshape array to vector (N,)
-    coefficient=np.ravel(md.friction.coefficient)
     N          =np.ravel(N)
-    r          =np.ravel(r)
-    s          =np.ravel(s)
-    ub         =np.ravel(ub)
-    ubx        =np.ravel(ubx)
-    uby        =np.ravel(uby)
+
+    if isinstance(md.friction,friction):
+        #compute exponents
+        s=averaging(md,1/md.friction.p,0)
+        r=averaging(md,md.friction.q/md.friction.p,0)
+        coefficient=np.ravel(md.friction.coefficient)
+        r          =np.ravel(r)
+        s          =np.ravel(s)
+
+        alpha2 = (N**r)*(md.friction.coefficient**2)*(ub**(s-1))
+    elif isinstance(md.friction,frictionschoof):
+        if np.any(N<0):
+            %NOTE: Sometimes, N negative values gives image number in alpha2. To prevent the image value in alpha2, we use small values.
+        m=averaging(md,md.friction.m,0)
+        C=averaging(md,md.friction.C,0)
+        Cmax=averaging(md,md.friction.Cmax,0)
+
+        alpha2 = (C**2 * ub**(m-1))/(1 + (C**2/(Cmax*N))**(1/m)*ub)**m
+    else:
+        raise Exception('not supported yet')
 
-    alpha2 = (N**r)*(md.friction.coefficient**2)*(ub**(s-1))
     b  =  alpha2*ub
     bx = -alpha2*ubx
     by = -alpha2*uby