Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26615) @@ -83,10 +83,16 @@ /*Get moving front parameters*/ - int calvinglaw; - iomodel->FindConstant(&calvinglaw,"md.calving.law"); - switch(calvinglaw){ - case DefaultCalvingEnum: - iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum); - break; + bool isstochastic; + int calvinglaw; + iomodel->FindConstant(&calvinglaw,"md.calving.law"); + iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing"); + switch(calvinglaw){ + case DefaultCalvingEnum: + iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum); + if(isstochastic){ + iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum); + } + break; case CalvingLevermannEnum: iomodel->FetchDataToInput(inputs,elements,"md.calving.coeff",CalvinglevermannCoeffEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26615) @@ -99,24 +99,47 @@ xDelete(phi_basin); }/*}}}*/ -void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); - int basinid,M,N; - IssmDouble beta0_basin,beta1_basin,noise_basin; +void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/ + + const int numvertices = this->GetNumberOfVertices(); + int basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type; + IssmDouble beta0_basin,beta1_basin,noiseterm; IssmDouble* phi_basin = xNew(arorder); IssmDouble* varlist = xNew(numvertices); IssmDouble* valuesautoregression = NULL; - - /*Get Basin ID and Basin coefficients*/ - if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + Input* noiseterm_input = NULL; + + /*Get field-specific enums*/ + switch(enum_type){ + case(SMBautoregressionEnum): + arenum_type = SmbValuesAutoregressionEnum; + basinenum_type = SmbBasinsIdEnum; + noiseenum_type = SmbAutoregressionNoiseEnum; + outenum_type = SmbMassBalanceEnum; + break; + case(FrontalForcingsRignotAutoregressionEnum): + arenum_type = ThermalforcingValuesAutoregressionEnum; + basinenum_type = FrontalForcingsBasinIdEnum; + noiseenum_type = ThermalforcingAutoregressionNoiseEnum; + outenum_type = FrontalForcingsThermalForcingEnum; + break; + } + + /*Get noise and autoregressive terms*/ + this->GetInputValue(&basinid,basinenum_type); + if(isfieldstochastic){ + noiseterm_input = this->GetInput(noiseenum_type); + Gauss* gauss = this->NewGauss(); + noiseterm_input->GetInputValue(&noiseterm,gauss); + delete gauss; + } + else noiseterm = 0.0; + this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M); + + /*Get basin coefficients*/ for(int ii=0;iiinputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M); - - /*If not AR model timestep: take the old values of variable*/ + + /*If not AR model timestep: take the old values of variable*/ if(isstepforar==false){ for(int i=0;i(numvertices*arorder); /*Assign newest values and shift older values*/ for(int i=0;iinputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); + this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder); xDelete(temparray); } /*Add input to element*/ - if(enum_type==SMBautoregressionEnum) this->AddInput(SmbMassBalanceEnum,varlist,P1Enum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum); + this->AddInput(outenum_type,varlist,P1Enum); /*Cleanup*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26614) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26615) @@ -69,5 +69,5 @@ bool AnyFSet(void); void AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type); - void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type); + void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type); void ComputeLambdaS(void); void ComputeNewDamage(); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26615) @@ -80,5 +80,6 @@ /*Load parameters*/ bool isstochastic; - int M,N,Nphi,arorder,numbasins,my_rank; + bool istfstochastic = false; + int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); @@ -87,5 +88,4 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* noiseterms = xNewZeroInit(numbasins); femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); @@ -94,5 +94,4 @@ femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); if(isstochastic){ @@ -102,8 +101,7 @@ femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum); _assert_(M==numbasins); - } + if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true; } + xDelete(stochasticfields); } /*Time elapsed with respect to AR model initial time*/ @@ -113,5 +111,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum); } @@ -120,4 +118,3 @@ xDelete(beta1); xDelete(phi); - xDelete(noiseterms); }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26615) @@ -500,28 +500,28 @@ bool isstochasticforcing; - parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); - if(isstochasticforcing){ - int num_fields; - char** fields; - parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum)); - iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields"); - if(num_fields<1) _error_("no stochasticforcing fields found"); - int* stochasticforcing_enums = xNew(num_fields); - for(int i=0;i(fields[i]); - } - xDelete(fields); - parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields)); - xDelete(stochasticforcing_enums); - + parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); + if(isstochasticforcing){ + int num_fields,stochastic_dim; + char** fields; + parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum)); + iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields"); + if(num_fields<1) _error_("no stochasticforcing fields found"); + int* stochasticforcing_enums = xNew(num_fields); + for(int i=0;i(fields[i]); + } + xDelete(fields); + parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields)); + xDelete(stochasticforcing_enums); parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions"); - parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance"); - parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N)); - xDelete(transparam); - } + iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions"); + parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance"); + parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N)); + xDelete(transparam); + } /*Deal with mass flux segments: {{{*/ Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26615) @@ -36,5 +36,5 @@ if(randomflag) fixedseed=-1; else fixedseed = reCast((time-starttime)/dt); - /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ + /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ IssmDouble* temparray = NULL; multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed); @@ -46,15 +46,66 @@ int i=0; for(int j=0;j(dimensions[j]); for(int k=0;kparameters->SetParam(noisefield,dimensions[j],enum_type); - i=i+dimensions[j]; + + int dimensionid; + + /*Deal with the autoregressive models*/ + if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){ + switch(fields[j]){ + case SMBautoregressionEnum: + dimenum_type = SmbBasinsIdEnum; + noiseenum_type = SmbAutoregressionNoiseEnum; + break; + case FrontalForcingsRignotAutoregressionEnum: + dimenum_type = FrontalForcingsBasinIdEnum; + noiseenum_type = ThermalforcingAutoregressionNoiseEnum; + break; + } + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + int numvertices = element->GetNumberOfVertices(); + IssmDouble* noise_element = xNew(numvertices); + element->GetInputValue(&dimensionid,dimenum_type); + for(int i=0;iAddInput(noiseenum_type,noise_element,P0Enum); + xDelete(noise_element); + } + } + else{ + switch(fields[j]){ + case SMBautoregressionEnum: + case FrontalForcingsRignotAutoregressionEnum: + /*Already done above*/ + break; + case DefaultCalvingEnum: + /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/ + femmodel->inputs->DeleteInput(CalvingCalvingrateEnum); + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + int numvertices = element->GetNumberOfVertices(); + IssmDouble baselinecalvingrate; + IssmDouble calvingrate_tot[numvertices]; + Input* baselinecalvingrate_input = NULL; + baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input); + element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum); + Gauss* gauss = element->NewGauss(); + for(int i=0;iGaussVertex(i); + baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss); + calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]); + } + element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum); + delete gauss; + } + break; + default: + _error_("Field "<(noisefield); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26615) @@ -177,66 +177,62 @@ void Smbautoregressionx(FemModel* femmodel){/*{{{*/ - /*Get time parameters*/ - IssmDouble time,dt,starttime,tstep_ar; - femmodel->parameters->FindParam(&time,TimeEnum); - femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + /*Get time parameters*/ + IssmDouble time,dt,starttime,tstep_ar; + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum); - /*Initialize module at first time step*/ - if(time<=starttime+dt){SmbautoregressionInitx(femmodel);} - /*Determine if this is a time step for the AR model*/ - bool isstepforar = false; - - #ifndef _HAVE_AD_ + /*Initialize module at first time step*/ + if(time<=starttime+dt){SmbautoregressionInitx(femmodel);} + /*Determine if this is a time step for the AR model*/ + bool isstepforar = false; + + #ifndef _HAVE_AD_ if((fmod(time,tstep_ar)parameters->FindParam(&numbasins,SmbNumBasinsEnum); - femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); - IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; - - femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - - /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ - IssmDouble* noiseterms = xNewZeroInit(numbasins); - femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); + #else + _error_("not implemented yet"); + #endif + + /*Load parameters*/ + bool isstochastic; + bool issmbstochastic = false; + int M,N,Nphi,arorder,numbasins,my_rank; + femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); + femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); + IssmDouble tinit_ar; + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + + femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); + femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + + femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); if(isstochastic){ - int numstochasticfields; + int numstochasticfields; int* stochasticfields; femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum); femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); for(int i=0;iparameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum); _assert_(M==numbasins); - } - } - xDelete(stochasticfields); - } - /*Time elapsed with respect to AR model initial time*/ - IssmDouble telapsed_ar = time-tinit_ar; - - /*Loop over each element to compute SMB at vertices*/ - for(Object* &object:femmodel->elements->objects){ - Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum); - } - - /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); - xDelete(noiseterms); + if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true; + } + xDelete(stochasticfields); + } + /*Time elapsed with respect to AR model initial time*/ + IssmDouble telapsed_ar = time-tinit_ar; + + /*Loop over each element to compute SMB at vertices*/ + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); + } + + /*Cleanup*/ + xDelete(beta0); + xDelete(beta1); + xDelete(phi); }/*}}}*/ void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26614) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26615) @@ -398,4 +398,5 @@ syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum syn keyword cConstant StochasticForcingCovarianceEnum +syn keyword cConstant StochasticForcingDefaultDimensionEnum syn keyword cConstant StochasticForcingDimensionsEnum syn keyword cConstant StochasticForcingFieldsEnum @@ -428,5 +429,4 @@ syn keyword cConstant SmbAdThreshEnum syn keyword cConstant SmbAutoregressionInitialTimeEnum -syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbAutoregressionTimestepEnum syn keyword cConstant SmbAutoregressiveOrderEnum @@ -504,5 +504,4 @@ syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum syn keyword cConstant StressbalanceShelfDampeningEnum -syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalIsdrainicecolumnEnum syn keyword cConstant ThermalIsdynamicbasalspcEnum @@ -602,4 +601,5 @@ syn keyword cConstant BaseSlopeXEnum syn keyword cConstant BaseSlopeYEnum +syn keyword cConstant BaselineCalvingCalvingrateEnum syn keyword cConstant BedEnum syn keyword cConstant BedGRDEnum @@ -895,4 +895,5 @@ syn keyword cConstant SmbAdiffiniEnum syn keyword cConstant SmbAiniEnum +syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbBasinsIdEnum syn keyword cConstant SmbBMaxEnum @@ -997,4 +998,5 @@ syn keyword cConstant SolidearthExternalDisplacementUpRateEnum syn keyword cConstant SolidearthExternalGeoidRateEnum +syn keyword cConstant StochasticForcingDefaultIdEnum syn keyword cConstant StrainRateeffectiveEnum syn keyword cConstant StrainRateparallelEnum @@ -1032,4 +1034,5 @@ syn keyword cConstant TemperaturePicardEnum syn keyword cConstant TemperatureSEMICEnum +syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalforcingValuesAutoregressionEnum syn keyword cConstant ThermalSpctemperatureEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26614) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26615) @@ -392,4 +392,5 @@ SolidearthSettingsOceanAreaScalingEnum, StochasticForcingCovarianceEnum, + StochasticForcingDefaultDimensionEnum, StochasticForcingDimensionsEnum, StochasticForcingFieldsEnum, @@ -422,5 +423,4 @@ SmbAdThreshEnum, SmbAutoregressionInitialTimeEnum, - SmbAutoregressionNoiseEnum, SmbAutoregressionTimestepEnum, SmbAutoregressiveOrderEnum, @@ -498,5 +498,4 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, - ThermalforcingAutoregressionNoiseEnum, ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, @@ -598,4 +597,5 @@ BaseSlopeXEnum, BaseSlopeYEnum, + BaselineCalvingCalvingrateEnum, BedEnum, BedGRDEnum, @@ -891,4 +891,5 @@ SmbAdiffiniEnum, SmbAiniEnum, + SmbAutoregressionNoiseEnum, SmbBasinsIdEnum, SmbBMaxEnum, @@ -994,4 +995,5 @@ SolidearthExternalDisplacementUpRateEnum, SolidearthExternalGeoidRateEnum, + StochasticForcingDefaultIdEnum, StrainRateeffectiveEnum, StrainRateparallelEnum, @@ -1029,4 +1031,5 @@ TemperaturePicardEnum, TemperatureSEMICEnum, + ThermalforcingAutoregressionNoiseEnum, ThermalforcingValuesAutoregressionEnum, ThermalSpctemperatureEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26615) @@ -400,4 +400,5 @@ case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling"; case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance"; + case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension"; case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions"; case StochasticForcingFieldsEnum : return "StochasticForcingFields"; @@ -430,5 +431,4 @@ case SmbAdThreshEnum : return "SmbAdThresh"; case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime"; - case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep"; case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder"; @@ -506,5 +506,4 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; - case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; @@ -604,4 +603,5 @@ case BaseSlopeXEnum : return "BaseSlopeX"; case BaseSlopeYEnum : return "BaseSlopeY"; + case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate"; case BedEnum : return "Bed"; case BedGRDEnum : return "BedGRD"; @@ -897,4 +897,5 @@ case SmbAdiffiniEnum : return "SmbAdiffini"; case SmbAiniEnum : return "SmbAini"; + case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbBasinsIdEnum : return "SmbBasinsId"; case SmbBMaxEnum : return "SmbBMax"; @@ -999,4 +1000,5 @@ case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate"; case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate"; + case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId"; case StrainRateeffectiveEnum : return "StrainRateeffective"; case StrainRateparallelEnum : return "StrainRateparallel"; @@ -1034,4 +1036,5 @@ case TemperaturePicardEnum : return "TemperaturePicard"; case TemperatureSEMICEnum : return "TemperatureSEMIC"; + case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression"; case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26614) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26615) @@ -170,6 +170,12 @@ syn keyword juliaConstC FrictionVoidRatioEnum syn keyword juliaConstC FrontalForcingsBasinIcefrontAreaEnum +syn keyword juliaConstC FrontalForcingsAutoregressionInitialTimeEnum +syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum +syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum +syn keyword juliaConstC FrontalForcingsBeta0Enum +syn keyword juliaConstC FrontalForcingsBeta1Enum syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum syn keyword juliaConstC FrontalForcingsParamEnum +syn keyword juliaConstC FrontalForcingsPhiEnum syn keyword juliaConstC GrdModelEnum syn keyword juliaConstC GroundinglineFrictionInterpolationEnum @@ -384,4 +390,11 @@ syn keyword juliaConstC SolidearthSettingsGrdOceanEnum syn keyword juliaConstC SolidearthSettingsOceanAreaScalingEnum +syn keyword juliaConstC StochasticForcingCovarianceEnum +syn keyword juliaConstC StochasticForcingDefaultDimensionEnum +syn keyword juliaConstC StochasticForcingDimensionsEnum +syn keyword juliaConstC StochasticForcingFieldsEnum +syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum +syn keyword juliaConstC StochasticForcingNumFieldsEnum +syn keyword juliaConstC StochasticForcingRandomflagEnum syn keyword juliaConstC RotationalPolarMoiEnum syn keyword juliaConstC SolidearthSettingsReltolEnum @@ -414,5 +427,4 @@ syn keyword juliaConstC SmbBeta0Enum syn keyword juliaConstC SmbBeta1Enum -syn keyword juliaConstC SmbCovmatEnum syn keyword juliaConstC SmbDesfacEnum syn keyword juliaConstC SmbDpermilEnum @@ -424,4 +436,5 @@ syn keyword juliaConstC SmbDtEnum syn keyword juliaConstC SmbEnum +syn keyword juliaConstC SmbEIdxEnum syn keyword juliaConstC SmbFEnum syn keyword juliaConstC SmbInitDensityScalingEnum @@ -432,4 +445,5 @@ syn keyword juliaConstC SmbIsdelta18oEnum syn keyword juliaConstC SmbIsdensificationEnum +syn keyword juliaConstC SmbIsdeltaLWupEnum syn keyword juliaConstC SmbIsfirnwarmingEnum syn keyword juliaConstC SmbIsgraingrowthEnum @@ -447,5 +461,4 @@ syn keyword juliaConstC SmbPfacEnum syn keyword juliaConstC SmbPhiEnum -syn keyword juliaConstC SmbRandomflagEnum syn keyword juliaConstC SmbRdlEnum syn keyword juliaConstC SmbRequestedOutputsEnum @@ -460,4 +473,5 @@ syn keyword juliaConstC SmbT0dryEnum syn keyword juliaConstC SmbT0wetEnum +syn keyword juliaConstC SmbTeThreshEnum syn keyword juliaConstC SmbTdiffEnum syn keyword juliaConstC SmbThermoDeltaTScalingEnum @@ -580,4 +594,5 @@ syn keyword juliaConstC BaseSlopeXEnum syn keyword juliaConstC BaseSlopeYEnum +syn keyword juliaConstC BaselineCalvingCalvingrateEnum syn keyword juliaConstC BedEnum syn keyword juliaConstC BedGRDEnum @@ -873,4 +888,5 @@ syn keyword juliaConstC SmbAdiffiniEnum syn keyword juliaConstC SmbAiniEnum +syn keyword juliaConstC SmbAutoregressionNoiseEnum syn keyword juliaConstC SmbBasinsIdEnum syn keyword juliaConstC SmbBMaxEnum @@ -894,4 +910,5 @@ syn keyword juliaConstC SmbDiniEnum syn keyword juliaConstC SmbDlwrfEnum +syn keyword juliaConstC SmbDulwrfValueEnum syn keyword juliaConstC SmbDswrfEnum syn keyword juliaConstC SmbDswdiffrfEnum @@ -974,4 +991,5 @@ syn keyword juliaConstC SolidearthExternalDisplacementUpRateEnum syn keyword juliaConstC SolidearthExternalGeoidRateEnum +syn keyword juliaConstC StochasticForcingDefaultIdEnum syn keyword juliaConstC StrainRateeffectiveEnum syn keyword juliaConstC StrainRateparallelEnum @@ -1009,4 +1027,6 @@ syn keyword juliaConstC TemperaturePicardEnum syn keyword juliaConstC TemperatureSEMICEnum +syn keyword juliaConstC ThermalforcingAutoregressionNoiseEnum +syn keyword juliaConstC ThermalforcingValuesAutoregressionEnum syn keyword juliaConstC ThermalSpctemperatureEnum syn keyword juliaConstC ThicknessAbsGradientEnum @@ -1255,4 +1275,5 @@ syn keyword juliaConstC FrontalForcingsDefaultEnum syn keyword juliaConstC FrontalForcingsRignotEnum +syn keyword juliaConstC FrontalForcingsRignotAutoregressionEnum syn keyword juliaConstC FsetEnum syn keyword juliaConstC FullMeltOnPartiallyFloatingEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26614) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26615) @@ -409,4 +409,5 @@ else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum; else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum; + else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum; else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum; else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum; @@ -439,5 +440,4 @@ else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum; - else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum; else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum; else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum; @@ -518,5 +518,4 @@ else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum; else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum; - else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum; else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum; else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; @@ -616,4 +615,5 @@ else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum; else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum; + else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum; else if (strcmp(name,"Bed")==0) return BedEnum; else if (strcmp(name,"BedGRD")==0) return BedGRDEnum; @@ -918,4 +918,5 @@ else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum; else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; + else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum; else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum; else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbT")==0) return SmbTEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; - else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; + if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; + else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; @@ -1023,4 +1024,5 @@ else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum; else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; + else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; @@ -1058,4 +1060,5 @@ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum; + else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum; else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum; else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; @@ -1118,11 +1121,11 @@ else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; - else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; - else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; - else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; + if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; + else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; + else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; + else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum; else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; @@ -1241,11 +1244,11 @@ else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; - else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; - else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum; - else if (strcmp(name,"Channel")==0) return ChannelEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; + if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; + else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum; + else if (strcmp(name,"Channel")==0) return ChannelEnum; + else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum; else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; @@ -1364,11 +1367,11 @@ else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; else if (strcmp(name,"Inputs")==0) return InputsEnum; - else if (strcmp(name,"Internal")==0) return InternalEnum; - else if (strcmp(name,"Intersect")==0) return IntersectEnum; - else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"J")==0) return JEnum; + if (strcmp(name,"Internal")==0) return InternalEnum; + else if (strcmp(name,"Intersect")==0) return IntersectEnum; + else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; + else if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum; else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum; @@ -1487,11 +1490,11 @@ else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; - else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; + else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26614) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26615) @@ -8,5 +8,6 @@ isstochasticforcing = 0; fields = NaN; - dimensions = NaN; + defaultdimension = 0; + default_id = NaN; covariance = NaN; randomflag = 1; @@ -34,25 +35,59 @@ num_fields = numel(self.fields); - size_tot = sum(self.dimensions); + + %Check that covariance matrix is positive definite + try + chol(self.covariance); + catch + error('md.stochasticforcing.covariance is not positive definite'); + end + + %Check that all fields agree with the corresponding md class and if any field needs the default params + checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality + for field=self.fields + %Checking agreement of classes + if(contains(field,'SMB')) + if~(isequal(class(md.smb),char(field))) + error('md.smb does not agree with stochasticforcing field %s', char(field)); + end + end + if(contains(field,'frontalforcings')) + if~(isequal(class(md.frontalforcings),char(field))) + error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); + end + end + %Checking for specific dimensions + if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression')) + checkdefaults = true; %field with non-specific dimensionality + end + end + + %Retrieve sum of all the field dimensionalities + size_tot = self.defaultdimension*num_fields; + indSMBar = -1; %about to check for index of SMBautoregression + indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression + if any(contains(self.fields,'SMBautoregression')) + size_tot = size_tot-self.defaultdimension+md.smb.num_basins; + indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021) + end + if any(contains(self.fields,'FrontalForcingsRignotAutoregression')) + size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins; + indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021) + end + + if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency + if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar))))~=0) + error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance'); + end + end md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]); - md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021) - md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field + md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); - - %Check that all fields agree with the corresponding md class - for field=self.fields - if(contains(field,'SMB')) - if~(isequal(class(md.smb),char(field))) - error('md.smb does not agree with stochasticforcing field %s', char(field)); - end - end - if(contains(field,'frontalforcings')) - if~(isequal(class(md.frontalforcings),char(field))) - error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); - end - end - end + if(checkdefaults) %need to check the defaults + md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]); + end end % }}} function disp(self) % {{{ @@ -60,5 +95,6 @@ fielddisplay(self,'isstochasticforcing','is stochasticity activated?'); fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}'); - fielddisplay(self,'dimensions','dimensionality of each field'); + fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)'); + fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'); fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)'); fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); @@ -77,10 +113,25 @@ return else + + %Retrieve dimensionality of each field + dimensions = self.defaultdimension*ones(1,num_fields); + ind = 1; + for field=self.fields + %Checking for specific dimensions + if(strcmp(field,'SMBautoregression')) + dimensions(ind) = md.smb.num_basins; + end + if(strcmp(field,'FrontalForcingsRignotAutoregression')) + dimensions(ind) = md.frontalforcings.num_basins; + end + ind = ind+1; + end + %Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts + scaledfields = {'DefaultCalving','SMBautoregression'}; %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 if any(strcmp(scaledfields,self.fields(i))) - inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))]; + inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))]; for row=inds %scale rows corresponding to scaled field tempcovariance(row,:) = 1./yts*tempcovariance(row,:); @@ -93,5 +144,7 @@ WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray'); - WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat'); + WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat'); + WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed + WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer'); WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat'); WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean'); @@ -105,6 +158,7 @@ list = {... - 'SMBautoregression',... - 'FrontalForcingsRignotAutoregression' + 'DefaultCalving',... + 'FrontalForcingsRignotAutoregression',... + 'SMBautoregression' }; end % }}} Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26614) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26615) @@ -15,8 +15,9 @@ def __init__(self, *args): # {{{ self.isstochasticforcing = 0 - self.fields = np.nan - self.dimensions = np.nan - self.covariance = np.nan - self.randomflag = 1 + self.fields = np.nan + self.defaultdimension = 0 + self.default_id = np.nan + self.covariance = np.nan + self.randomflag = 1 if len(args) == 0: @@ -29,7 +30,10 @@ s += '{}\n'.format(fielddisplay(self, 'isstochasticforcing', 'is stochasticity activated?')) s += '{}\n'.format(fielddisplay(self, 'fields', 'fields with stochasticity applied, ex: [\'SMBautoregression\'], or [\'FrontalForcingsRignotAutoregression\']')) + s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)')) + s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)')) s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)')) s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)')) s += 'Available fields:\n' + s += ' DefaultCalving\n' s += ' SMBautoregression\n' s += ' FrontalForcingsRignotAutoregression (thermal forcing)\n' @@ -51,13 +55,13 @@ num_fields = numel(self.fields) - size_tot = np.sum(self.dimensions) - md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) # VV check here 'cell' (19Oct2021) - md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list - md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix - md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) + #Check that covariance matrix is positive definite + try: + np.linalg.cholesky(self.covariance); + except: + error('md.stochasticforcing.covariance is not positive definite'); # Check that all fields agree with the corresponding md class + checkdefaults = False for field in self.fields: if (contains(field, 'SMB')): @@ -67,4 +71,30 @@ if not (type(md.frontalforcings) == field): error('md.frontalforcings does not agree with stochasticforcing field {}'.format(field)) + #Checking for specific dimensions + if (field!='SMBautoregression' or field!='FrontalForcingsRignotAutoregression'): + checkdefaults = True #field with non-specific dimensionality + + #Retrieve sum of all the field dimensionalities + size_tot = self.defaultdimension*num_fields; + indSMBar = -1; #about to check for index of SMBautoregression + indTFar = -1; #about to check for index of FrontalForcingsRignotAutoregression + if ('SMBautoregression' in self.fields): + size_tot = size_tot-self.defaultdimension+md.smb.num_basins + indSMBar = np.where(self.fields=='SMBautoregression')[0][0] + if ('FrontalForcingsRignotAutoregression' in self.fields): + size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins + indSMBar = np.where(self.fields=='FrontalForcingsRignotAutoregression')[0][0] + if (indSMBar!=-1 and indTFar!=-1): + if((md.smb.ar_timestep!=md.frontalforcings.ar_timestep) and np.any(self.covariance[np.sum(self.dimensions[0:indSMBar]).astype(int):np.sum(self.dimensions[0:indSMBar+1]).astype(int),np.sum(self.dimensions[0:indTFar]).astype(int):np.sum(self.dimensions[0:indTFar+1]).astype(int)]!=0)): + error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance') + + md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) + #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list + md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix + md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) + if (checkdefaults): + md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements,1]) return md # }}} @@ -83,17 +113,28 @@ return md else: + #Retrieve dimensionality of each field + dimensions = self.defaultdimension*np.ones(num_fields) + for ind,field in enumerate(self.fields): + #Checking for specific dimensions + if (field=='SMBautoregression'): + dimensions[ind] = md.smb.num_basins + if (field=='FrontalForcingsRignotAutoregression'): + dimensions[ind] = md.frontalforcings.num_basins + # Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts - tempcovariance = self.covariance + scaledfields = ['DefaultCalving','SMBautoregression'] # list of fields that need scaling * 1/yts + tempcovariance = np.copy(self.covariance) for i in range(num_fields): if self.fields[i] in scaledfields: inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1]))) for row in inds: # scale rows corresponding to scaled field - tempcovariance[row, :] = 1 / yts * self.covariance[row, :] + tempcovariance[row, :] = 1 / yts * tempcovariance[row, :] for col in inds: # scale columns corresponding to scaled field - tempcovariance[:, col] = 1 / yts * self.covariance[:, col] + tempcovariance[:, col] = 1 / yts * tempcovariance[:, col] WriteData(fid, prefix, 'data', num_fields, 'name', 'md.stochasticforcing.num_fields', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray') - WriteData(fid, prefix, 'object', self, 'fieldname','dimensions', 'format', 'IntMat') + WriteData(fid, prefix, 'data', 'dimensions', 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat') + WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed! + WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer') WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean') @@ -104,5 +145,6 @@ """ return [ - 'SMBautoregression', - 'FrontalForcingsRignotAutoregression' + 'DefaultCalving', + 'FrontalForcingsRignotAutoregression', + 'SMBautoregression' ] Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26614) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26615) @@ -49,5 +49,4 @@ md.stochasticforcing.isstochasticforcing = 1; md.stochasticforcing.fields = [{'SMBautoregression'}]; -md.stochasticforcing.dimensions = [md.smb.num_basins]; %dimension of each field md.stochasticforcing.covariance = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields md.stochasticforcing.randomflag = 0; %fixed random seeds Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26614) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26615) @@ -56,5 +56,4 @@ md.stochasticforcing.isstochasticforcing = 1 md.stochasticforcing.fields = ['SMBautoregression'] -md.stochasticforcing.dimensions = [md.smb.num_basins] # dimension of each field md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields md.stochasticforcing.randomflag = 0 # fixed random seeds Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26614) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26615) @@ -39,5 +39,4 @@ md.stochasticforcing.isstochasticforcing = 1; md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'}]; -md.stochasticforcing.dimensions = [md.frontalforcings.num_basins]; %dimension of each field md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields md.stochasticforcing.randomflag = 0; %determines true/false randomness Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26614) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26615) @@ -49,5 +49,4 @@ md.stochasticforcing.isstochasticforcing = 1 md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression'] -md.stochasticforcing.dimensions = [md.frontalforcings.num_basins] # dimension of each field md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields md.stochasticforcing.randomflag = 0 # determines true/false randomness