In case this helps, here is some more info, assuming that you are using a weertman friction law and want to try out different exponents m
:
1. Assuming you are happy with a run with m=1
, you can run a stress balance:
%make sure we use linearize = 2 (for element-wise friction coefficient)
md.friction.linearize = 2;
%Set friction to m=1
m = 1;
md.friction.m = m*ones(md.mesh.numberofelements,1); %m = 1;
md.verbose.convergence = 1;
md = solve(md,'sb');
2. Get basal stress and basal speed so that we can rescale the friction coefficient for other values of m
%Get basal sliding speed and basal stress from run with m=1 (all these vectors are of size numberofelements)
vx = mean(md.results.StressbalanceSolution.Vx(md.mesh.elements),2)/md.constants.yts;
vy = mean(md.results.StressbalanceSolution.Vy(md.mesh.elements),2)/md.constants.yts;
vb = sqrt(vx.^2+vy.^2);
tau_b = md.friction.C.^2 .* vb.^(1./m);
3. Rescale basal friction for m=100 (or any other values of m)
%Change fricion to m=8 or 100, or anything else
m = 100;
md.friction.m(:) = m;
md.friction.C = sqrt(tau_b./((vb+eps).^(1/m)));
4. Solve a stress balance with m=100, it should converged in 1 single iteration!
%Solve with m=8/100/etc should only take 1 iteration since we start from same stress and converged initial guess
md.initialization.vx = md.results.StressbalanceSolution.Vx;
md.initialization.vy = md.results.StressbalanceSolution.Vy;
md = solve(md,'sb');
I hope this helps!