Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 26525) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 26526) @@ -247,4 +247,5 @@ ./modules/ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.cpp \ ./modules/Solverx/Solverx.cpp \ + ./modules/StochasticForcingx/StochasticForcingx.cpp \ ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \ ./cores/ProcessArguments.cpp \ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26526) @@ -147,5 +147,4 @@ case FrontalForcingsRignotAutoregressionEnum: /*Retrieve autoregressive parameters*/ - parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.randomflag",FrontalForcingsRandomflagEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsAutoregressiveOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_initialtime",FrontalForcingsAutoregressionInitialTimeEnum)); @@ -159,7 +158,4 @@ iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.phi"); parameters->AddObject(new DoubleMatParam(FrontalForcingsPhiEnum,transparam,M,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.covmat"); - parameters->AddObject(new DoubleMatParam(FrontalForcingsCovmatEnum,transparam,M,N)); xDelete(transparam); /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26526) @@ -60,14 +60,14 @@ return false; }/*}}}*/ -void Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin,int enum_type){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); +void Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type){/*{{{*/ + +const int numvertices = this->GetNumberOfVertices(); int basinid; - IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; + IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; IssmDouble* phi_basin = xNew(arorder); IssmDouble* varspin = xNewZeroInit(numvertices*arorder); /*Get Basin ID and Basin coefficients*/ - if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); + if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); for(int i=0;i(numvertices*arorder); for(int i=0;i(oldvarspin); + /*Adjust older values in varspin*/ + for(int i=0;i<(arorder-1)*numvertices;i++) varspin[i+numvertices]=oldvarspin[i]; + + xDelete(oldvarspin); } - if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); + if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); - /*Cleanup and return*/ - xDelete(varspin); - xDelete(phi_basin); + xDelete(varspin); + 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(); + + const int numvertices = this->GetNumberOfVertices(); int basinid,M,N; - IssmDouble beta0_basin,beta1_basin,noise_basin; + IssmDouble beta0_basin,beta1_basin,noise_basin; IssmDouble* phi_basin = xNew(arorder); IssmDouble* varlist = xNew(numvertices); @@ -113,5 +111,5 @@ if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); - for(int ii=0;ii(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); + if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); xDelete(temparray); } - - //if(this->lid==55){_printf_("varlist: "<AddInput(SmbMassBalanceEnum,varlist,P1Enum); if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26526) @@ -127,5 +127,5 @@ /*parameters: */ bool isstressbalance,ismasstransport,isoceantransport,issmb,isthermal,isgroundingline,isesa,issampling;; - bool isslc,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,save_results; + bool isslc,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,isstochasticforcing,save_results; int step,sb_coupling_frequency; int domaintype,numoutputs; @@ -150,8 +150,13 @@ femmodel->parameters->FindParam(&issampling,TransientIssamplingEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); + femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); #if defined(_HAVE_OCEAN_ ) if(isoceancoupling) OceanExchangeDatax(femmodel,false); #endif + + if(isstochasticforcing){ + StochasticForcingx(femmodel); + } if(isthermal && domaintype==Domain3DEnum){ Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26526) @@ -39,5 +39,4 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* covmat = NULL; femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum); @@ -46,30 +45,10 @@ femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum); _assert_(M==numbasins); femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins); - /*AR model spin-up*/ - int nspin{2*arorder+5}; - IssmDouble* noisespin = xNewZeroInit(numbasins*nspin); - my_rank=IssmComm::GetRank(); - if(my_rank==0){ - bool randomflag{}; - femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum); - int fixedseed; - for(int i=0;i(temparray); - } - } - ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Initialize ThermalforcingValuesAutoregressionEnum*/ + /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum*/ + int nspin{2*arorder+5}; for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); //generate element object - element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,FrontalForcingsRignotAutoregressionEnum); + element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum); } @@ -78,6 +57,4 @@ xDelete(beta1); xDelete(phi); - xDelete(noisespin); - xDelete(covmat); }/*}}}*/ void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/ @@ -102,4 +79,5 @@ /*Load parameters*/ + bool isstochastic; int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); @@ -109,33 +87,28 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* covmat = NULL; - IssmDouble* noiseterms = xNew(numbasins); - femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); + IssmDouble* noiseterms = xNewZeroInit(numbasins); + + femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum); _assert_(M==numbasins); femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum); _assert_(M==numbasins); femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins); + /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ + femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); + if(isstochastic){ + 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,ThermalforcingAutoregressionNoiseEnum); _assert_(M==numbasins); + } + } + } /*Time elapsed with respect to AR model initial time*/ IssmDouble telapsed_ar = time-tinit_ar; - /*Before looping through elements: compute noise term specific to each basin from covmat*/ - my_rank=IssmComm::GetRank(); - if(my_rank==0){ - bool randomflag{}; - femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum); - /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/ - int fixedseed; - if(randomflag) fixedseed=-1; - else fixedseed = reCast((time-starttime)/dt); - /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ - IssmDouble* temparray = NULL; - multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed); - for(int i=0;i(temparray); - } - ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Loop over each element to compute SMB at vertices*/ + /*Loop over each element to compute Thermal Forcing at vertices*/ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); @@ -148,4 +121,3 @@ xDelete(phi); xDelete(noiseterms); - xDelete(covmat); }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26526) @@ -90,4 +90,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.issampling",TransientIssamplingEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.isstochasticforcing",StochasticForcingIsStochasticForcingEnum)); /*For stress balance only*/ @@ -399,5 +400,4 @@ /*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)); @@ -412,7 +412,4 @@ parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N)); xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.covmat"); - parameters->AddObject(new DoubleMatParam(SmbCovmatEnum,transparam,M,N)); - xDelete(transparam); break; case SMBgembEnum: @@ -500,4 +497,30 @@ if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs)); iomodel->DeleteData(&requestedoutputs,numoutputs,"md.damage.requested_outputs"); + } + + 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->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); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26526) @@ -157,5 +157,4 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* covmat = NULL; femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum); @@ -164,39 +163,15 @@ 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(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins); - /*AR model spin-up*/ - int nspin{2*arorder+5}; - IssmDouble* noisespin = xNewZeroInit(numbasins*nspin); - my_rank=IssmComm::GetRank(); - if(my_rank==0){ - bool randomflag{}; - femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum); - int fixedseed; - for(int i=0;i(numbasins); - IssmDouble* temparray = NULL; - /*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); - } - } - ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Initialize SmbValuesAutoregressionEnum*/ + /*AR model spin-up with 0 noise to initialize SmbValuesAutoregressionEnum*/ + int nspin{2*arorder+5}; for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); //generate element object - element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,SMBautoregressionEnum); - } - + element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,SMBautoregressionEnum); + } /*Cleanup*/ xDelete(beta0); xDelete(beta1); xDelete(phi); - xDelete(noisespin); - xDelete(covmat); }/*}}}*/ void Smbautoregressionx(FemModel* femmodel){/*{{{*/ @@ -221,4 +196,5 @@ /*Load parameters*/ + bool isstochastic; int M,N,Nphi,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); @@ -228,31 +204,27 @@ IssmDouble* beta1 = NULL; IssmDouble* phi = NULL; - IssmDouble* covmat = NULL; - IssmDouble* noiseterms = xNew(numbasins); + 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(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins); - + + /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/ + IssmDouble* noiseterms = xNewZeroInit(numbasins); + femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); + if(isstochastic){ + 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; - - /*Before looping through elements: compute noise term specific to each basin from covmat*/ - my_rank=IssmComm::GetRank(); - if(my_rank==0){ - bool randomflag{}; - femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum); - /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/ - int fixedseed; - if(randomflag) fixedseed=-1; - else fixedseed = reCast((time-starttime)/dt); - /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/ - IssmDouble* temparray = NULL; - multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed); - for(int i=0;i(temparray); - } - ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); /*Loop over each element to compute SMB at vertices*/ @@ -267,5 +239,4 @@ xDelete(phi); xDelete(noiseterms); - xDelete(covmat); }/*}}}*/ void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 26525) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 26526) @@ -91,4 +91,5 @@ #include "./Scotchx/Scotchx.h" #include "./Shp2Kmlx/Shp2Kmlx.h" +#include "./StochasticForcingx/StochasticForcingx.h" #include "./SurfaceMassBalancex/SurfaceMassBalancex.h" #include "./Solverx/Solverx.h" Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26525) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26526) @@ -182,9 +182,7 @@ syn keyword cConstant FrontalForcingsBeta0Enum syn keyword cConstant FrontalForcingsBeta1Enum -syn keyword cConstant FrontalForcingsCovmatEnum syn keyword cConstant FrontalForcingsNumberofBasinsEnum syn keyword cConstant FrontalForcingsParamEnum syn keyword cConstant FrontalForcingsPhiEnum -syn keyword cConstant FrontalForcingsRandomflagEnum syn keyword cConstant GrdModelEnum syn keyword cConstant GroundinglineFrictionInterpolationEnum @@ -399,4 +397,10 @@ syn keyword cConstant SolidearthSettingsGrdOceanEnum syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum +syn keyword cConstant StochasticForcingCovarianceEnum +syn keyword cConstant StochasticForcingDimensionsEnum +syn keyword cConstant StochasticForcingFieldsEnum +syn keyword cConstant StochasticForcingIsStochasticForcingEnum +syn keyword cConstant StochasticForcingNumFieldsEnum +syn keyword cConstant StochasticForcingRandomflagEnum syn keyword cConstant RotationalPolarMoiEnum syn keyword cConstant SolidearthSettingsReltolEnum @@ -424,4 +428,5 @@ syn keyword cConstant SmbAdThreshEnum syn keyword cConstant SmbAutoregressionInitialTimeEnum +syn keyword cConstant SmbAutoregressionNoiseEnum syn keyword cConstant SmbAutoregressionTimestepEnum syn keyword cConstant SmbAutoregressiveOrderEnum @@ -429,5 +434,4 @@ syn keyword cConstant SmbBeta0Enum syn keyword cConstant SmbBeta1Enum -syn keyword cConstant SmbCovmatEnum syn keyword cConstant SmbDesfacEnum syn keyword cConstant SmbDpermilEnum @@ -462,5 +466,4 @@ syn keyword cConstant SmbPfacEnum syn keyword cConstant SmbPhiEnum -syn keyword cConstant SmbRandomflagEnum syn keyword cConstant SmbRdlEnum syn keyword cConstant SmbRequestedOutputsEnum @@ -498,4 +501,5 @@ syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum syn keyword cConstant StressbalanceShelfDampeningEnum +syn keyword cConstant ThermalforcingAutoregressionNoiseEnum syn keyword cConstant ThermalIsdrainicecolumnEnum syn keyword cConstant ThermalIsdynamicbasalspcEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26525) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26526) @@ -176,9 +176,7 @@ FrontalForcingsBeta0Enum, FrontalForcingsBeta1Enum, - FrontalForcingsCovmatEnum, FrontalForcingsNumberofBasinsEnum, FrontalForcingsParamEnum, FrontalForcingsPhiEnum, - FrontalForcingsRandomflagEnum, GrdModelEnum, GroundinglineFrictionInterpolationEnum, @@ -393,4 +391,10 @@ SolidearthSettingsGrdOceanEnum, SolidearthSettingsOceanAreaScalingEnum, + StochasticForcingCovarianceEnum, + StochasticForcingDimensionsEnum, + StochasticForcingFieldsEnum, + StochasticForcingIsStochasticForcingEnum, + StochasticForcingNumFieldsEnum, + StochasticForcingRandomflagEnum, RotationalPolarMoiEnum, SolidearthSettingsReltolEnum, @@ -418,10 +422,10 @@ SmbAdThreshEnum, SmbAutoregressionInitialTimeEnum, - SmbAutoregressionTimestepEnum, + SmbAutoregressionNoiseEnum, + SmbAutoregressionTimestepEnum, SmbAutoregressiveOrderEnum, SmbAveragingEnum, SmbBeta0Enum, SmbBeta1Enum, - SmbCovmatEnum, SmbDesfacEnum, SmbDpermilEnum, @@ -456,5 +460,4 @@ SmbPfacEnum, SmbPhiEnum, - SmbRandomflagEnum, SmbRdlEnum, SmbRequestedOutputsEnum, @@ -492,4 +495,5 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, + ThermalforcingAutoregressionNoiseEnum, ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26526) @@ -184,9 +184,7 @@ case FrontalForcingsBeta0Enum : return "FrontalForcingsBeta0"; case FrontalForcingsBeta1Enum : return "FrontalForcingsBeta1"; - case FrontalForcingsCovmatEnum : return "FrontalForcingsCovmat"; case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; case FrontalForcingsParamEnum : return "FrontalForcingsParam"; case FrontalForcingsPhiEnum : return "FrontalForcingsPhi"; - case FrontalForcingsRandomflagEnum : return "FrontalForcingsRandomflag"; case GrdModelEnum : return "GrdModel"; case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation"; @@ -401,4 +399,10 @@ case SolidearthSettingsGrdOceanEnum : return "SolidearthSettingsGrdOcean"; case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling"; + case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance"; + case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions"; + case StochasticForcingFieldsEnum : return "StochasticForcingFields"; + case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing"; + case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields"; + case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag"; case RotationalPolarMoiEnum : return "RotationalPolarMoi"; case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol"; @@ -426,4 +430,5 @@ case SmbAdThreshEnum : return "SmbAdThresh"; case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime"; + case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise"; case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep"; case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder"; @@ -431,5 +436,4 @@ case SmbBeta0Enum : return "SmbBeta0"; case SmbBeta1Enum : return "SmbBeta1"; - case SmbCovmatEnum : return "SmbCovmat"; case SmbDesfacEnum : return "SmbDesfac"; case SmbDpermilEnum : return "SmbDpermil"; @@ -464,5 +468,4 @@ case SmbPfacEnum : return "SmbPfac"; case SmbPhiEnum : return "SmbPhi"; - case SmbRandomflagEnum : return "SmbRandomflag"; case SmbRdlEnum : return "SmbRdl"; case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; @@ -500,4 +503,5 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; + case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise"; case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26525) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26526) @@ -187,9 +187,7 @@ else if (strcmp(name,"FrontalForcingsBeta0")==0) return FrontalForcingsBeta0Enum; else if (strcmp(name,"FrontalForcingsBeta1")==0) return FrontalForcingsBeta1Enum; - else if (strcmp(name,"FrontalForcingsCovmat")==0) return FrontalForcingsCovmatEnum; else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; else if (strcmp(name,"FrontalForcingsPhi")==0) return FrontalForcingsPhiEnum; - else if (strcmp(name,"FrontalForcingsRandomflag")==0) return FrontalForcingsRandomflagEnum; else if (strcmp(name,"GrdModel")==0) return GrdModelEnum; else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; @@ -260,10 +258,10 @@ else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum; else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum; + else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum; + else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum; - else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum; - else if (strcmp(name,"InversionType")==0) return InversionTypeEnum; + if (strcmp(name,"InversionType")==0) return InversionTypeEnum; else if (strcmp(name,"Ivins")==0) return IvinsEnum; else if (strcmp(name,"IsSlcCoupling")==0) return IsSlcCouplingEnum; @@ -383,10 +381,10 @@ else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum; else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum; + else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum; + else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum; - else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum; - else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum; + if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum; else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum; else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; @@ -410,4 +408,10 @@ else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum; else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum; + else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum; + else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum; + else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum; + else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum; + else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum; + else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum; else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum; else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; @@ -435,4 +439,5 @@ 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; @@ -440,5 +445,4 @@ else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum; else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum; - else if (strcmp(name,"SmbCovmat")==0) return SmbCovmatEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; @@ -473,5 +477,4 @@ 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; @@ -503,13 +506,14 @@ else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum; else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; - else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum; - else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum; - else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum; + if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum; + else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum; + else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum; + else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum; 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; @@ -625,12 +629,12 @@ else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; - else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; + else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"Converged")==0) return ConvergedEnum; @@ -748,12 +752,12 @@ else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum; else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum; - else if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum; else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"Ice")==0) return IceEnum; + else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"Input")==0) return InputEnum; @@ -871,12 +875,12 @@ else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum; else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum; - else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum; else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum; else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum; else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum; + else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum; else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum; else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum; @@ -994,12 +998,12 @@ else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; - else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbValuesAutoregression")==0) return SmbValuesAutoregressionEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SmbV")==0) return SmbVEnum; + else if (strcmp(name,"SmbV")==0) return SmbVEnum; else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; @@ -1117,12 +1121,12 @@ else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum; - else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; + else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; @@ -1240,12 +1244,12 @@ else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; else if (strcmp(name,"Closed")==0) return ClosedEnum; - else if (strcmp(name,"Colinear")==0) return ColinearEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Colinear")==0) return ColinearEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Contact")==0) return ContactEnum; else if (strcmp(name,"Contour")==0) return ContourEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Contours")==0) return ContoursEnum; + else if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum; @@ -1363,12 +1367,12 @@ else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum; else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; - else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; @@ -1486,12 +1490,12 @@ else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; - else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; + else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26525) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26526) @@ -13,7 +13,5 @@ ar_timestep = 0; phi = NaN; - covmat = NaN; basin_id = NaN; - randomflag = 1; steps_per_step = 1; averaging = 0; @@ -57,12 +55,7 @@ disp(' smb.phi (lag coefficients) not specified: order of autoregressive model set to 0'); end - if isnan(self.covmat) - self.covmat = 1e-21*eye(self.num_basins); %no stochasticity and no covariance - disp(' smb.covmat not specified: stochasticity set to 0'); - end end % }}} function self = setdefaultparameters(self) % {{{ self.ar_order = 0.0; %autoregression model of order 0 - self.randomflag = 1; end % }}} function md = checkconsistency(self,md,solution,analyses) @@ -77,6 +70,4 @@ md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); - md = checkfield(md,'fieldname','smb.covmat','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.num_basins]); - md = checkfield(md,'fieldname','smb.randomflag','numel',[1],'values',[0 1]); end md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); @@ -94,6 +85,4 @@ fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); 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'); @@ -117,6 +106,4 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',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'); Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 26525) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 26526) @@ -13,7 +13,5 @@ ar_timestep = 0; phi = NaN; - covmat = NaN; basin_id = NaN; - randomflag = 1; subglacial_discharge = NaN; end @@ -46,5 +44,4 @@ subglacial_discharge = NaN; self.ar_order = 0.0; %autoregression model of order 0 - self.randomflag = 1; end % }}} @@ -62,6 +59,4 @@ md = checkfield(md,'fieldname','frontalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','frontalforcings.phi','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); - md = checkfield(md,'fieldname','frontalforcings.covmat','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.num_basins]); - md = checkfield(md,'fieldname','frontalforcings.randomflag','numel',[1],'values',[0 1]); end % }}} @@ -77,6 +72,4 @@ fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); 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 [∘C^2]'); - fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -92,6 +85,4 @@ WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta1','format','DoubleMat','name','md.frontalforcings.beta1','scale',1./yts,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','covmat','format','DoubleMat','name','md.frontalforcings.covmat'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','randomflag','format','Boolean'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 26525) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 26526) @@ -55,4 +55,5 @@ miscellaneous = 0; private = 0; + stochasticforcing= 0; %}}} @@ -191,4 +192,6 @@ %2021 February 17 if isa(md.sampling,'double'); md.sampling=sampling(); end + %VV + if ~isa(md.stochasticforcing,'stochasticforcing'); md.stochasticforcing=stochasticforcing(); end end% }}} end @@ -250,4 +253,5 @@ disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay')); disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); + disp(sprintf('%19s: %-22s -- %s','stochasticforcing',['[1x1 ' class(self.stochasticforcing) ']'],'stochasticity applied to model forcings')); end % }}} function md = setdefaultparameters(md,planet) % {{{ @@ -297,4 +301,5 @@ md.miscellaneous = miscellaneous(); md.private = private(); + md.stochasticforcing= stochasticforcing(); end %}}} Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26526) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26526) @@ -0,0 +1,99 @@ +%STOCHASTICFORCING class definition +% +% Usage: +% stochasticforcing=stochasticforcing(); + +classdef stochasticforcing + properties (SetAccess=public) + isstochasticforcing = 0; + fields = NaN; + dimensions = NaN; + covariance = NaN; + randomflag = 1; + end + methods + function self = stochasticforcing(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = extrude(self,md) % {{{ + %Nothing for now + end % }}} + function self = setdefaultparameters(self) % {{{ + self.isstochasticforcing = 0; %stochasticforcing is turned off by default + self.randomflag = 1; %true randomness is implemented by default + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %Early return + if ~self.isstochasticforcing, return; end + + num_fields = numel(self.fields); + size_tot = 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',[1,num_fields]); %specific dimension for each field + 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 + end % }}} + function disp(self) % {{{ + disp(sprintf(' stochasticforcing parameters:')); + 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,'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)'); + disp('Available fields:'); + disp(' SMBautoregression'); + disp(' FrontalForcingsRignotAutoregression (thermal forcing)'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + num_fields = numel(self.fields); + %Scaling covariance matrix (scale column-by-column and row-by-row) + scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts + 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))]; + for row=inds %scale rows corresponding to scaled field + self.covariance(row,:) = 1./yts*self.covariance(row,:); + end + for col=inds %scale columns corresponding to scaled field + self.covariance(:,col) = 1./yts*self.covariance(:,col); + end + end + end + + WriteData(fid,prefix,'object',self,'fieldname','isstochasticforcing','format','Boolean'); + if ~self.isstochasticforcing + return + else + 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,'object',self,'fieldname','covariance','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean'); + end + end % }}} + end +end Index: /issm/trunk-jpl/src/m/miscellaneous/supportedstochforcings.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/supportedstochforcings.m (revision 26526) +++ /issm/trunk-jpl/src/m/miscellaneous/supportedstochforcings.m (revision 26526) @@ -0,0 +1,8 @@ +function list = supportedstochforcings(), + % Defines list of fields supported + % by the class md.stochasticforcings + + list = {... + 'SMBautoregression',... + 'FrontalForcingsRignotAutoregression' + }; Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26525) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26526) @@ -45,6 +45,11 @@ md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] md.smb.phi = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]]; -md.smb.covmat = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; -md.smb.randomflag = 0; %fixed random seeds + +%Stochastic forcing +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 md=solve(md,'Transient'); Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26525) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26526) @@ -35,6 +35,11 @@ md.frontalforcings.ar_timestep = 2; %timestep of the autoregressive model [yr] md.frontalforcings.phi = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters -md.frontalforcings.covmat = 1e-4*[[1.5,0.5];[0.5,0.4]]; -md.frontalforcings.randomflag = 0; + +%Stochastic forcing +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 md.transient.ismovingfront = 1;