Changeset 26487


Ignore:
Timestamp:
10/19/21 06:50:29 (3 years ago)
Author:
vverjans
Message:

CHG: adding Nightly runs for frontalforcingsrignot and SMBautoregression, changed frontalforcingrignot basin_id variable to be element-specific

Location:
issm/trunk-jpl
Files:
4 added
10 edited

Legend:

Unmodified
Added
Removed

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

    r26468 r26487  
    9393                        break;
    9494                case FrontalForcingsRignotEnum:
    95                         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin",FrontalForcingsBasinIdEnum);
     95                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum);
    9696                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
    9797                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
    … …  
    141141                        break;
    142142                case FrontalForcingsRignotEnum:
    143                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.numberofbasins",FrontalForcingsNumberofBasinsEnum));
     143                        parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum));
    144144                        break;
    145145                default:

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

    r26483 r26487  
    33783378
    33793379} /*}}}*/
     3380void       Element::RignotMeltParameterization(){/*{{{*/
     3381
     3382        const int numvertices = this->GetNumberOfVertices();
     3383   IssmDouble A, B, alpha, beta;
     3384   IssmDouble bed,qsg,qsg_basin,TF,yts;
     3385   int numbasins,basinid;
     3386   IssmDouble* basin_icefront_area=NULL;
     3387
     3388   /* Coefficients */
     3389   A    = 3e-4;
     3390   B    = 0.15;
     3391   alpha = 0.39;
     3392   beta = 1.18;
     3393
     3394   /*Get inputs*/
     3395        this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
     3396        Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
     3397   Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);     _assert_(qsg_input);
     3398   Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
     3399
     3400        this->FindParam(&yts, ConstantsYtsEnum);
     3401   this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
     3402   this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
     3403   IssmDouble meltrates[numvertices]; 
     3404
     3405   /* Start looping on the number of vertices: */
     3406   GaussTria gauss;
     3407   for(int iv=0;iv<numvertices;iv++){
     3408      gauss.GaussVertex(iv);
     3409
     3410      /* Get variables */
     3411      bed_input->GetInputValue(&bed,&gauss);
     3412      qsg_input->GetInputValue(&qsg,&gauss);
     3413      TF_input->GetInputValue(&TF,&gauss);
     3414
     3415      if(basin_icefront_area[basinid]==0.) meltrates[iv]=0.;
     3416      else{
     3417         /* change the unit of qsg (m^3/d -> m/d) with ice front area */
     3418         qsg_basin=qsg/basin_icefront_area[basinid];
     3419
     3420         /* calculate melt rates */
     3421         meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
     3422      }
     3423
     3424      if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
     3425      if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
     3426   }
     3427
     3428   /*Add input*/
     3429   this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum);
     3430
     3431   /*Cleanup and return*/
     3432   xDelete<IssmDouble>(basin_icefront_area);
     3433}
     3434/*}}}*/
    33803435void       Element::SetBoolInput(Inputs* inputs,int enum_in,bool value){/*{{{*/
    33813436

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r26478 r26487  
    167167                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    168168                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
     169                void               RignotMeltParameterization();
    169170                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);
    170171                void               ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum);
    … …  
    345346                virtual void       ResetFSBasalBoundaryCondition()=0;
    346347                virtual void       ResetHooks()=0;
    347                 virtual void       RignotMeltParameterization(void){_error_("not implemented yet");};
    348348                virtual void       SetElementInput(int enum_in,IssmDouble value){_error_("not implemented yet");};
    349349                virtual void       SetElementInput(int enum_in,IssmDouble value, int type)=0;

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

    r26468 r26487  
    34023402        if(this->hneighbors) this->hneighbors->reset();
    34033403
    3404 }
    3405 /*}}}*/
    3406 void       Penta::RignotMeltParameterization(){/*{{{*/
    3407 
    3408         if(!this->IsOnBase()) return;
    3409 
    3410         IssmDouble A, B, alpha, beta;
    3411         IssmDouble bed,qsg,qsg_basin,TF,yts;
    3412         int numbasins;
    3413         IssmDouble basinid[NUMVERTICES];
    3414         IssmDouble* basin_icefront_area=NULL;
    3415 
    3416         /* Coefficients */
    3417         A    = 3e-4;
    3418         B    = 0.15;
    3419         alpha = 0.39;
    3420         beta = 1.18;
    3421 
    3422         /*Get inputs*/
    3423         Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
    3424         Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);               _assert_(qsg_input);
    3425         Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
    3426         Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
    3427 
    3428         this->FindParam(&yts, ConstantsYtsEnum);
    3429         this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    3430         this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
    3431 
    3432         IssmDouble meltrates[NUMVERTICES2D];  //frontal melt-rate
    3433 
    3434         /* Start looping on the number of vertices: */
    3435         GaussPenta gauss;
    3436         for(int iv=0;iv<NUMVERTICES2D;iv++){
    3437                 gauss.GaussVertex(iv);
    3438 
    3439                 /* Get variables */
    3440                 bed_input->GetInputValue(&bed,&gauss);
    3441                 qsg_input->GetInputValue(&qsg,&gauss);
    3442                 TF_input->GetInputValue(&TF,&gauss);
    3443 
    3444                 if(basinid[iv]==0 || basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
    3445                 else{
    3446                         /* change the unit of qsg (m^3/d -> m/d) with ice front area */
    3447                         qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
    3448 
    3449                         /* calculate melt rates */
    3450                         meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
    3451                 }
    3452 
    3453                 if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
    3454                 if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
    3455         }
    3456 
    3457         /*Add input*/
    3458         this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum);
    3459 
    3460         this->InputExtrude(CalvingMeltingrateEnum,-1);
    3461 
    3462         /*Cleanup and return*/
    3463         xDelete<IssmDouble>(basin_icefront_area);
    34643404}
    34653405/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r26296 r26487  
    171171                void           ResetFSBasalBoundaryCondition(void);
    172172                void           ResetHooks();
    173                 void                            RignotMeltParameterization();
    174173                void           SetElementInput(int enum_in,IssmDouble value);
    175174                void           SetElementInput(int enum_in,IssmDouble value,int type){_error_("not implemented yet");};

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

    r26468 r26487  
    47174717        if(this->hneighbors) this->hneighbors->reset();
    47184718
    4719 }
    4720 /*}}}*/
    4721 void       Tria::RignotMeltParameterization(){/*{{{*/
    4722 
    4723    IssmDouble A, B, alpha, beta;
    4724         IssmDouble bed,qsg,qsg_basin,TF,yts;
    4725         int numbasins;
    4726         IssmDouble basinid[NUMVERTICES];
    4727         IssmDouble* basin_icefront_area=NULL;
    4728 
    4729         /* Coefficients */
    4730         A    = 3e-4;
    4731         B    = 0.15;
    4732         alpha = 0.39;
    4733         beta = 1.18;
    4734 
    4735         /*Get inputs*/
    4736         Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
    4737         Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);               _assert_(qsg_input);
    4738         Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
    4739         Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
    4740 
    4741         this->FindParam(&yts, ConstantsYtsEnum);
    4742         this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    4743         this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
    4744 
    4745         IssmDouble meltrates[NUMVERTICES];  //frontal melt-rate
    4746 
    4747         /* Start looping on the number of vertices: */
    4748         GaussTria gauss;
    4749         for(int iv=0;iv<NUMVERTICES;iv++){
    4750                 gauss.GaussVertex(iv);
    4751 
    4752                 /* Get variables */
    4753                 bed_input->GetInputValue(&bed,&gauss);
    4754                 qsg_input->GetInputValue(&qsg,&gauss);
    4755                 TF_input->GetInputValue(&TF,&gauss);
    4756 
    4757                 if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
    4758                 else{
    4759                         /* change the unit of qsg (m^3/d -> m/d) with ice front area */
    4760                         qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
    4761 
    4762                         /* calculate melt rates */
    4763                         meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
    4764                 }
    4765 
    4766                 if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
    4767                 if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
    4768         }
    4769 
    4770         /*Add input*/
    4771         this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum);
    4772 
    4773         /*Cleanup and return*/
    4774         xDelete<IssmDouble>(basin_icefront_area);
    47754719}
    47764720/*}}}*/

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

    r26296 r26487  
    132132                void        ResetHooks();
    133133                void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
    134                 void        RignotMeltParameterization();
    135134                void        SetElementInput(int enum_in,IssmDouble value);
    136135                void        SetElementInput(int enum_in,IssmDouble value,int type);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r26468 r26487  
    16021602void FemModel::IcefrontAreax(){/*{{{*/
    16031603
    1604         int numvertices      = 6;
    1605         int numbasins;
    1606         IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
     1604        int numbasins,BasinId;
    16071605        this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    1608         IssmDouble* basin_icefront_area           = xNewZeroInit<IssmDouble>(numbasins);
    1609 
    1610         for(int basin=1;basin<numbasins+1;basin++){
     1606        IssmDouble* basin_icefront_area = xNewZeroInit<IssmDouble>(numbasins);
     1607
     1608        for(int basin=0;basin<numbasins;basin++){
    16111609                IssmDouble local_icefront_area = 0;
    16121610                IssmDouble total_icefront_area;
    16131611
     1612                /*Add contribution of each element to icefront area of the basin*/
    16141613                for(Object* & object : this->elements->objects){
    16151614                        Element* element= xDynamicCast<Element*>(object);
    1616                         element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
    1617                         for(int j=0;j<numvertices;j++){
    1618                                 if(BasinId[j]==basin){
    1619                                         local_icefront_area+=element->GetIcefrontArea();
    1620                                         break;
    1621                                 }
    1622                         }
     1615                        element->GetInputValue(&BasinId,FrontalForcingsBasinIdEnum);
     1616                        if(BasinId==basin) local_icefront_area+=element->GetIcefrontArea();
    16231617                }
    16241618                ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    16251619                ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    16261620
    1627                 basin_icefront_area[basin-1]=total_icefront_area;
    1628         }
    1629 
     1621                basin_icefront_area[basin]=total_icefront_area;
     1622        }
     1623       
    16301624        this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
    16311625
    16321626        xDelete<IssmDouble>(basin_icefront_area);
    1633         xDelete<IssmDouble>(BasinId);
    16341627}/*}}}*/
    16351628void FemModel::IcefrontMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r26485 r26487  
    243243                bool randomflag{};
    244244                femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
    245                 _printf_("randomflag: "<<randomflag<<'\n');
    246245                /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
    247246                int fixedseed;

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

    r23673 r26487  
    66classdef frontalforcingsrignot
    77        properties (SetAccess=public)
    8                 basin= NaN;
    9                 numberofbasins=0;
    10                 subglacial_discharge= NaN;
    11                 thermalforcing=NaN;
     8                basin_id             = NaN;
     9                num_basins           = 0;
     10                subglacial_discharge = NaN;
     11                thermalforcing       = NaN;
    1212        end
    1313        methods
    … …  
    3535                function self = setdefaultparameters(self) % {{{
    3636
    37                         basin=NaN;
    38                         numberofbasins=0;
    39                         subglacial_discharge=NaN;
    40                         thermalforcing=NaN;
     37                        basin_id             = NaN;
     38                        num_basins           = 0;
     39                        subglacial_discharge = NaN;
     40                        thermalforcing       = NaN;
    4141                end % }}}
    4242                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    4444                        if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end
    4545
    46                         md = checkfield(md,'fieldname','frontalforcings.basin','>',0,'nan',1,'Inf',1);
    47                         md = checkfield(md,'fieldname','frontalforcings.numberofbasins','numel',[1]);
     46                        md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0);
     47                        md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]);
    4848                        md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'nan',1,'Inf',1,'timeseries',1);
    4949                        md = checkfield(md,'fieldname','frontalforcings.thermalforcing','nan',1,'Inf',1,'timeseries',1);
    … …  
    5252                function disp(self) % {{{
    5353                        disp(sprintf('   Frontalforcings parameters:'));
    54                         fielddisplay(self,'basin','basin ID for vertices');
    55                         fielddisplay(self,'numberofbasins', 'number of basins');
     54                        fielddisplay(self,'basin_id','basin ID for elements');
     55                        fielddisplay(self,'num_basins', 'number of basins');
    5656                        fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]');
    5757                        fielddisplay(self,'thermalforcing','thermal forcing [∘C]');
    … …  
    5959                function marshall(self,prefix,md,fid) % {{{
    6060                        WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',2,'format','Integer');
    61                         WriteData(fid,prefix,'object',self,'fieldname','basin','format','DoubleMat','mattype',1);
    62                         WriteData(fid,prefix,'object',self,'fieldname','numberofbasins','format','Integer');
     61                        WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.frontalforcings.basin_id','format','IntMat','mattype',2); %0-indexed
     62                        WriteData(fid,prefix,'object',self,'fieldname','num_basins','format','Integer');
    6363                        WriteData(fid,prefix,'object',self,'fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    6464                        WriteData(fid,prefix,'object',self,'fieldname','thermalforcing','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
Note: See TracChangeset for help on using the changeset viewer.