Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27120) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27121) @@ -114,6 +114,4 @@ break; case CalvingParameterizationEnum: - iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum); - iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum); iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); break; @@ -221,10 +219,10 @@ break; case CalvingParameterizationEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.calving.scale_theta",CalvingScaleThetaEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.calving.amp_alpha",CalvingAmpAlphaEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.calving.midp",CalvingMidpointEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinearlaw",CalvingNonlinearLawEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.theta",CalvingThetaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.alpha",CalvingAlphaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.xoffset",CalvingXoffsetEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.yoffset",CalvingYoffsetEnum)); break; default: Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27120) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27121) @@ -376,4 +376,5 @@ IssmDouble time; IssmDouble coeff, indrate; + IssmDouble bed, bedrate = 1.0; /*Retrieve all inputs and parameters we will need*/ @@ -382,4 +383,5 @@ parameters->FindParam(&indrate,CalvingTestIndependentRateEnum,time); + Input *bs_input = this->GetInput(BedEnum);_assert_(bs_input); Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); @@ -399,9 +401,12 @@ lsf_slopey_input->GetInputValue(&dphidy,&gauss); + bs_input->GetInputValue(&bed,&gauss); + bedrate = (bed>0)?0.0:1.0; + vel=sqrt(vx*vx + vy*vy) + 1e-14; dphi=sqrt(dphidx*dphidx+dphidy*dphidy)+ 1e-14; - calvingratex[iv]= coeff*vx + indrate*dphidx/dphi; - calvingratey[iv]= coeff*vy + indrate*dphidy/dphi; + calvingratex[iv]= coeff*vx + bedrate*indrate*dphidx/dphi; + calvingratey[iv]= coeff*vy + bedrate*indrate*dphidy/dphi; calvingrate[iv] = sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); } @@ -835,10 +840,9 @@ IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; IssmDouble sigma_vm[NUMVERTICES]; - IssmDouble B,sigma_max,sigma_max_floating,sigma_max_grounded,n; + IssmDouble B, n; IssmDouble epse_2,groundedice,bed,sealevel; - IssmDouble arate, rho_ice, rho_water, thickness, paramX, Hab; - int use_parameter=0; - int nonlinear_law=0; - IssmDouble theta, alpha, midp, gamma; + IssmDouble arate, rho_ice, rho_water, thickness; + int use_parameter=-1; + IssmDouble gamma, theta, alpha, xoffset, yoffset; /* Get node coordinates and dof list: */ @@ -852,6 +856,4 @@ Input *bs_input = this->GetInput(BedEnum); _assert_(bs_input); Input *H_input = this->GetInput(ThicknessEnum); _assert_(H_input); - Input *smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input); - Input *smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input); Input *n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input *sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); @@ -864,8 +866,8 @@ /* Use which parameter */ this->FindParam(&use_parameter, CalvingUseParamEnum); - this->FindParam(&theta, CalvingScaleThetaEnum); - this->FindParam(&alpha, CalvingAmpAlphaEnum); - this->FindParam(&midp, CalvingMidpointEnum); - this->FindParam(&nonlinear_law, CalvingNonlinearLawEnum); + this->FindParam(&theta, CalvingThetaEnum); + this->FindParam(&alpha, CalvingAlphaEnum); + this->FindParam(&xoffset, CalvingXoffsetEnum); + this->FindParam(&yoffset, CalvingYoffsetEnum); /* Start looping on the number of vertices: */ @@ -881,7 +883,4 @@ gr_input->GetInputValue(&groundedice,gauss); bs_input->GetInputValue(&bed,gauss); - H_input->GetInputValue(&thickness,gauss); - smax_fl_input->GetInputValue(&sigma_max_floating,gauss); - smax_gr_input->GetInputValue(&sigma_max_grounded,gauss); vel=sqrt(vx*vx+vy*vy)+1.e-14; sl_input->GetInputValue(&sealevel,gauss); @@ -904,30 +903,16 @@ sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); - /*Tensile stress threshold*/ - if(groundedice<0) - sigma_max = sigma_max_floating; - else - sigma_max = sigma_max_grounded; - switch (use_parameter) { case 0: - /* bed elevation */ - paramX = bed; + /* 0 Linear: f(x) = y_{o} + \alpha (x+x_{o}) */ + gamma = yoffset = alpha * (bed+xoffset); break; case 1: - /* Height above floatation */ - if (bed>sealevel) paramX = 0.0; - else paramX = thickness - (rho_water/rho_ice) * (sealevel-bed); - break; - case 2: - /* Thicknese */ - paramX = thickness; - break; - case 4: - /* bed elevation+linear curve fitting */ - paramX = bed; + /* 1 tanh: f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */ + gamma = yoffset - 0.5*theta*tanh(alpha*(bed+xoffset)); break; case -1: - /* use nothing, just the arate*/ + /* nothing, just the arate*/ + gamma = 1; break; default: @@ -936,29 +921,9 @@ /* Compute the hyperbolic tangent function */ - if ((use_parameter>-0.5) & (use_parameter<3)) { - gamma = 0.5*theta*(1.0-tanh((paramX+midp)/alpha))+(1.0-theta); - if (gamma<0.0) gamma =0.0; - } - else if (use_parameter>=3) { - gamma = alpha*paramX + theta; - if (gamma > 1.0) gamma = 1.0; - if (gamma < 0.0) gamma = 0.0; - } - else gamma = 1; + if (gamma > 1.0) gamma = 1.0; + if (gamma < 0.0) gamma = 0.0; /*-------------------------------------------*/ - if (nonlinear_law) { - /*This von Mises type has too strong positive feedback with vel included - * calvingrate[iv] = (arate+sigma_vm[iv]*vel/sigma_max)*gamma; - */ - Hab = thickness - (rho_water/rho_ice) * (sealevel-bed); - if (Hab < 0.) Hab = 0.; - if (bed > sealevel) Hab = 0.; - - calvingrate[iv] = (arate+Hab/sigma_max)*gamma; - } - else { - calvingrate[iv] = arate*gamma; - } + calvingrate[iv] = arate*gamma; } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27120) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27121) @@ -109,8 +109,8 @@ CalvingTestIndependentRateEnum, CalvingUseParamEnum, - CalvingScaleThetaEnum, - CalvingAmpAlphaEnum, - CalvingMidpointEnum, - CalvingNonlinearLawEnum, + CalvingThetaEnum, + CalvingAlphaEnum, + CalvingXoffsetEnum, + CalvingYoffsetEnum, ConfigurationTypeEnum, ConstantsGEnum, Index: /issm/trunk-jpl/src/m/classes/calvingparameterization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingparameterization.m (revision 27120) +++ /issm/trunk-jpl/src/m/classes/calvingparameterization.m (revision 27121) @@ -6,12 +6,10 @@ classdef calvingparameterization properties (SetAccess=public) - stress_threshold_groundedice = 0.; - stress_threshold_floatingice = 0.; min_thickness = 0.; use_param = 0; - scale_theta = 0.; - amp_alpha = 0; - midp = 0; - nonlinearlaw = 0; + theta = 0.; + alpha = 0; + xoffset = 0; + yoffset = 0; end methods @@ -37,24 +35,20 @@ end % }}} function self = setdefaultparameters(self) % {{{ - - %Default sigma max - self.stress_threshold_groundedice = 1e6; - self.stress_threshold_floatingice = 150e3; - %For now we turn this off by setting the threshold to 0 self.min_thickness = 0.; %parameters for the spatial temporal seperation - %The coefficient follows: \gamma = \frac{\theta}{2}(1-\tanh(\frac{b+p}{\alpha}))+(1-\theta) - % 0 - Use bed elevation, 1 - use heigh above floatation + %The coefficient follows: gamma= f(x) + % 0 - f(x) = y_{o} + \alpha (x+x_{o}) + % 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) self.use_param = 0; - % between 0 and 1, larger theta means more reduction for shallower ice - self.scale_theta = 0; - % alpha in the denominator - self.amp_alpha = 0; - % mid-point of this step function - self.midp = 0; - % if use a nonlinear calving law - self.nonlinearlaw = 0; + % the amplifier + self.theta = 0; + % the slope alpha + self.alpha = 0; + % offset in x-axis + self.xoffset = 0; + % offset in y-axis + self.yoffset = 0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -62,36 +56,29 @@ if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end - md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1); - md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1); md = checkfield(md,'fieldname','calving.min_thickness','>=',0,'NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1, 2, 3, 4]); - md = checkfield(md,'fieldname','calving.scale_theta','NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','calving.amp_alpha','<>',0,'NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','calving.midp','NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','calving.nonlinearlaw','values',[0, 1]); + md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1]); + md = checkfield(md,'fieldname','calving.theta','NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.alpha','NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.xoffset','NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.yoffset','NaN',1,'Inf',1,'numel',1); end % }}} function disp(self) % {{{ disp(sprintf(' Calving test parameters:')); - fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]'); - fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]'); fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]'); - fielddisplay(self,'use_param','0 - Use bed elevation, 1 - use heigh above floatation, 2 - use ice thickness'); - fielddisplay(self,'scale_theta','larger than 0, larger theta means more reduction for shallower ice'); - fielddisplay(self,'amp_alpha','alpha'); - fielddisplay(self,'midp','mid-point'); - fielddisplay(self,'nonlinearlaw','use a nonlinear law'); - + fielddisplay(self,'use_param','-1 - just use frontal ablation rate, 0 - f(x) = y_{o} + \alpha (x+x_{o}), 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o}))'); + fielddisplay(self,'theta','the amplifier'); + fielddisplay(self,'alpha','the slope'); + fielddisplay(self,'xoffset','offset in x-axis'); + fielddisplay(self,'yoffset','offset in y-axis'); end % }}} function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; WriteData(fid,prefix,'name','md.calving.law','data',9,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','use_param','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','scale_theta','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','amp_alpha','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','midp','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','nonlinearlaw','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','theta','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','xoffset','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','yoffset','format','Double'); end % }}} end