Changeset 22488 for issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

Timestamp:

03/01/18 14:48:25 (7 years ago)

Author:

youngmc3

Message:

CHG: updated crevasse-depth calving law

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -302 +302 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  vx,vy,vel;
-        IssmDouble  rheology_B,critical_fraction,water_height;
-        IssmDouble  bathymetry,Ho,thickness,float_depth,groundedice;
+        IssmDouble  critical_fraction,water_height;
+        IssmDouble  bed,Ho,thickness,float_depth;
         IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
         IssmDouble  strainparallel, straineffective;
-        IssmDouble  yts;
+        IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
 
         /* Get node coordinates and dof list: */
@@ -313 +313 @@
         /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
         this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
-        this->parameters->FindParam(&yts,ConstantsYtsEnum);
                 
         IssmDouble rho_ice        = this->GetMaterialParameter(MaterialsRhoIceEnum);
@@ -321 +320 @@
         IssmDouble rheology_n     = this->GetMaterialParameter(MaterialsRheologyNEnum);
 
-        Input*   B_input       = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
-        Input*   H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
-        Input*   bed           = inputs->GetInput(BedEnum); _assert_(bed);
-        Input*   surface       = inputs->GetInput(SurfaceEnum); _assert_(surface);
-        Input*  strainrateparallel  = inputs->GetInput(StrainRateparallelEnum);  _assert_(strainrateparallel);
-        Input*  strainrateeffective = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective);
-        Input*  vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
-        Input*  vy_input = inputs->GetInput(VxEnum); _assert_(vy_input);
-        Input*  gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
-        Input*   waterheight_input = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input);
-
+        Input*   H_input                 = inputs->GetInput(ThicknessEnum); _assert_(H_input);
+        Input*   bed_input               = inputs->GetInput(BedEnum); _assert_(bed_input);
+        Input*   surface_input           = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+        Input*  strainrateparallel_input  = inputs->GetInput(StrainRateparallelEnum);  _assert_(strainrateparallel_input);
+        Input*  strainrateeffective_input = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective_input);
+        Input*  vx_input                  = inputs->GetInput(VxEnum); _assert_(vx_input);
+        Input*  vy_input                  = inputs->GetInput(VxEnum); _assert_(vy_input);
+        Input*   waterheight_input       = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input);
+        Input*   s_xx_input              = inputs->GetInput(DeviatoricStressxxEnum);     _assert_(s_xx_input);
+        Input*   s_xy_input              = inputs->GetInput(DeviatoricStressxyEnum);     _assert_(s_xy_input);
+        Input*   s_yy_input              = inputs->GetInput(DeviatoricStressyyEnum);     _assert_(s_yy_input);
+        
         /*Loop over all elements of this partition*/
         GaussTria* gauss=new GaussTria();
@@ -337 +337 @@
                 gauss->GaussVertex(iv);
         
-                B_input->GetInputValue(&rheology_B,gauss);
                 H_input->GetInputValue(&thickness,gauss);
-                bed->GetInputValue(&bathymetry,gauss);
-                surface->GetInputValue(&float_depth,gauss);
-                strainrateparallel->GetInputValue(&strainparallel,gauss);
-                strainrateeffective->GetInputValue(&straineffective,gauss);
+                bed_input->GetInputValue(&bed,gauss);
+                surface_input->GetInputValue(&float_depth,gauss);
+                strainrateparallel_input->GetInputValue(&strainparallel,gauss);
+                strainrateeffective_input->GetInputValue(&straineffective,gauss);
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
-                gr_input->GetInputValue(&groundedice,gauss);
                 waterheight_input->GetInputValue(&water_height,gauss);
+                s_xx_input->GetInputValue(&s_xx,gauss);
+                s_xy_input->GetInputValue(&s_xy,gauss);
+                s_yy_input->GetInputValue(&s_yy,gauss);
+                
                 vel=sqrt(vx*vx+vy*vy)+1.e-14;
 
-                Ho = thickness - (rho_seawater/rho_ice) * (-bathymetry);
+                s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
+                s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
+                if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
+                
+                Ho = thickness - (rho_seawater/rho_ice) * (-bed);
                 if(Ho<0.)  Ho=0.;
 
                 /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
                 /*surface crevasse*/
-                surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
-                //surface_crevasse[iv] = 2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice * constant_g);
+                //surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
+                surface_crevasse[iv] = s1 / (rho_ice*constant_g);
                 if (surface_crevasse[iv]<0.) {
                         surface_crevasse[iv]=0.;
                         water_height = 0.;
                 }
-                if (surface_crevasse[iv]<water_height){
-                        water_height = surface_crevasse[iv];
-                }
+                //if (surface_crevasse[iv]<water_height){
+                //      water_height = surface_crevasse[iv];
+                //}
                 
                 /*basal crevasse*/
-                basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
-                //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice*constant_g) - Ho);
+                //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
+                basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
                 if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
-                if (bathymetry>0.) basal_crevasse[iv] = 0.; 
+                if (bed>0.) basal_crevasse[iv] = 0.; 
         
                 H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
@@ -374 +380 @@
                 
                 crevasse_depth[iv] = max(H_surf,H_surfbasal);
+        }
         
-                /*Assign values */
-                if(crevasse_depth[iv]>=0. && bathymetry<=0.){
-                 calvingrate[iv] = vel+3000./yts;
-                }       
-                else
-                 calvingrate[iv]=0.;
-        }
-
         this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
         this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum));
         this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum));
-        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
 
         delete gauss;