Changeset 16416

Timestamp:

10/15/13 15:34:56 (11 years ago)

Author:

cborstad

Message:

NEW: functions for calculating ice shelf damage and backstress from inversion (B) results.

Location:

issm/trunk-jpl/src/m/mech

Files:

• analyticaldamage.m (modified) (6 diffs)
• backstressfrominversion.m (added)
• damagefrominversion.m (added)
• thomasparams.m (added)

issm/trunk-jpl/src/m/mech/analyticaldamage.m

@@ -13 +13 @@
 %                                                               'Thomas' for a 2D ice shelf, taking into account full strain rate tensor (default)
 %               - 'smoothing'   : the amount of smoothing to be applied to the strain rate data.
-%                                                               Type 'help averaging' for more information on its usage.
-%               - 'sigmab'              : a compressive backstress term to be subtracted from the driving stress 
-%                                                               in the damage calculation
+%                                                               Type 'help averaging' for more information on its
+%                                                               usage. Defaults to 0.
+%               - 'sigmab'              : a-priori backstress in opposition to the driving stress.
+%                                                               Defaults to 0 everywhere. 
+
+%               - 'coordsys'    : coordinate system for calculating the strain rate
+%                                                       components. Must be one of: 
+%                               'longitudinal': x axis aligned along a flowline at every point (default)
+%                               'principal': x axis aligned along maximum principal strain rate
+%                                       at every point
+%                               'xy': x and y axes same as in polar stereographic projection 
 %
 %   Return values:
 %               'damage' which is truncated in the range [0,1-1e-9]
 %
-%          'B' is the rigidity, which is equal to md.materials.rheology_B in areas outside
-%               those defined by 'mask.'  Within areas defined by 'mask,' where negative damage 
-%               is inferred, 'B' is updated to make damage equal to zero.  
+%          'B' is the ice rigidity calculated assuming D=0 everywhere. 
 %
-%               'backstress' is the inferred backstress necessary to balance the analytical solution
-%               (keeping damage within its appropriate limits, e.g. D in [0,1]).
+%               'backstress' is the inferred backstress necessary to balance the
+%               analytical solution (keeping damage within its appropriate limits, e.g. D
+%               in [0,1]).
 %
 %   Usage:
@@ -31 +38 @@
 %
 %   Example:
-%      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3);
+%      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'sigmab',10e3,'coordsys','longitudinal');
 
 % check inputs
@@ -53 +60 @@
 smoothing = getfieldvalue(options,'smoothing',0);
 sigmab = getfieldvalue(options,'sigmab',0);
+coordsys = getfieldvalue(options,'coordsys','longitudinal');
 if length(sigmab)==1,
         sigmab=sigmab*ones(md.mesh.numberofelements,1);
 end
 
-% average element strain rates onto vertices
-e1=averaging(md,md.results.strainrate.principalvalue1,smoothing)/md.constants.yts; % convert to s^-1
-e2=averaging(md,md.results.strainrate.principalvalue2,smoothing)/md.constants.yts;
-exx=averaging(md,md.results.strainrate.xx,smoothing)/md.constants.yts;
-eyy=averaging(md,md.results.strainrate.yy,smoothing)/md.constants.yts;
-exy=averaging(md,md.results.strainrate.xy,smoothing)/md.constants.yts;
-
-% checks: any of e1 or e2 equal to zero?
-pos=find(e1==0);
-if any(pos==1)
-        disp('WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1');
-        e1(pos)=1e-13;
-end
-pos=find(e2==0);
-if any(pos==1)
-        disp('WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1');
-        e2(pos)=1e-13;
-end
-
-%% old method using principal strain rates {{{
-%% ex=maximum principal tensile strain rate
-%ex=e1;
-%a=e2./e1;
-%pos=find(e1<0 & e2>0); % longitudinal compression and lateral tension
-%a(pos)=e1(pos)./e2(pos);
-%ex(pos)=e2(pos);
-%pos2=find(e1<0 & e2<0 & abs(e1)<abs(e2)); % lateral and longitudinal compression
-%a(pos2)=e1(pos2)./e2(pos2);
-%ex(pos2)=e2(pos2);
-%pos3=find(e1>0 & e2>0 & abs(e1)<abs(e2)); % lateral and longitudinal tension 
-%a(pos3)=e1(pos3)./e2(pos3);
-%ex(pos3)=e2(pos3);
-%id=find(e1<0 & e2<0);
-%a(id)=-a(id); % where both strain rates are compressive, enforce negative alpha
-%
-%% }}}
-
-% new method using longitudinal strain rates defined by observed velocity vector
-velangle=atan(md.initialization.vy./md.initialization.vx);
-ex=0.5*(exx+eyy)+0.5*(exx-eyy).*cos(2*velangle)+exy.*sin(2*velangle);
-ey=exx+eyy-ex; % trace of strain rate tensor is invariant
-exy=-0.5*(exx-eyy).*sin(2*velangle)+exy.*cos(2*velangle);
-a=ey./ex;
-b=exy./ex;
-pos=find(ex<0 & ey<0);
-%length(pos)
-a(pos)=-a(pos);
-
-% a < -1 in areas of strong lateral compression or longitudinal compression
-% and theta is undefined at a = -2
-pos=find(abs((abs(a)-2))<1e-3);
-a(pos)=-2+1e-3;
+[a,b,theta,ex]=thomasparams(md,'eq',eq,'smoothing',smoothing,'coordsys',coordsys);
 
 % spreading stress
@@ -120 +77 @@
 n=averaging(md,md.materials.rheology_n,0);
 
-switch eq
-        case 'Weertman1D'
-                theta=1./8*ones(md.mesh.numberofvertices,1);
-                a=zeros(md.mesh.numberofvertices,1);
-        case 'Weertman2D'
-                theta=1./9*ones(md.mesh.numberofvertices,1);
-                a=ones(md.mesh.numberofvertices,1);
-        case 'Thomas'
-                theta=((1+a+a.^2+b.^2).^((n-1)/2))./(abs(2+a).^n);
-        otherwise
-                error('argument passed to "eq" not valid.  Type "help analyticaldamage" for usage');
-end
-
 D=1-(1+a+a.^2+b.^2).^((n-1)./(2*n))./abs(ex).^(1./n).*(T-sigmab)./B./(2+a)./sign(ex);
-
-backstress=zeros(md.mesh.numberofvertices,1);
 
 % D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
@@ -141 +83 @@
 D(pos)=0;
 
-% backstress to bring damage to zero
+backstress=zeros(md.mesh.numberofvertices,1);
+
+% backstress to bring D down to one 
 backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*sign(ex(pos)).*(2+a(pos)).*abs(ex(pos)).^(1./n(pos))./(1+a(pos)+a(pos).^2).^((n(pos)-1)/2./n(pos));
 
@@ -154 +98 @@
 backstress(pos)=0;
 
-% increased rigidity to bring negative damage to zero
-B(pos)=sign(ex(pos))./(2+a(pos)).*(1+a(pos)+a(pos).^2).^((n(pos)-1)/2./n(pos)).*T(pos)./abs(ex(pos)).^(1./n(pos));
+% rigidity from Thomas relation for D=0 and backstress=0
+B=sign(ex)./(2+a).*(1+a+a.^2).^((n-1)/2./n).*T./(abs(ex).^(1./n));
+pos=find(B<0);
+B(pos)=md.materials.rheology_B(pos);
 
 damage=D;