Changeset 27252

Timestamp:

08/31/22 12:50:05 (3 years ago)

Author:

Cheng Gong

Message:

CHG: use max Vel to normalize the parameterization of the frontal ablation rate

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (5 diffs)
• c/cores/movingfront_core.cpp (modified) (2 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• m/classes/calvingparameterization.m (modified) (4 diffs)

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

@@ -226 +226 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.yoffset",CalvingYoffsetEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_lowerbound",CalvingVelLowerboundEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_threshold",CalvingVelThresholdEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_upperbound",CalvingVelUpperboundEnum));
                         break;

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

@@ -845 +845 @@
         int                     use_parameter=-1;
         IssmDouble  gamma, theta, alpha, xoffset, yoffset;
-        IssmDouble  vel_lower, vel_upper, vrate, truncateVrate;
+        IssmDouble  vel_lower, vel_upper, vel_threshold, vrate, truncateVrate;
 
         /* Get node coordinates and dof list: */
@@ -873 +873 @@
         this->FindParam(&vel_lower, CalvingVelLowerboundEnum);
         this->FindParam(&vel_upper, CalvingVelUpperboundEnum);
+        this->FindParam(&vel_threshold, CalvingVelThresholdEnum);
 
         /* Start looping on the number of vertices: */
@@ -891 +892 @@
                 arate_input->GetInputValue(&arate,gauss);
                 vrate = 1.0;
-                if (vel < vel_upper) vrate = vel / vel_upper;
+                if (vel < vel_threshold) vrate = vel / vel_threshold;
 
                 /*Compute strain rate and viscosity: */
@@ -924 +925 @@
                         case 3:
                                 /* 3 tanh(normal vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
-                                _error_("The normalized velocity is not supported yet!");
-                                gamma = yoffset -  0.5*theta*tanh(alpha*(vel+xoffset));
+                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / (vel_upper - vel_lower);
+                                gamma = yoffset -  0.5*theta*tanh(alpha*(truncateVrate+xoffset));
                                 break;
                         case 4:
-                                /* 4 tanh(truncated vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
-                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / vel_upper;
+                                /* 4 tanh(normal vel), cross (0,0): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
+                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / (vel_upper - vel_lower);
                                 gamma = 0.5*theta*(tanh(alpha*xoffset) - tanh(alpha*(truncateVrate+xoffset)));
                                 break;
@@ -946 +947 @@
 
                 /*-------------------------------------------*/
-                if (use_parameter < 3) {
-                        calvingrate[iv] = arate*gamma*vrate;
-                } 
-                else {
-                        calvingrate[iv] = arate*gamma;
-                }
+                calvingrate[iv] = arate*gamma*vrate;
         }
         /*Add input*/

issm/trunk-jpl/src/c/cores/movingfront_core.cpp

@@ -20 +20 @@
         int* extrapol_vars=NULL;
         Analysis  *analysis=NULL;
+        IssmDouble maxVel;
 
         /* recover parameters */
@@ -50 +51 @@
         /* smoothen slope of lsf for computation of normal on ice domain*/
         levelsetfunctionslope_core(femmodel);
+
+        /* compute the maximal velocity over the whole domain */
+        femmodel->MaxVelx(&maxVel);
+        femmodel->parameters->SetParam(maxVel, CalvingVelThresholdEnum);
 
         /* start the work from here */

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

@@ -119 +119 @@
         CalvingYoffsetEnum,
         CalvingVelLowerboundEnum,
+        CalvingVelThresholdEnum,
         CalvingVelUpperboundEnum,
         ConfigurationTypeEnum,

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

@@ -12 +12 @@
                 xoffset = 0;
                 yoffset = 0;
-                vel_upperbound = 6000;
+                vel_upperbound = 8000;
+                vel_threshold = 6000;
                 vel_lowerbound = 0;
         end
@@ -68 +69 @@
                         md = checkfield(md,'fieldname','calving.yoffset','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.vel_lowerbound','NaN',1,'Inf',1,'numel',1);
+                        md = checkfield(md,'fieldname','calving.vel_threshold','NaN',1,'Inf',1,'numel',1);
                         md = checkfield(md,'fieldname','calving.vel_upperbound','NaN',1,'Inf',1,'numel',1);
                 end % }}}
@@ -79 +81 @@
                         fielddisplay(self,'yoffset','offset in y-axis');
                         fielddisplay(self,'vel_lowerbound','lowerbound of ice velocity to reduce the calving rate [m/a]');
+                        fielddisplay(self,'vel_threshold','threshold of ice velocity to reduce the calving rate [m/a]');
                         fielddisplay(self,'vel_upperbound','upperbound of ice velocity to reduce the calving rate [m/a]');
                 end % }}}
@@ -91 +94 @@
                         WriteData(fid,prefix,'object',self,'fieldname','yoffset','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','vel_lowerbound','format','Double','scale', 1./yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','vel_threshold','format','Double','scale', 1./yts);
                         WriteData(fid,prefix,'object',self,'fieldname','vel_upperbound','format','Double','scale', 1./yts);
                 end % }}}