Changeset 17367

Timestamp:

02/28/14 12:03:48 (11 years ago)

Author:

cborstad

Message:

CHG: added a model field to choose the type of scalar equivalent stress for a damage evolution solution

Location:

issm/trunk-jpl/src

Files:

• c/analyses/DamageEvolutionAnalysis.cpp (modified) (4 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• m/classes/damage.m (modified) (5 diffs)
• m/classes/damage.py (modified) (7 diffs)
• m/classes/model.m (modified) (1 diff)
• m/enum/DamageEquivStressEnum.m (added)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/mech/mechanicalproperties.m (modified) (2 diffs)
• m/mech/steadystateiceshelftemp.m (modified) (1 diff)
• m/mech/thomasparams.m (modified) (6 diffs)
• m/plot/colormaps/getcolormap.m (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

@@ -31 +31 @@
                 parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
                 parameters->AddObject(iomodel->CopyConstantObject(DamageHealingEnum));
+                parameters->AddObject(iomodel->CopyConstantObject(DamageEquivStressEnum));
         }
         xDelete<char>(law);
@@ -395 +396 @@
         IssmDouble c1,c2,c3,healing,stress_threshold;
         IssmDouble s_xx,s_xy,s_yy;
-        IssmDouble J2s,Xis,Psi,PosPsi,NegPsi;
+        IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
         IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
+        int equivstress;
 
         /*Fetch number of vertices and allocate output*/
@@ -418 +420 @@
         Input* damage_input    = element->GetInput(DamageDbarEnum); _assert_(damage_input);
 
-        /*Damage evolution z mapping: */
+        /*retrieve the desired type of equivalent stress*/
+        element->FindParam(&equivstress,DamageEquivStressEnum);
+
+        /*Calculate damage evolution source term: */
         Gauss* gauss=element->NewGauss();
         for (int iv=0;iv<numvertices;iv++){
@@ -432 +437 @@
                 s_yy=sigma_yy/(1.-damage);
 
-                J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
-                Xis=sqrt(3.0)*J2s;
-                Psi=Xis-stress_threshold;
+                if(equivstress==1){ /* von Mises */
+                        J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
+                        Chi=sqrt(3.0)*J2s;
+                }
+                Psi=Chi-stress_threshold;
                 PosPsi=max(Psi,0.);
                 NegPsi=max(-Psi,0.);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -179 +179 @@
         DamageSpcdamageEnum,
         DamageMaxDamageEnum,
+        DamageEquivStressEnum,
         MaterialsRhoIceEnum,
         MaterialsRhoWaterEnum,
@@ -671 +672 @@
         OptionStructEnum,
         /*}}}*/
-        /*Rheology law (move too Material) {{{*/
+        /*Rheology law (move to Material) {{{*/
         PatersonEnum,
         ArrheniusEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -187 +187 @@
                 case DamageSpcdamageEnum : return "DamageSpcdamage";
                 case DamageMaxDamageEnum : return "DamageMaxDamage";
+                case DamageEquivStressEnum : return "DamageEquivStress";
                 case MaterialsRhoIceEnum : return "MaterialsRhoIce";
                 case MaterialsRhoWaterEnum : return "MaterialsRhoWater";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -190 +190 @@
               else if (strcmp(name,"DamageSpcdamage")==0) return DamageSpcdamageEnum;
               else if (strcmp(name,"DamageMaxDamage")==0) return DamageMaxDamageEnum;
+              else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
               else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
@@ -259 +260 @@
               else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
               else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
-              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
+              if (strcmp(name,"Surface")==0) return SurfaceEnum;
+              else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
               else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
               else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
@@ -382 +383 @@
               else if (strcmp(name,"Vertices")==0) return VerticesEnum;
               else if (strcmp(name,"Results")==0) return ResultsEnum;
-              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum;
+              if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
+              else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum;
               else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
               else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
@@ -505 +506 @@
               else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
               else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
-              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
+              if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
+              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
               else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
               else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
@@ -628 +629 @@
               else if (strcmp(name,"Separate")==0) return SeparateEnum;
               else if (strcmp(name,"Sset")==0) return SsetEnum;
-              else if (strcmp(name,"Verbose")==0) return VerboseEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
+              if (strcmp(name,"Verbose")==0) return VerboseEnum;
+              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
               else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
               else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;

issm/trunk-jpl/src/m/classes/damage.m

@@ -26 +26 @@
                 c4                  = NaN;
                 healing             = NaN;
+                equiv_stress              = NaN;
         end
         methods
@@ -76 +77 @@
                         obj.c3=0;
                         obj.c4=0;
+                        obj.equiv_stress=0;
 
                 end % }}}
@@ -98 +100 @@
                                 md = checkfield(md,'fieldname','damage.c4','>=',0);
                                 md = checkfield(md,'fieldname','damage.stress_threshold','>=',0);
+                                md = checkfield(md,'fieldname','damage.equiv_stress','numel',[1],'values',[0]);
                         elseif strcmpi(obj.law,'undamaged'),
                                 if (solution==DamageEvolutionSolutionEnum),
@@ -128 +131 @@
                                 fielddisplay(obj,'healing','damage healing parameter 1');
                                 fielddisplay(obj,'stress_threshold','damage stress threshold [Pa]');
+                                fielddisplay(obj,'equiv_stress','0: von Mises');
                         end
 
@@ -151 +155 @@
                                 WriteData(fid,'object',obj,'fieldname','stress_threshold','format','Double');
                                 WriteData(fid,'object',obj,'fieldname','healing','format','Double');
+                                WriteData(fid,'object',obj,'fieldname','equiv_stress','format','Integer');
                         end
 

issm/trunk-jpl/src/m/classes/damage.py

@@ -35 +35 @@
                 self.c4                 = float('NaN')
                 self.healing                            = float('NaN')
+                self.equiv_stress       = float('NaN')
 
                 if not len(args):
@@ -62 +63 @@
                         s+="%s\n" % fielddisplay(self,"c4","damage parameter 4 ")
                         s+="%s\n" % fielddisplay(self,"stress_threshold","damage stress threshold [Pa]")
+                        s+="%s\n" % fielddisplay(self,"equiv_stresss","0: von Mises")
 
                 return s
@@ -95 +97 @@
                 self.c4=0
                 self.healing=0
+                self.equiv_stress=0
 
         # }}}
@@ -104 +107 @@
                 md = checkfield(md,'fieldname','damage.spcdamage','forcing',1)
                         
-                md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2]);
+                md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2])
                 md = checkfield(md,'fieldname','damage.maxiter','>=0',0);
                 md = checkfield(md,'fieldname','damage.penalty_factor','>=0',0);
@@ -111 +114 @@
 
                 if self.law == 'pralong':
-                        md = checkfield(md,'fieldname','damage.healing','>=',0);
+                        md = checkfield(md,'fieldname','damage.healing','>=',0)
                         md = checkfield(md,'fieldname','damage.c1','>=',0)
                         md = checkfield(md,'fieldname','damage.c2','>=',0)
@@ -117 +120 @@
                         md = checkfield(md,'fieldname','damage.c4','>=',0)
                         md = checkfield(md,'fieldname','damage.stress_threshold','>=',0)
+                        md = checkfield(md,'fieldname','damage.equiv_stress','numel',[1],'values',[0])
                 elif strcmpi(self.law,'undamaged'):
                         if (solution==DamageEvolutionSolutionEnum):
@@ -142 +146 @@
                         WriteData(fid,'object',self,'fieldname','c4','format','Double')
                         WriteData(fid,'object',self,'fieldname','stress_threshold','format','Double')
+                        WriteData(fid,'object',self,'fieldname','equiv_stress','format','Integer')
         # }}}

issm/trunk-jpl/src/m/classes/model.m

@@ -1107 +1107 @@
                         disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings'));
                         disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(obj.materials) ']'],'material properties'));
-                        disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(obj.damage) ']'],'damage propagation laws'));
+                        disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(obj.damage) ']'],'parameters for damage evolution solution'));
                         disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties'));
                         disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(obj.flowequation) ']'],'flow equations'));

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -179 +179 @@
 def DamageSpcdamageEnum(): return StringToEnum("DamageSpcdamage")[0]
 def DamageMaxDamageEnum(): return StringToEnum("DamageMaxDamage")[0]
+def DamageEquivStressEnum(): return StringToEnum("DamageEquivStress")[0]
 def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
 def MaterialsRhoWaterEnum(): return StringToEnum("MaterialsRhoWater")[0]

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

@@ -15 +15 @@
 %some checks
 if length(vx)~=md.mesh.numberofvertices | length(vy)~=md.mesh.numberofvertices,
-        error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!'])
+        %error(['the input velocity should be of size ' num2str(md.mesh.numberofvertices) '!'])
 end
 if ~(md.mesh.dimension==2)
@@ -139 +139 @@
 deviatoricstress.effectivevalue=1/sqrt(2)*sqrt(stress.xx.^2+stress.yy.^2+2*stress.xy.^2);
 md.results.deviatoricstress=deviatoricstress;
+
+viscosity=struct('nu',[]);
+viscosity.nu=nu;
+md.results.viscosity=viscosity;

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

@@ -46 +46 @@
 
 %calculate depth-averaged temperature (in Celsius)
-%temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1));
-temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki)));
+temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1));
+%temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki)));
 
 %temperature should not be less than surface temp

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

@@ -1 @@
-function [alpha,beta,theta,ex]=thomasparams(md,varargin)
+function [alpha,beta,theta,ex,sigxx]=thomasparams(md,varargin)
 %THOMASPARAMS - compute Thomas' geometric parameters for an ice shelf 
 %
@@ -35 +35 @@
 %               the equivalent stress
 %
-%               'exx' is the longitudinal strain rate along a coordinate system defined
-%               by 'coordsys' 
+%               'exx' is the strain rate along a coordinate system defined by 'coordsys' 
 %
-%   Usage: [alpha,beta,theta,exx]=ThomasParams(md,options)
+%               'sigxx' is the deviatoric stress along a coordinate system defined by 'coordsys' 
 %
-%   Example: [alpha,beta,theta,exx]=ThomasParams(md,'eq','Thomas','smoothing',2,'coordsys','longitudinal')
+%   Usage: [alpha,beta,theta,exx,sigxx]=ThomasParams(md,options)
+%
+%   Example: [alpha,beta,theta,exx,sigxx]=ThomasParams(md,'eq','Thomas','smoothing',2,'coordsys','longitudinal')
 
 %some checks
@@ -82 +83 @@
 end
 
+% rheology
+n=averaging(md,md.materials.rheology_n,0);
+B=md.materials.rheology_B;
+
 switch coordsys
         case 'principal'
@@ -98 +103 @@
                 id=find(e1<0 & e2<0);
                 a(id)=-a(id); % where both strain rates are compressive, enforce negative alpha
+                sigxx=(abs(ex)./((1+a+a.^2).^((n-1)/2))).^(1./n).*B;
                 %mask=ismember(1:md.mesh.numberofvertices,id);
                 %plotmodel(md,'data',ex,'mask',mask)
@@ -114 +120 @@
                 a=ey./ex;
                 b=exy./ex;
-                pos=find(ex<0 & ey<0);
-                %length(pos)
-                a(pos)=-a(pos);
+                %pos=find(ex<0 & ey<0);
+                %a(pos)=-a(pos);
+                %sigxx=(abs(ex)./((1+a+a.^2+b.^2).^((n-1)/2))).^(1./n).*B;
+                sigxx=abs(ex).^(1./n-1).*ex./((1+a+a.^2+b.^2).^((n-1)./(2*n))).*B;
         otherwise
                 error('argument passed to "coordsys" not valid');
@@ -128 +135 @@
 end
 a(pos)=-2+1e-3;
-
-% rheology
-n=averaging(md,md.materials.rheology_n,0);
 
 switch eq

issm/trunk-jpl/src/m/plot/colormaps/getcolormap.m

@@ -41 +41 @@
 elseif strcmpi(map,'Rignot'),
         alpha=getfieldvalue(options,'alpha',1);
-        map = hsv;
+        map = hsv(128);
         map = rgb2hsv(map);
         map(:,2) = max(min( (0.1+map(:,1)).^(1/alpha) ,1),0);