Changeset 26483


Ignore:
Timestamp:
10/15/21 12:36:16 (3 years ago)
Author:
vverjans
Message:

NEW: randomflag allows to use fixed random seed in random.cpp, zero-initialization of smbspin for SmbAutoregression

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

Legend:

Unmodified
Added
Removed

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

    r26482 r26483  
    36093609        IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; 
    36103610        IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
    3611    IssmDouble* smbspin     = xNew<IssmDouble>(numvertices*arorder);
     3611   IssmDouble* smbspin     = xNewZeroInit<IssmDouble>(numvertices*arorder);
    36123612
    36133613        /*Get Basin ID and Basin coefficients*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r26477 r26483  
    399399         /*Add parameters that are not in standard nbvertices format*/
    400400         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum));
     401         parameters->AddObject(iomodel->CopyConstantObject("md.smb.randomflag",SmbRandomflagEnum));
    401402         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbAutoregressiveOrderEnum));
    402403         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum));

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

    r26482 r26483  
    171171        my_rank=IssmComm::GetRank();
    172172        if(my_rank==0){
     173                bool randomflag{};
     174                femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
     175                int fixedseed;
    173176                for(int i=0;i<nspin;i++){
    174177                        IssmDouble* temparray = xNew<IssmDouble>(numbasins);
    175                         multivariateNormal(&temparray,numbasins,0.0,covmat);
     178                        /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
     179                        if(randomflag) fixedseed=-1;
     180                        else fixedseed = i;
     181                        multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
    176182                        for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
    177183                        xDelete<IssmDouble>(temparray);
    … …  
    194200}/*}}}*/
    195201void Smbautoregressionx(FemModel* femmodel){/*{{{*/
     202
    196203        /*Get time parameters*/
    197204        IssmDouble time,dt,starttime,tstep_ar;
    … …  
    234241        my_rank=IssmComm::GetRank();
    235242        if(my_rank==0){
    236                 multivariateNormal(&noiseterms,numbasins,0.0,covmat);
     243                bool randomflag{};
     244                femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
     245                _printf_("randomflag: "<<randomflag<<'\n');
     246                /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
     247                int fixedseed;
     248                if(randomflag) fixedseed=-1;
     249                else fixedseed = static_cast<int>((time-starttime)/dt);
     250                multivariateNormal(&noiseterms,numbasins,0.0,covmat,fixedseed);
    237251        }
    238252        ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26445 r26483  
    415415syn keyword cConstant SmbAccurefEnum
    416416syn keyword cConstant SmbAdThreshEnum
     417syn keyword cConstant SmbAutoregressionInitialTimeEnum
     418syn keyword cConstant SmbAutoregressionTimestepEnum
     419syn keyword cConstant SmbAutoregressiveOrderEnum
    417420syn keyword cConstant SmbAveragingEnum
     421syn keyword cConstant SmbBeta0Enum
     422syn keyword cConstant SmbBeta1Enum
     423syn keyword cConstant SmbCovmatEnum
    418424syn keyword cConstant SmbDesfacEnum
    419425syn keyword cConstant SmbDpermilEnum
    … …  
    444450syn keyword cConstant SmbIsturbulentfluxEnum
    445451syn keyword cConstant SmbKEnum
     452syn keyword cConstant SmbNumBasinsEnum
    446453syn keyword cConstant SmbNumRequestedOutputsEnum
    447454syn keyword cConstant SmbPfacEnum
     455syn keyword cConstant SmbPhiEnum
     456syn keyword cConstant SmbRandomflagEnum
    448457syn keyword cConstant SmbRdlEnum
    449458syn keyword cConstant SmbRequestedOutputsEnum
    … …  
    540549syn keyword cConstant AdjointpEnum
    541550syn keyword cConstant AdjointxEnum
     551syn keyword cConstant AdjointxBaseEnum
     552syn keyword cConstant AdjointxShearEnum
    542553syn keyword cConstant AdjointyEnum
     554syn keyword cConstant AdjointyBaseEnum
     555syn keyword cConstant AdjointyShearEnum
    543556syn keyword cConstant AdjointzEnum
    544557syn keyword cConstant AirEnum
    … …  
    863876syn keyword cConstant SmbAdiffiniEnum
    864877syn keyword cConstant SmbAiniEnum
     878syn keyword cConstant SmbBasinsIdEnum
    865879syn keyword cConstant SmbBMaxEnum
    866880syn keyword cConstant SmbBMinEnum
    … …  
    948962syn keyword cConstant SmbTmeanEnum
    949963syn keyword cConstant SmbTzEnum
     964syn keyword cConstant SmbValuesAutoregressionEnum
    950965syn keyword cConstant SmbVEnum
    951966syn keyword cConstant SmbVmeanEnum
    … …  
    14111426syn keyword cConstant SamplingSolutionEnum
    14121427syn keyword cConstant SIAApproximationEnum
     1428syn keyword cConstant SMBautoregressionEnum
    14131429syn keyword cConstant SMBcomponentsEnum
    14141430syn keyword cConstant SMBd18opddEnum

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

    r26477 r26483  
    448448        SmbPfacEnum,
    449449        SmbPhiEnum,
     450        SmbRandomflagEnum,
    450451        SmbRdlEnum,
    451452        SmbRequestedOutputsEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26477 r26483  
    456456                case SmbPfacEnum : return "SmbPfac";
    457457                case SmbPhiEnum : return "SmbPhi";
     458                case SmbRandomflagEnum : return "SmbRandomflag";
    458459                case SmbRdlEnum : return "SmbRdl";
    459460                case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r26477 r26483  
    465465              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    466466              else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
     467              else if (strcmp(name,"SmbRandomflag")==0) return SmbRandomflagEnum;
    467468              else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
    468469              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    … …  
    505506              else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
    506507              else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
    507               else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     511              if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
     512              else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    512513              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    513514              else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
    … …  
    628629              else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
    629630              else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum;
    630               else if (strcmp(name,"DamageF")==0) return DamageFEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
     634              if (strcmp(name,"DamageF")==0) return DamageFEnum;
     635              else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
    635636              else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
    636637              else if (strcmp(name,"DeltaIceThickness")==0) return DeltaIceThicknessEnum;
    … …  
    751752              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    752753              else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum;
    753               else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
     757              if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
     758              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
    758759              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
    759760              else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum;
    … …  
    874875              else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
    875876              else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
    876               else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
     880              if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     881              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    881882              else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
    882883              else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum;
    … …  
    997998              else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum;
    998999              else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum;
    999               else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
     1003              if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
     1004              else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
    10041005              else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    10051006              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    … …  
    11201121              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
    11211122              else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    1122               else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
     1126              if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
     1127              else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
    11271128              else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
    11281129              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    … …  
    12431244              else if (strcmp(name,"DataSet")==0) return DataSetEnum;
    12441245              else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
    1245               else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
     1249              if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
     1250              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    12501251              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    12511252              else if (strcmp(name,"Dense")==0) return DenseEnum;
    … …  
    13661367              else if (strcmp(name,"LoveLi")==0) return LoveLiEnum;
    13671368              else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
    1368               else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"MINI")==0) return MINIEnum;
     1372              if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
     1373              else if (strcmp(name,"MINI")==0) return MINIEnum;
    13731374              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    13741375              else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
    … …  
    14891490              else if (strcmp(name,"Seg")==0) return SegEnum;
    14901491              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    1491               else if (strcmp(name,"Segment")==0) return SegmentEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     1495              if (strcmp(name,"Segment")==0) return SegmentEnum;
     1496              else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    14961497              else if (strcmp(name,"Separate")==0) return SeparateEnum;
    14971498              else if (strcmp(name,"Seq")==0) return SeqEnum;

  • issm/trunk-jpl/src/c/shared/Random/random.cpp

    r26482 r26483  
    2020/*}}}*/
    2121
    22 void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev) { /*{{{*/
     22void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seedfixed=-1) { /*{{{*/
    2323
     24        unsigned seed;
    2425        /*Random seed using time_since_epoch*/
    25         unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count();
    26    std::default_random_engine generator(seed);
     26        if(seedfixed==-1) seed = std::chrono::steady_clock::now().time_since_epoch().count();
     27        /*Seed fixed by input argument*/
     28        else seed = seedfixed;
     29        std::default_random_engine generator(seed);
    2730        /*Normal Probability Distribution*/
    2831   std::normal_distribution<IssmPDouble> normdistri(mean,sdev);
    2932        *prand = normdistri(generator);
    3033} /*}}}*/
    31 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
     34void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
    3235   
    3336        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    … …  
    3538   IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
    3639
    37    for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
    38    CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
     40        /*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
     41        for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed);
     42        CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
    3943   
    4044        /*Matrix by vector multiplication*/
    … …  
    5155   xDelete<IssmDouble>(Lchol);
    5256} /*}}}*/
    53 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
     57void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
    5458       
    5559        IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    5660        IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
    5761        IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
    58         for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
    5962
     63        /*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
     64        for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed);
    6065        CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
    6166

  • issm/trunk-jpl/src/c/shared/Random/random.h

    r26477 r26483  
    66#define _RANDOM_H_
    77
    8 void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev);
    9 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix);
    10 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix);
     8void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev, int seedfixed=-1);
     9void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1);
     10void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1);
    1111
    1212#endif //ifndef _RANDOM_H_

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

    r26477 r26483  
    1414                phi               = NaN;
    1515                covmat            = NaN;       
    16                 basin_id          = NaN;       
     16                basin_id          = NaN;
     17                randomflag        = 1;
    1718                steps_per_step    = 1;
    1819                averaging         = 0;
    … …  
    6263                end % }}}
    6364                function self = setdefaultparameters(self) % {{{
    64                         self.ar_order    = 0.0; %autoregression model of order 0
     65                        self.ar_order    = 0.0; %autoregression model of order 0
     66                        self.randomflag  = 1;
    6567                end % }}}
    66                 function md = checkconsistency(self,md,solution,analyses) % {{{
     68                function md = checkconsistency(self,md,solution,analyses)
    6769
    6870                        if ismember('MasstransportAnalysis',analyses),
    … …  
    7678                                md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
    7779                                md = checkfield(md,'fieldname','smb.covmat','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.num_basins]);
     80                                md = checkfield(md,'fieldname','smb.randomflag','numel',[1],'values',[0 1]);
    7881                        end
    7982                        md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
    8083                        md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
    8184                        md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
    82                 end % }}}
     85                end
    8386                function disp(self) % {{{
    8487                        disp(sprintf('   surface forcings parameters:'));
    … …  
    9295                        fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
    9396                        fielddisplay(self,'covmat','inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]');
     97                        fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    9498                        fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
    9599                        fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
    … …  
    114118                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',md.constants.yts);
    115119                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',md.constants.yts);
     120                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','randomflag','format','Boolean');
    116121                        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
    117122                        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
Note: See TracChangeset for help on using the changeset viewer.