Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.m =================================================================== --- /issm/trunk-jpl/src/m/mech/analyticaldamage.m (revision 16064) +++ /issm/trunk-jpl/src/m/mech/analyticaldamage.m (revision 16065) @@ -20,8 +20,4 @@ % '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,8 +28,8 @@ % % 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,7 +110,4 @@ 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,109 +133,27 @@ 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; Index: /issm/trunk-jpl/src/m/mech/analyticaldamage.py =================================================================== --- /issm/trunk-jpl/src/m/mech/analyticaldamage.py (revision 16064) +++ /issm/trunk-jpl/src/m/mech/analyticaldamage.py (revision 16065) @@ -26,8 +26,4 @@ '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,8 +34,8 @@ 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,5 +98,5 @@ 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,7 +106,4 @@ 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,110 +128,27 @@ 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