Changeset 26443

Timestamp:

09/20/21 13:21:48 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: rewriting/reorganizing crevasse depth law

File:

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

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

@@ -405 +405 @@
 
         IssmDouble  calvingrate[NUMVERTICES];
-        IssmDouble  vx,vy,vel;
+        IssmDouble  vx,vy;
         IssmDouble  water_height, bed,Ho,thickness,surface;
         IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
-        IssmDouble  B, strainparallel, straineffective,n;
+        IssmDouble  strainparallel, straineffective,B,n;
         IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
-        int crevasse_opening_stress;
+        int         crevasse_opening_stress;
 
         /*retrieve the type of crevasse_opening_stress*/
@@ -442 +442 @@
                 bed_input->GetInputValue(&bed,&gauss);
                 surface_input->GetInputValue(&surface,&gauss);
-                strainrateparallel_input->GetInputValue(&strainparallel,&gauss);
-                strainrateeffective_input->GetInputValue(&straineffective,&gauss);
+
                 vx_input->GetInputValue(&vx,&gauss);
                 vy_input->GetInputValue(&vy,&gauss);
@@ -450 +449 @@
                 s_xy_input->GetInputValue(&s_xy,&gauss);
                 s_yy_input->GetInputValue(&s_yy,&gauss);
-                B_input->GetInputValue(&B,&gauss);
-                n_input->GetInputValue(&n,&gauss);
-
-                vel=sqrt(vx*vx+vy*vy)+1.e-14;
-
-                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.;
-
-                if(crevasse_opening_stress==0){         /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
-                        surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / n)-1)) / (rho_ice * constant_g);
-                        basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/n)-1)) / (rho_ice*constant_g) - Ho);
-                }
-                else if(crevasse_opening_stress==1){     /*Benn2017,Todd2018: maximum principal stress */
-                        surface_crevasse[iv] = s1 / (rho_ice*constant_g);
-                        basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
-                }
-
-                /* some constraints */
-                if (surface_crevasse[iv]<0.) {
+
+      /*Get longitudinal or maximum Eigen stress*/
+                if(crevasse_opening_stress==0){
+         /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
+                        strainrateparallel_input->GetInputValue(&strainparallel,&gauss);
+                        strainrateeffective_input->GetInputValue(&straineffective,&gauss);
+                        B_input->GetInputValue(&B,&gauss);
+                        n_input->GetInputValue(&n,&gauss);
+         s1 =  B * strainparallel * pow(straineffective, (1./n)-1);
+                }
+                else if(crevasse_opening_stress==1){
+         /*Benn2017,Todd2018: maximum principal stress */
+                        Matrix2x2Eigen(&s1,&s2,NULL,NULL,s_xx,s_xy,s_yy);
+                        if(fabs(s2)>fabs(s1)){
+            stmp=s2; s2=s1; s1=stmp;
+         }
+                }
+      else{
+         _error_("not supported");
+      }
+
+      /*Surface crevasse: sigma'_xx - rho_i g d + rho_fw g d_w = 0*/
+      surface_crevasse[iv] = s1 / (rho_ice*constant_g) + (rho_freshwater/rho_ice)*water_height;
+                if(surface_crevasse[iv]<0.){
                         surface_crevasse[iv]=0.;
                         water_height = 0.;
                 }
-                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
-                if (bed>0.) basal_crevasse[iv] = 0.;
-
-                //if (surface_crevasse[iv]<water_height){
-                //      water_height = surface_crevasse[iv];
-                //}
-
-                /* add water in surface crevasse */
-                surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */
-                crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */
-
+
+      /*Basal crevasse: sigma'_xx - rho_i g (H-d) - rho_w g (b+d) = 0*/
+      if(bed>0.){
+         basal_crevasse[iv] = 0.;
+      }
+      else{
+         Ho = thickness - (rho_seawater/rho_ice) * (-bed);
+         if(Ho<0.)  Ho=0.;
+         basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
+         if(basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
+      }
+
+      /*Total crevasse depth (surface + basal)*/
+                crevasse_depth[iv]   = surface_crevasse[iv] + basal_crevasse[iv];
         }