Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26483) @@ -3609,5 +3609,5 @@ IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; IssmDouble* phi_basin = xNew(arorder); - IssmDouble* smbspin = xNew(numvertices*arorder); + IssmDouble* smbspin = xNewZeroInit(numvertices*arorder); /*Get Basin ID and Basin coefficients*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26483) @@ -399,4 +399,5 @@ /*Add parameters that are not in standard nbvertices format*/ parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.randomflag",SmbRandomflagEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbAutoregressiveOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum)); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26483) @@ -171,7 +171,13 @@ my_rank=IssmComm::GetRank(); if(my_rank==0){ + bool randomflag{}; + femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum); + int fixedseed; for(int i=0;i(numbasins); - multivariateNormal(&temparray,numbasins,0.0,covmat); + /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/ + if(randomflag) fixedseed=-1; + else fixedseed = i; + multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed); for(int j=0;j(temparray); @@ -194,4 +200,5 @@ }/*}}}*/ void Smbautoregressionx(FemModel* femmodel){/*{{{*/ + /*Get time parameters*/ IssmDouble time,dt,starttime,tstep_ar; @@ -234,5 +241,12 @@ my_rank=IssmComm::GetRank(); if(my_rank==0){ - multivariateNormal(&noiseterms,numbasins,0.0,covmat); + bool randomflag{}; + femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum); + _printf_("randomflag: "<((time-starttime)/dt); + multivariateNormal(&noiseterms,numbasins,0.0,covmat,fixedseed); } ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26482) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26483) @@ -415,5 +415,11 @@ syn keyword cConstant SmbAccurefEnum syn keyword cConstant SmbAdThreshEnum +syn keyword cConstant SmbAutoregressionInitialTimeEnum +syn keyword cConstant SmbAutoregressionTimestepEnum +syn keyword cConstant SmbAutoregressiveOrderEnum syn keyword cConstant SmbAveragingEnum +syn keyword cConstant SmbBeta0Enum +syn keyword cConstant SmbBeta1Enum +syn keyword cConstant SmbCovmatEnum syn keyword cConstant SmbDesfacEnum syn keyword cConstant SmbDpermilEnum @@ -444,6 +450,9 @@ syn keyword cConstant SmbIsturbulentfluxEnum syn keyword cConstant SmbKEnum +syn keyword cConstant SmbNumBasinsEnum syn keyword cConstant SmbNumRequestedOutputsEnum syn keyword cConstant SmbPfacEnum +syn keyword cConstant SmbPhiEnum +syn keyword cConstant SmbRandomflagEnum syn keyword cConstant SmbRdlEnum syn keyword cConstant SmbRequestedOutputsEnum @@ -540,5 +549,9 @@ syn keyword cConstant AdjointpEnum syn keyword cConstant AdjointxEnum +syn keyword cConstant AdjointxBaseEnum +syn keyword cConstant AdjointxShearEnum syn keyword cConstant AdjointyEnum +syn keyword cConstant AdjointyBaseEnum +syn keyword cConstant AdjointyShearEnum syn keyword cConstant AdjointzEnum syn keyword cConstant AirEnum @@ -863,4 +876,5 @@ syn keyword cConstant SmbAdiffiniEnum syn keyword cConstant SmbAiniEnum +syn keyword cConstant SmbBasinsIdEnum syn keyword cConstant SmbBMaxEnum syn keyword cConstant SmbBMinEnum @@ -948,4 +962,5 @@ syn keyword cConstant SmbTmeanEnum syn keyword cConstant SmbTzEnum +syn keyword cConstant SmbValuesAutoregressionEnum syn keyword cConstant SmbVEnum syn keyword cConstant SmbVmeanEnum @@ -1411,4 +1426,5 @@ syn keyword cConstant SamplingSolutionEnum syn keyword cConstant SIAApproximationEnum +syn keyword cConstant SMBautoregressionEnum syn keyword cConstant SMBcomponentsEnum syn keyword cConstant SMBd18opddEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26482) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26483) @@ -448,4 +448,5 @@ SmbPfacEnum, SmbPhiEnum, + SmbRandomflagEnum, SmbRdlEnum, SmbRequestedOutputsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26483) @@ -456,4 +456,5 @@ case SmbPfacEnum : return "SmbPfac"; case SmbPhiEnum : return "SmbPhi"; + case SmbRandomflagEnum : return "SmbRandomflag"; case SmbRdlEnum : return "SmbRdl"; case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26483) @@ -465,4 +465,5 @@ else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum; + else if (strcmp(name,"SmbRandomflag")==0) return SmbRandomflagEnum; else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum; - else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; + if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; + else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum; else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum; else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum; - else if (strcmp(name,"DamageF")==0) return DamageFEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; + if (strcmp(name,"DamageF")==0) return DamageFEnum; + else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum; else if (strcmp(name,"DeltaIceThickness")==0) return DeltaIceThicknessEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum; - else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; + if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; + else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; - else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; + if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum; else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum; else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum; - else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; + if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; + else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; - else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; + if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; - else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; + if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; + else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"Dense")==0) return DenseEnum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"LoveLi")==0) return LoveLiEnum; else if (strcmp(name,"LoveLr")==0) return LoveLrEnum; - else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"MINI")==0) return MINIEnum; + if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum; + else if (strcmp(name,"MINI")==0) return MINIEnum; else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum; @@ -1489,9 +1490,9 @@ else if (strcmp(name,"Seg")==0) return SegEnum; else if (strcmp(name,"SegInput")==0) return SegInputEnum; - else if (strcmp(name,"Segment")==0) return SegmentEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; + if (strcmp(name,"Segment")==0) return SegmentEnum; + else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; else if (strcmp(name,"Separate")==0) return SeparateEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Random/random.cpp (revision 26482) +++ /issm/trunk-jpl/src/c/shared/Random/random.cpp (revision 26483) @@ -20,14 +20,17 @@ /*}}}*/ -void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev) { /*{{{*/ +void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seedfixed=-1) { /*{{{*/ + unsigned seed; /*Random seed using time_since_epoch*/ - unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); - std::default_random_engine generator(seed); + if(seedfixed==-1) seed = std::chrono::steady_clock::now().time_since_epoch().count(); + /*Seed fixed by input argument*/ + else seed = seedfixed; + std::default_random_engine generator(seed); /*Normal Probability Distribution*/ std::normal_distribution normdistri(mean,sdev); *prand = normdistri(generator); } /*}}}*/ -void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/ +void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/ IssmPDouble* sampleStandardNormal = xNew(dim); @@ -35,6 +38,7 @@ IssmDouble* Lchol = xNewZeroInit(dim*dim); - for(int i=0;i(Lchol); } /*}}}*/ -void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/ +void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/ IssmPDouble* sampleStandardNormal = xNew(dim); IssmDouble* sampleMultivariateNormal = xNew(dim); IssmDouble* Lchol = xNewZeroInit(dim*dim); - for(int i=0;i=',1,'numel',[1]); md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]); md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); - end % }}} + end function disp(self) % {{{ disp(sprintf(' surface forcings parameters:')); @@ -92,4 +95,5 @@ fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); fielddisplay(self,'covmat','inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]'); + fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); @@ -114,4 +118,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','randomflag','format','Boolean'); WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer'); WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');