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

Timestamp:

10/20/17 16:28:00 (7 years ago)

Author:

youngmc3

Message:

CHG: adding calving parameterizations (crevasse-depth and height-above-buoyancy laws), minor changes for calvingdev and levermann

File:

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

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

@@ -216 +216 @@
         IssmDouble  calvingrate[NUMVERTICES];
         IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
-        IssmDouble  sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded;
+        IssmDouble  sigma_vm[NUMVERTICES];
+        IssmDouble  sigma_max,sigma_max_floating,sigma_max_grounded;
         IssmDouble  epse_2,groundedice,bed;
 
@@ -260 +261 @@
                 /*Calculate sigma_vm*/
                 epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
-                sigma_vm  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
+                sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
 
                 /*OLD (keep for a little bit)*/
@@ -279 +280 @@
                 }
                 else{
-                        calvingratex[iv]=vx*sigma_vm/sigma_max;
-                        calvingratey[iv]=vy*sigma_vm/sigma_max;
+                        calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
+                        calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
                 }
                 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
@@ -289 +290 @@
         this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
         this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(SigmaVMEnum,&sigma_vm[0],P1Enum));
 
         /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
+void       Tria::CalvingCrevasseDepth(){/*{{{*/
+        
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  vx,vy,vel;
+        IssmDouble  critical_fraction,water_height;
+        IssmDouble  bathymetry,Ho,thickness,float_depth,groundedice;
+        IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
+        IssmDouble  strainparallel, straineffective;
+        IssmDouble  yts;
+
+        /* Get node coordinates and dof list: */
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+                
+        /*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);
+        IssmDouble rho_seawater   = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+        IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
+        IssmDouble constant_g     = this->GetMaterialParameter(ConstantsGEnum);
+        IssmDouble rheology_B     = this->GetMaterialParameter(MaterialsRheologyBbarEnum);
+        IssmDouble rheology_n     = this->GetMaterialParameter(MaterialsRheologyNEnum);
+        
+        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);
+
+        /*Loop over all elements of this partition*/
+        GaussTria* gauss=new GaussTria();
+        for (int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+        
+                H_input->GetInputValue(&thickness,gauss);
+                bed->GetInputValue(&bathymetry,gauss);
+                surface->GetInputValue(&float_depth,gauss);
+                strainrateparallel->GetInputValue(&strainparallel,gauss);
+                strainrateeffective->GetInputValue(&straineffective,gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                gr_input->GetInputValue(&groundedice,gauss);
+                waterheight_input->GetInputValue(&water_height,gauss);
+                vel=sqrt(vx*vx+vy*vy)+1.e-14;
+
+                Ho = thickness - (rho_seawater/rho_ice) * (-bathymetry);
+                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);
+                if (surface_crevasse[iv]<0.) {
+                        surface_crevasse[iv]=0.;
+                        water_height = 0.;
+                }
+                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);
+                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
+                if (bathymetry>0.) basal_crevasse[iv] = 0.; 
+        
+                H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
+                H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness);
+                
+                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;
 }
@@ -331 +426 @@
 
                 /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
-                calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
-                if(calvingrate[iv]<0. || bed>0.){
+                if(strainparallel>0. && strainperpendicular>0. && bed<=0.){
+                        calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
+                }
+                else
                         calvingrate[iv]=0.;
-                }
+                
                 calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
                 calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);