Changeset 27121


Ignore:
Timestamp:
06/29/22 08:39:32 (3 years ago)
Author:
Cheng Gong
Message:

CHG: rewrite calving parameterization with tanh and linear functions

Location:
issm/trunk-jpl/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r27102 r27121  
    114114                        break;
    115115                case CalvingParameterizationEnum:
    116                         iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
    117                         iomodel->FetchDataToInput(inputs,elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
    118116                        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
    119117                        break;
    … …  
    221219                        break;
    222220                case CalvingParameterizationEnum:
     221                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum));
    223222                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum));
    224                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.use_param",CalvingUseParamEnum));
    225                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.scale_theta",CalvingScaleThetaEnum));
    226                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.amp_alpha",CalvingAmpAlphaEnum));
    227                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.midp",CalvingMidpointEnum));
    228                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.nonlinearlaw",CalvingNonlinearLawEnum));
     223                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.theta",CalvingThetaEnum));
     224                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.alpha",CalvingAlphaEnum));
     225                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.xoffset",CalvingXoffsetEnum));
     226                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.yoffset",CalvingYoffsetEnum));
    229227                        break;
    230228                default:

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

    r27102 r27121  
    376376        IssmDouble  time;
    377377        IssmDouble  coeff, indrate;
     378        IssmDouble  bed, bedrate = 1.0;
    378379
    379380        /*Retrieve all inputs and parameters we will need*/
    … …  
    382383        parameters->FindParam(&indrate,CalvingTestIndependentRateEnum,time);
    383384
     385        Input *bs_input = this->GetInput(BedEnum);_assert_(bs_input);
    384386        Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
    385387        Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
    … …  
    399401      lsf_slopey_input->GetInputValue(&dphidy,&gauss);
    400402
     403                bs_input->GetInputValue(&bed,&gauss);
     404                bedrate = (bed>0)?0.0:1.0;
     405
    401406      vel=sqrt(vx*vx + vy*vy) + 1e-14;
    402407      dphi=sqrt(dphidx*dphidx+dphidy*dphidy)+ 1e-14;
    403408
    404                 calvingratex[iv]= coeff*vx + indrate*dphidx/dphi;
    405                 calvingratey[iv]= coeff*vy + indrate*dphidy/dphi;
     409                calvingratex[iv]= coeff*vx + bedrate*indrate*dphidx/dphi;
     410                calvingratey[iv]= coeff*vy + bedrate*indrate*dphidy/dphi;
    406411                calvingrate[iv] = sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    407412        }
    … …  
    835840        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
    836841        IssmDouble  sigma_vm[NUMVERTICES];
    837         IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
     842        IssmDouble  B, n;
    838843        IssmDouble  epse_2,groundedice,bed,sealevel;
    839         IssmDouble  arate, rho_ice, rho_water, thickness, paramX, Hab;
    840         int                     use_parameter=0;
    841         int                     nonlinear_law=0;
    842         IssmDouble  theta, alpha, midp, gamma;
     844        IssmDouble  arate, rho_ice, rho_water, thickness;
     845        int                     use_parameter=-1;
     846        IssmDouble  gamma, theta, alpha, xoffset, yoffset;
    843847
    844848        /* Get node coordinates and dof list: */
    … …  
    852856        Input *bs_input      = this->GetInput(BedEnum);                               _assert_(bs_input);
    853857        Input *H_input       = this->GetInput(ThicknessEnum);                         _assert_(H_input);
    854         Input *smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
    855         Input *smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
    856858        Input *n_input       = this->GetInput(MaterialsRheologyNEnum);                _assert_(n_input);
    857859        Input *sl_input      = this->GetInput(SealevelEnum);                          _assert_(sl_input);
    … …  
    864866        /* Use which parameter  */
    865867        this->FindParam(&use_parameter, CalvingUseParamEnum);
    866         this->FindParam(&theta, CalvingScaleThetaEnum);
    867         this->FindParam(&alpha, CalvingAmpAlphaEnum);
    868         this->FindParam(&midp, CalvingMidpointEnum);
    869         this->FindParam(&nonlinear_law, CalvingNonlinearLawEnum);
     868        this->FindParam(&theta, CalvingThetaEnum);
     869        this->FindParam(&alpha, CalvingAlphaEnum);
     870        this->FindParam(&xoffset, CalvingXoffsetEnum);
     871        this->FindParam(&yoffset, CalvingYoffsetEnum);
    870872
    871873        /* Start looping on the number of vertices: */
    … …  
    881883                gr_input->GetInputValue(&groundedice,gauss);
    882884                bs_input->GetInputValue(&bed,gauss);
    883                 H_input->GetInputValue(&thickness,gauss);
    884                 smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
    885                 smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
    886885                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    887886                sl_input->GetInputValue(&sealevel,gauss);
    … …  
    904903                sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
    905904
    906                 /*Tensile stress threshold*/
    907                 if(groundedice<0)
    908                  sigma_max = sigma_max_floating;
    909                 else
    910                  sigma_max = sigma_max_grounded;
    911 
    912905                switch (use_parameter) {
    913906                        case 0:
    914                                 /* bed elevation */
    915                                 paramX = bed;
     907                                /* 0 Linear: f(x) = y_{o} + \alpha (x+x_{o}) */
     908                                gamma = yoffset = alpha * (bed+xoffset);
    916909                                break;
    917910                        case 1:
    918                                 /* Height above floatation */
    919                                 if (bed>sealevel)       paramX = 0.0;
    920                                 else paramX = thickness - (rho_water/rho_ice) * (sealevel-bed);
    921                                 break;
    922                         case 2:
    923                                 /* Thicknese */
    924                                 paramX = thickness;
    925                                 break;
    926                         case 4:
    927                                 /* bed elevation+linear curve fitting */
    928                                 paramX = bed;
     911                                /* 1 tanh: f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
     912                                gamma = yoffset -  0.5*theta*tanh(alpha*(bed+xoffset));
    929913                                break;
    930914                        case -1:
    931                                 /* use nothing, just the arate*/
     915                                /* nothing, just the arate*/
     916                                gamma = 1;
    932917                                break;
    933918                        default:
    … …  
    936921
    937922                /* Compute the hyperbolic tangent function */
    938                 if ((use_parameter>-0.5) & (use_parameter<3)) {
    939                         gamma = 0.5*theta*(1.0-tanh((paramX+midp)/alpha))+(1.0-theta);
    940                         if (gamma<0.0) gamma =0.0;
    941                 }
    942                 else if (use_parameter>=3) {
    943                         gamma = alpha*paramX + theta;
    944                         if (gamma > 1.0) gamma = 1.0;
    945                         if (gamma < 0.0) gamma = 0.0;
    946                 }
    947                 else gamma = 1;
     923                if (gamma > 1.0) gamma = 1.0;
     924                if (gamma < 0.0) gamma = 0.0;
    948925
    949926                /*-------------------------------------------*/
    950                 if (nonlinear_law) {
    951                         /*This von Mises type has too strong positive feedback with vel included
    952                          * calvingrate[iv] = (arate+sigma_vm[iv]*vel/sigma_max)*gamma;
    953                          */
    954                         Hab = thickness - (rho_water/rho_ice) * (sealevel-bed);
    955                         if (Hab < 0.) Hab = 0.;
    956                         if (bed > sealevel) Hab = 0.;
    957 
    958                         calvingrate[iv] = (arate+Hab/sigma_max)*gamma;
    959                 }
    960                 else {
    961                         calvingrate[iv] = arate*gamma;
    962                 }
     927                calvingrate[iv] = arate*gamma;
    963928        }
    964929

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

    r27031 r27121  
    109109        CalvingTestIndependentRateEnum,
    110110        CalvingUseParamEnum,
    111         CalvingScaleThetaEnum,
    112         CalvingAmpAlphaEnum,
    113         CalvingMidpointEnum,
    114         CalvingNonlinearLawEnum,
     111        CalvingThetaEnum,
     112        CalvingAlphaEnum,
     113        CalvingXoffsetEnum,
     114        CalvingYoffsetEnum,
    115115        ConfigurationTypeEnum,
    116116        ConstantsGEnum,

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

    r27009 r27121  
    66classdef calvingparameterization
    77        properties (SetAccess=public)
    8                 stress_threshold_groundedice = 0.;
    9                 stress_threshold_floatingice = 0.;
    108                min_thickness = 0.;
    119                use_param = 0;
    12                 scale_theta = 0.;
    13                 amp_alpha = 0;
    14                 midp = 0;
    15                 nonlinearlaw = 0;
     10                theta = 0.;
     11                alpha = 0;
     12                xoffset = 0;
     13                yoffset = 0;
    1614        end
    1715        methods
    … …  
    3735                end % }}}
    3836                function self = setdefaultparameters(self) % {{{
    39 
    40                         %Default sigma max
    41                         self.stress_threshold_groundedice = 1e6;
    42                         self.stress_threshold_floatingice = 150e3;
    43 
    4437                        %For now we turn this off by setting the threshold to 0
    4538                        self.min_thickness = 0.;
    4639
    4740                        %parameters for the spatial temporal seperation
    48                         %The coefficient follows: \gamma = \frac{\theta}{2}(1-\tanh(\frac{b+p}{\alpha}))+(1-\theta)
    49                         % 0 - Use bed elevation, 1 - use heigh above floatation
     41                        %The coefficient follows: gamma= f(x)
     42                        % 0 - f(x) = y_{o} + \alpha (x+x_{o})
     43                        % 1 - f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o}))
    5044                        self.use_param = 0;
    51                         % between 0 and 1, larger theta means more reduction for shallower ice
    52                         self.scale_theta = 0;
    53                         % alpha in the denominator
    54                         self.amp_alpha = 0;
    55                         % mid-point of this step function
    56                         self.midp = 0;
    57                         % if use a nonlinear calving law
    58                         self.nonlinearlaw = 0;
     45                        % the amplifier
     46                        self.theta = 0;
     47                        % the slope alpha
     48                        self.alpha = 0;
     49                        % offset in x-axis
     50                        self.xoffset = 0;
     51                        % offset in y-axis
     52                        self.yoffset = 0;
    5953                end % }}}
    6054                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    6256                        if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end
    6357
    64                         md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1);
    65                         md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1);
    6658                        md = checkfield(md,'fieldname','calving.min_thickness','>=',0,'NaN',1,'Inf',1,'numel',1);
    67                         md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1, 2, 3, 4]);
    68                         md = checkfield(md,'fieldname','calving.scale_theta','NaN',1,'Inf',1,'numel',1);
    69                         md = checkfield(md,'fieldname','calving.amp_alpha','<>',0,'NaN',1,'Inf',1,'numel',1);
    70                         md = checkfield(md,'fieldname','calving.midp','NaN',1,'Inf',1,'numel',1);
    71                         md = checkfield(md,'fieldname','calving.nonlinearlaw','values',[0, 1]);
     59                        md = checkfield(md,'fieldname','calving.use_param','values',[-1, 0, 1]);
     60                        md = checkfield(md,'fieldname','calving.theta','NaN',1,'Inf',1,'numel',1);
     61                        md = checkfield(md,'fieldname','calving.alpha','NaN',1,'Inf',1,'numel',1);
     62                        md = checkfield(md,'fieldname','calving.xoffset','NaN',1,'Inf',1,'numel',1);
     63                        md = checkfield(md,'fieldname','calving.yoffset','NaN',1,'Inf',1,'numel',1);
    7264                end % }}}
    7365                function disp(self) % {{{
    7466                        disp(sprintf('   Calving test parameters:'));
    75                         fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]');
    76                         fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]');
    7767                        fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]');
    78                         fielddisplay(self,'use_param','0 - Use bed elevation, 1 - use heigh above floatation, 2 - use ice thickness');
    79                         fielddisplay(self,'scale_theta','larger than 0, larger theta means more reduction for shallower ice');
    80                         fielddisplay(self,'amp_alpha','alpha');
    81                         fielddisplay(self,'midp','mid-point');
    82                         fielddisplay(self,'nonlinearlaw','use a nonlinear law');
    83 
     68                        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}))');
     69                        fielddisplay(self,'theta','the amplifier');
     70                        fielddisplay(self,'alpha','the slope');
     71                        fielddisplay(self,'xoffset','offset in x-axis');
     72                        fielddisplay(self,'yoffset','offset in y-axis');
    8473                end % }}}
    8574                function marshall(self,prefix,md,fid) % {{{
    8675                        yts=md.constants.yts;
    8776                        WriteData(fid,prefix,'name','md.calving.law','data',9,'format','Integer');
    88                         WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1);
    89                         WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1);
    9077                        WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double');
    9178                        WriteData(fid,prefix,'object',self,'fieldname','use_param','format','Integer');
    92                         WriteData(fid,prefix,'object',self,'fieldname','scale_theta','format','Double');
    93                         WriteData(fid,prefix,'object',self,'fieldname','amp_alpha','format','Double');
    94                         WriteData(fid,prefix,'object',self,'fieldname','midp','format','Double');
    95                         WriteData(fid,prefix,'object',self,'fieldname','nonlinearlaw','format','Integer');
     79                        WriteData(fid,prefix,'object',self,'fieldname','theta','format','Double');
     80                        WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Double');
     81                        WriteData(fid,prefix,'object',self,'fieldname','xoffset','format','Double');
     82                        WriteData(fid,prefix,'object',self,'fieldname','yoffset','format','Double');
    9683                end % }}}
    9784        end
Note: See TracChangeset for help on using the changeset viewer.