Changeset 16065

Timestamp:

09/04/13 12:10:42 (12 years ago)

Author:

cborstad

Message:

CHG: cleanup of damage function

Location:

issm/trunk-jpl/src/m/mech

Files:

• analyticaldamage.m (modified) (4 diffs)
• analyticaldamage.py (modified) (5 diffs)

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

@@ -20 +20 @@
 %               'damage' which is truncated in the range [0,1-1e-9]
 %
-%               'mask' is a mask defining areas where damage was calculated as negative prior
-%               to truncation, indicating that the ice is too warm or a back stress has not
-%               been accounted for.
-%
 %          '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 
@@ -32 +28 @@
 %
 %   Usage:
-%      [damage,mask,B,backstress]=analyticaldamage(md,options)
+%      [damage,B,backstress]=analyticaldamage(md,options)
 %
 %   Example:
-%      [damage,mask,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3);
+%      [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3);
 
 % check inputs
@@ -114 +110 @@
 a(pos)=-2+1e-3;
 
-%plotmodel(md,'data',a,'caxis',[-2,1])
-%plotmodel(md,'data',ex,'caxis',[-0.5e-9,0.5e-9])
-
 % spreading stress
 rhoi=md.materials.rho_ice;
@@ -140 +133 @@
 end
 
-%D2=1-((theta./ex).^(1./n)).*(T-sigmab)./B;
-%length(find(D2<0))
-%length(find(ex<0))
-%length(find(D2>1))
-%length(find(theta<0))
-%pos=find(ex<0);
-%mask2=ismember(1:md.mesh.numberofvertices,pos);
 D=1-(1+a+a.^2+b.^2).^((n-1)./(2*n))./abs(ex).^(1./n).*(T-sigmab)./B./(2+a)./sign(ex);
-
-%pos=find(D2>1);
-%mask1=ismember(1:md.mesh.numberofvertices,pos);
-%pos=find(a<-2);
-%mask2=ismember(1:md.mesh.numberofvertices,pos);
-%pos=find(ex<-0);
-%mask3=ismember(1:md.mesh.numberofvertices,pos);
-%pos=find(((2+a).*sign(ex))<0);
-%mask4=ismember(1:md.mesh.numberofvertices,pos);
-%plotmodel(md,'nlines',2,'ncols',2,'data',D2,'data',a,'data',ex,'data',(2+a).*sign(ex),...
-%       'mask#1',mask1,'mask#2',mask2,'mask#3',mask3,'mask#4',mask4,...
-%       'expdisp#all','./Exp/LarsenC_DomainOutline_Bedmap2extent.exp',...
-%       'caxis#1',[0 1],'caxis#2',[-2,1],'caxis#3',[0 4e-9],'caxis#4',[0 1]);
 
 backstress=zeros(md.mesh.numberofvertices,1);
 
-% new corrections
 % D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
 pos=find(D>1);
 D(pos)=0;
+
 % backstress to bring damage to zero
 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));
-%backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos));
 
 pos=find(D<0);
 mask=ismember(1:md.mesh.numberofvertices,pos);
-%plotmodel(md,'data',D,'caxis',[0,1],'mask',mask)
 D(pos)=0;
+
 % backstress to bring negative damage to zero
 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));
 
 pos=find(backstress<0);
-%length(pos)
 backstress(pos)=0;
-
-%plotmodel(md,'data',backstress,'caxis',[0 2e5],'edgecolor','k')
 
 % 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));
-%B(pos)=T(pos)./(1-D(pos)).*(theta(pos)./ex(pos)).^(1./n(pos)); % same calculated both ways
 
-% enforce maximum B corresponding to -50 deg C ice
-%Bmax=paterson(273.15-50);
-%posB=find(B>Bmax);
-%B(posB)=Bmax;
-
-%plotmodel(md,'nlines',2,'ncols',1,'data',D,'data',D2,'caxis#all',[0,1])
-%plotmodel(md,'data',paterson(md.initialization.temperature)-B,'caxis',[-1e8 1e8])
-%plotmodel(md,'data',B,'caxis',[0.8e8 5e8])
-%plotmodel(md,'data',D2,'caxis',[0,1])
-%plotmodel(md,'data',(1-D2).*B,'caxis',[0.1e8 2e8])
-
-%plotmodel(md,'data',backstress2./T,'caxis',[1 2],'data',D2)
-%max(D2)
-%min(D2)
-
-%pos=find(D>1);
-%D(pos)=1-1e-9;
-%
-%pos=find(isnan(D));
-%D(pos)=0;
-%
-%% additional calculations here where negative damage is inferred
-%pos=find(D<0);
-%D(pos)=0;
-%
-%% magnitude of back stress to bring negative damage to zero
-%backstress=zeros(md.mesh.numberofvertices,1);
-%backstress(pos)=T(pos)-(1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos));
-%
-%mask=ismember(1:md.mesh.numberofvertices,pos);
-%
-%%% enforce max back stress equal to spreading stress
-%%posS=find(backstress>T); 
-%%backstress(posS)=T(posS);
-%
-%% decreased thickness in areas of negative damage
-%dH=zeros(md.mesh.numberofvertices,1);
-%dH(pos)=1/C*((1-D(pos)).*B(pos).*(ex(pos)./theta(pos)).^(1./n(pos))-T(pos));
-%
-%% higher rigidity to eliminate areas of negative damage
-%B(pos)=T(pos)./(1-D(pos)).*(theta(pos)./ex(pos)).^(1./n(pos));
-%
-%
-%temp=[0:-0.1:-50]+273.15;
-%rig=paterson(temp);
-%newtemp=md.initialization.temperature;
-%
-%for i=1:length(pos)
-%       [mindiff,index]=min(abs(rig-B(pos(i))));
-%       newtemp(pos(i))=temp(index);
-%end
-%
-%dT=newtemp-md.initialization.temperature;
-
-%any(D<0)
-%any(D>1)
-%any(backstress<0)
 damage=D;

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

@@ -26 +26 @@
                         'damage' which is truncated in the range [0,1-1e-9]
         
-                        'mask' is a mask defining areas where damage was calculated as negative prior
-                        to truncation, indicating that the ice is too warm or a back stress has not
-                        been accounted for.
-        
                    '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 
@@ -38 +34 @@
         
            Usage:
-              [damage,mask,B,backstress]=analyticaldamage(md,options)
+              [damage,B,backstress]=analyticaldamage(md,options)
         
            Example:
-              [damage,mask,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3)
+              [damage,B,backstress]=analyticaldamage(md,'eq','Weertman2D','smoothing',2,'backstress',10e3)
         '''
 
@@ -102 +98 @@
         a=ey/ex
         b=exy/ex
-        pos=npy.nonzero(npy.logical_or(ex<0,ey<0))
+        pos=npy.nonzero(npy.logical_and(ex<0,ey<0))
         #length(pos)
         a[pos]=-a[pos]
@@ -110 +106 @@
         pos=npy.nonzero(abs((abs(a)-2))<1e-3)
         a[pos]=-2+1e-3
-        
-        #plotmodel(md,'data',a,'caxis',[-2,1])
-        #plotmodel(md,'data',ex,'caxis',[-0.5e-9,0.5e-9])
         
         # spreading stress
@@ -135 +128 @@
                 raise StandardError('argument passed to "eq" not valid.  Type "help analyticaldamage" for usage')
         
-        #D2=1-((theta./ex).^(1./n)).*(T-sigmab)./B;
-        #length(find(D2<0))
-        #length(find(ex<0))
-        #length(find(D2>1))
-        #length(find(theta<0))
-        #pos=find(ex<0);
-        #mask2=ismember(1:md.mesh.numberofvertices,pos);
         D=1-(1+a+a**2+b**2)**((n-1)/(2*n))/abs(ex)**(1./n)*(T-sigmab)/B/(2+a)/npy.sign(ex)
-        
-        #pos=find(D2>1);
-        #mask1=ismember(1:md.mesh.numberofvertices,pos);
-        #pos=find(a<-2);
-        #mask2=ismember(1:md.mesh.numberofvertices,pos);
-        #pos=find(ex<-0);
-        #mask3=ismember(1:md.mesh.numberofvertices,pos);
-        #pos=find(((2+a).*sign(ex))<0);
-        #mask4=ismember(1:md.mesh.numberofvertices,pos);
-        #plotmodel(md,'nlines',2,'ncols',2,'data',D2,'data',a,'data',ex,'data',(2+a).*sign(ex),...
-        #       'mask#1',mask1,'mask#2',mask2,'mask#3',mask3,'mask#4',mask4,...
-        #       'expdisp#all','./Exp/LarsenC_DomainOutline_Bedmap2extent.exp',...
-        #       'caxis#1',[0 1],'caxis#2',[-2,1],'caxis#3',[0 4e-9],'caxis#4',[0 1]);
         
         backstress=npy.zeros(md.mesh.numberofvertices,)
         
-        # new corrections
         # D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
         pos=npy.nonzero(D>1)
         D[pos]=0
+        
         # backstress to bring damage to zero
         backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
-        #backstress[pos]=T[pos]-(1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos]);
         
         pos=npy.nonzero(D<0)
-        #mask=ismember(1:md.mesh.numberofvertices,pos)
-        #plotmodel(md,'data',D,'caxis',[0,1],'mask',mask)
         D[pos]=0
+
         # backstress to bring negative damage to zero
         backstress[pos]=T[pos]-(1-D[pos])*B[pos]*npy.sign(ex[pos])*(2+a[pos])*abs(ex[pos])**(1./n[pos])/(1+a[pos]+a[pos]**2)**((n[pos]-1)/2./n[pos])
         
         pos=npy.nonzero(backstress<0)
-        #length[pos]
         backstress[pos]=0
-        
-        #plotmodel(md,'data',backstress,'caxis',[0 2e5],'edgecolor','k')
         
         # increased rigidity to bring negative damage to zero
         B[pos]=npy.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]);
-        #B[pos]=T[pos]./(1-D[pos]).*(theta[pos]./ex[pos]).^(1./n[pos]); # same calculated both ways
         
-        # enforce maximum B corresponding to -50 deg C ice
-        #Bmax=paterson(273.15-50);
-        #posB=find(B>Bmax);
-        #B(posB)=Bmax;
-        
-        #plotmodel(md,'nlines',2,'ncols',1,'data',D,'data',D2,'caxis#all',[0,1])
-        #plotmodel(md,'data',paterson(md.initialization.temperature)-B,'caxis',[-1e8 1e8])
-        #plotmodel(md,'data',B,'caxis',[0.8e8 5e8])
-        #plotmodel(md,'data',D2,'caxis',[0,1])
-        #plotmodel(md,'data',(1-D2).*B,'caxis',[0.1e8 2e8])
-        
-        #plotmodel(md,'data',backstress2./T,'caxis',[1 2],'data',D2)
-        #max(D2)
-        #min(D2)
-        
-        #pos=find(D>1);
-        #D[pos]=1-1e-9;
-        #
-        #pos=find(isnan(D));
-        #D[pos]=0;
-        #
-        ## additional calculations here where negative damage is inferred
-        #pos=find(D<0);
-        #D[pos]=0;
-        #
-        ## magnitude of back stress to bring negative damage to zero
-        #backstress=zeros(md.mesh.numberofvertices,1);
-        #backstress[pos]=T[pos]-(1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos]);
-        #
-        #mask=ismember(1:md.mesh.numberofvertices,pos);
-        #
-        ## enforce max back stress equal to spreading stress
-        ##posS=find(backstress>T); 
-        ##backstress(posS)=T(posS);
-        #
-        ## decreased thickness in areas of negative damage
-        #dH=zeros(md.mesh.numberofvertices,1);
-        #dH[pos]=1/C*((1-D[pos]).*B[pos].*(ex[pos]./theta[pos]).^(1./n[pos])-T[pos]);
-        #
-        ## higher rigidity to eliminate areas of negative damage
-        #B[pos]=T[pos]./(1-D[pos]).*(theta[pos]./ex[pos]).^(1./n[pos]);
-        #
-        #
-        #temp=[0:-0.1:-50]+273.15;
-        #rig=paterson(temp);
-        #newtemp=md.initialization.temperature;
-        #
-        #for i=1:length(pos)
-        #       [mindiff,index]=min(abs(rig-B(pos(i))));
-        #       newtemp(pos(i))=temp(index);
-        #end
-        #
-        #dT=newtemp-md.initialization.temperature;
-        
-        #any(D<0)
-        #any(D>1)
-        #any(backstress<0)
         damage=D