Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26749) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 26750) @@ -24,5 +24,5 @@ int smb_model; - bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming; + bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming,isstochastic; /*Update elements: */ @@ -38,7 +38,12 @@ /*Figure out smb model: */ iomodel->FindConstant(&smb_model,"md.smb.model"); + iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing"); switch(smb_model){ case SMBforcingEnum: iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.); + if(isstochastic){ + iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",BaselineSmbMassBalanceEnum,0.); + } break; case SMBgembEnum: Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26749) +++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26750) @@ -58,5 +58,5 @@ multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed); for(int i=0;i(temparray); + xDelete(temparray); } ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -142,4 +142,25 @@ } break; + case SMBforcingEnum: + /*Delete SmbMassBalanceEnum at previous time step (required if it is transient)*/ + femmodel->inputs->DeleteInput(SmbMassBalanceEnum); + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + int numvertices = element->GetNumberOfVertices(); + IssmDouble baselinesmb; + IssmDouble smb_tot[numvertices]; + Input* baselinesmb_input = NULL; + baselinesmb_input = element->GetInput(BaselineSmbMassBalanceEnum); _assert_(baselinesmb_input); + element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum); + Gauss* gauss = element->NewGauss(); + for(int i=0;iGaussVertex(i); + baselinesmb_input->GetInputValue(&baselinesmb,gauss); + smb_tot[i] = baselinesmb+noisefield[dimensionid]; + } + element->AddInput(SmbMassBalanceEnum,&smb_tot[0],P1DGEnum); + delete gauss; + } + break; //VV working(29Nov2021) case FrictionWaterPressureEnum: Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26749) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26750) @@ -604,4 +604,5 @@ syn keyword cConstant BaselineCalvingCalvingrateEnum syn keyword cConstant BaselineFrictionEffectivePressureEnum +syn keyword cConstant BaselineSmbMassBalanceEnum syn keyword cConstant BedEnum syn keyword cConstant BedGRDEnum @@ -1582,4 +1583,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1588,6 +1590,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1602,6 +1604,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1610,9 +1612,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1630,4 +1633,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1640,4 +1644,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1652,6 +1657,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1667,11 +1672,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1692,8 +1697,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26749) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26750) @@ -600,4 +600,5 @@ BaselineCalvingCalvingrateEnum, BaselineFrictionEffectivePressureEnum, + BaselineSmbMassBalanceEnum, BedEnum, BedGRDEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26749) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26750) @@ -606,4 +606,5 @@ case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate"; case BaselineFrictionEffectivePressureEnum : return "BaselineFrictionEffectivePressure"; + case BaselineSmbMassBalanceEnum : return "BaselineSmbMassBalance"; case BedEnum : return "Bed"; case BedGRDEnum : return "BedGRD"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26749) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26750) @@ -597,4 +597,5 @@ syn keyword juliaConstC BaselineCalvingCalvingrateEnum syn keyword juliaConstC BaselineFrictionEffectivePressureEnum +syn keyword juliaConstC BaselineSmbMassBalanceEnum syn keyword juliaConstC BedEnum syn keyword juliaConstC BedGRDEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26749) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26750) @@ -618,4 +618,5 @@ else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum; else if (strcmp(name,"BaselineFrictionEffectivePressure")==0) return BaselineFrictionEffectivePressureEnum; + else if (strcmp(name,"BaselineSmbMassBalance")==0) return BaselineSmbMassBalanceEnum; else if (strcmp(name,"Bed")==0) return BedEnum; else if (strcmp(name,"BedGRD")==0) return BedGRDEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum; else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum; - else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; + if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; + else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum; else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum; - else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; + if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; + else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum; else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum; - else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum; + if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; + else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum; else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum; else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; - else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; + if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; + else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; - else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; + if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; + else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; - else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"Boundary")==0) return BoundaryEnum; + if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + else if (strcmp(name,"Boundary")==0) return BoundaryEnum; else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; - else if (strcmp(name,"Incremental")==0) return IncrementalEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"Indexed")==0) return IndexedEnum; + if (strcmp(name,"Incremental")==0) return IncrementalEnum; + else if (strcmp(name,"Indexed")==0) return IndexedEnum; else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"ElementInput")==0) return ElementInputEnum; @@ -1489,9 +1490,9 @@ else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else if (strcmp(name,"Regular")==0) return RegularEnum; - else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; + else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum; else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum; Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26749) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26750) @@ -48,5 +48,11 @@ for field=self.fields %Checking agreement of classes - if(contains(field,'SMB')) + if(contains(field,'SMBautoregression')) + mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); + if~(isequal(class(md.smb),char(mdname))) + error('md.smb does not agree with stochasticforcing field %s', char(field)); + end + end + if(contains(field,'SMBforcing')) mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); if~(isequal(class(md.smb),char(mdname))) @@ -156,5 +162,5 @@ %Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = {'DefaultCalving','FloatingMeltRate','SMBautoregression'}; %list of fields that need scaling *1/yts + scaledfields = {'DefaultCalving','FloatingMeltRate','SMBautoregression','SMBforcing'}; %list of fields that need scaling *1/yts tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable for i=1:num_fields @@ -199,5 +205,6 @@ 'FrictionWaterPressure',... 'FrontalForcingsRignotAutoregression',... - 'SMBautoregression' + 'SMBautoregression',... + 'SMBforcing' }; structure.mdnames = {... @@ -206,5 +213,6 @@ 'friction',... 'frontalforcingsrignotautoregression',... - 'SMBautoregression' + 'SMBautoregression',... + 'SMBforcing' }; end % }}} Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26749) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26750) @@ -67,5 +67,9 @@ for field in self.fields: # Checking agreement of classes - if 'SMB' in field: + if 'SMBautoregression' in field: + mdname = structstoch[field] + if (type(md.smb).__name__ != mdname): + raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(field)) + if 'SMBforcing' in field: mdname = structstoch[field] if (type(md.smb).__name__ != mdname): @@ -144,5 +148,5 @@ # Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = ['DefaultCalving','FloatingMeltRate','SMBautoregression'] # list of fields that need scaling * 1/yts + scaledfields = ['DefaultCalving','FloatingMeltRate','SMBautoregression','SMBforcing'] # list of fields that need scaling * 1/yts tempcovariance = np.copy(self.covariance) for i in range(num_fields): @@ -184,5 +188,6 @@ 'FrictionWaterPressure': 'friction', 'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression', - 'SMBautoregression': 'SMBautoregression'} + 'SMBautoregression': 'SMBautoregression', + 'SMBforcing': 'SMBforcing'} return structure # }}}