Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26495) @@ -93,8 +93,11 @@ break; case FrontalForcingsRignotEnum: - iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum); - iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum); - iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum); - break; + /*Retrieve thermal forcing only in the case of non-autoregressive FrontalForcingsRignot*/ + iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum); + /* Do not break here, still retrieve basin_ID,subglacial_discharge, etc.*/ + case FrontalForcingsRignotAutoregressionEnum: + iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum); + break; default: _error_("Frontal forcings"<FindConstant(&melt_parameterization,"md.frontalforcings.parameterization"); + int M,N; + IssmDouble* transparam = NULL; switch(melt_parameterization){ case FrontalForcingsDefaultEnum: break; + 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)); + parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_timestep",FrontalForcingsAutoregressionTimestepEnum)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta0"); + parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta0Enum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta1"); + parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta1Enum,transparam,N)); + xDelete(transparam); + 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*/ case FrontalForcingsRignotEnum: parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum)); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26495) @@ -59,4 +59,100 @@ } 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(); + int basinid; + 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==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + for(int i=0;i(numvertices*arorder); + for(int i=0;i(oldvarspin); + } + 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); +}/*}}}*/ +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; + 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); + 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(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); + xDelete(temparray); + } + + //if(this->lid==55){_printf_("varlist: "<AddInput(SmbMassBalanceEnum,varlist,P1Enum); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum); + + /*Cleanup*/ + xDelete(phi_basin); + xDelete(varlist); + xDelete(valuesautoregression); }/*}}}*/ void Element::ComputeLambdaS(){/*{{{*/ @@ -3658,91 +3754,4 @@ } /*}}}*/ -void Element::SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); - int basinid; - IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; - IssmDouble* phi_basin = xNew(arorder); - IssmDouble* smbspin = xNewZeroInit(numvertices*arorder); - - /*Get Basin ID and Basin coefficients*/ - this->GetInputValue(&basinid,SmbBasinsIdEnum); - for(int i=0;i(numvertices*arorder); - for(int i=0;i(oldsmbspin); - } - this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder); - - /*Cleanup and return*/ - xDelete(smbspin); - xDelete(phi_basin); -}/*}}}*/ -void Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); - int basinid,M,N; - IssmDouble beta0_basin,beta1_basin,noise_basin; - IssmDouble* phi_basin = xNew(arorder); - IssmDouble* smblist = xNew(numvertices); - IssmDouble* smbvaluesautoregression = NULL; - - /*Get Basin ID and Basin coefficients*/ - this->GetInputValue(&basinid,SmbBasinsIdEnum); - for(int ii=0;iiinputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); - - /*If not AR model timestep: take the old SMB values*/ - if(isstepforar==false){ - /*VV do something with the lapse rate here if needed (12Oct2021)*/ - 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); - xDelete(temparray); - } - /*Add input to element*/ - this->AddInput(SmbMassBalanceEnum,smblist,P1Enum); - - /*Cleanup*/ - xDelete(phi_basin); - xDelete(smblist); - xDelete(smbvaluesautoregression); -}/*}}}*/ void Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26494) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26495) @@ -68,4 +68,6 @@ /*bool AnyActive(void);*/ 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,IssmDouble* noisespin,int enum_type); + void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type); void ComputeLambdaS(void); void ComputeNewDamage(); @@ -178,6 +180,4 @@ void SmbSemic(); int Sid(); - void SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin); - void Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms); void SmbGemb(IssmDouble timeinputs, int count); void StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26495) @@ -6,6 +6,7 @@ #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" +#include "../../shared/Random/random.h" -void FrontalForcingsx(FemModel* femmodel){ +void FrontalForcingsx(FemModel* femmodel){/*{{{*/ /*Recover melt_parameterization*/ @@ -17,4 +18,7 @@ case FrontalForcingsDefaultEnum: break; + case FrontalForcingsRignotAutoregressionEnum: + Thermalforcingautoregressionx(femmodel); + /*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/ case FrontalForcingsRignotEnum: femmodel->IcefrontAreax(); @@ -24,3 +28,124 @@ _error_("Frontal forcings "<parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); + femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* covmat = NULL; + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum); + 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); + + /*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*/ + 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); + } + + /*Cleanup*/ + xDelete(beta0); + xDelete(beta1); + xDelete(phi); + xDelete(noisespin); + xDelete(covmat); +}/*}}}*/ +void Thermalforcingautoregressionx(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); + femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum); + + /*Initialize module at first time step*/ + if(time<=starttime+dt){ThermalforcingautoregressionInitx(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,FrontalForcingsNumberofBasinsEnum); + femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); + IssmDouble tinit_ar; + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* covmat = NULL; + IssmDouble* noiseterms = xNew(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); + + /*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*/ + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum); + } + + /*Cleanup*/ + xDelete(beta0); + xDelete(beta1); + xDelete(phi); + xDelete(noiseterms); + xDelete(covmat); +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h (revision 26494) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h (revision 26495) @@ -7,4 +7,6 @@ /* local prototypes: */ void FrontalForcingsx(FemModel* femmodel); +void ThermalforcingautoregressionInitx(FemModel* femmodel); +void Thermalforcingautoregressionx(FemModel* femmodel); #endif Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26495) @@ -154,8 +154,8 @@ femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); - IssmDouble* beta0 = xNew(numbasins); - IssmDouble* beta1 = xNew(numbasins); - IssmDouble* phi = xNew(numbasins*arorder); - IssmDouble* covmat = xNew(numbasins*numbasins); + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* covmat = NULL; femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum); @@ -175,5 +175,6 @@ int fixedseed; for(int i=0;i(numbasins); + //IssmDouble* temparray = xNew(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; @@ -189,5 +190,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); //generate element object - element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin); + element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,SMBautoregressionEnum); } @@ -224,9 +225,9 @@ femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = xNew(numbasins); - IssmDouble* beta1 = xNew(numbasins); - IssmDouble* phi = xNew(numbasins*arorder); + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* covmat = NULL; IssmDouble* noiseterms = xNew(numbasins); - IssmDouble* covmat = xNew(numbasins*numbasins); femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); @@ -247,5 +248,9 @@ if(randomflag) fixedseed=-1; else fixedseed = reCast((time-starttime)/dt); - multivariateNormal(&noiseterms,numbasins,0.0,covmat,fixedseed); + /*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()); @@ -254,5 +259,5 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Smbautoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms); + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum); } Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26494) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26495) @@ -177,6 +177,14 @@ syn keyword cConstant FrictionVoidRatioEnum syn keyword cConstant FrontalForcingsBasinIcefrontAreaEnum +syn keyword cConstant FrontalForcingsAutoregressionInitialTimeEnum +syn keyword cConstant FrontalForcingsAutoregressionTimestepEnum +syn keyword cConstant FrontalForcingsAutoregressiveOrderEnum +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 @@ -1016,4 +1024,5 @@ syn keyword cConstant TemperaturePicardEnum syn keyword cConstant TemperatureSEMICEnum +syn keyword cConstant ThermalforcingValuesAutoregressionEnum syn keyword cConstant ThermalSpctemperatureEnum syn keyword cConstant ThicknessAbsGradientEnum @@ -1262,4 +1271,5 @@ syn keyword cConstant FrontalForcingsDefaultEnum syn keyword cConstant FrontalForcingsRignotEnum +syn keyword cConstant FrontalForcingsRignotAutoregressionEnum syn keyword cConstant FsetEnum syn keyword cConstant FullMeltOnPartiallyFloatingEnum @@ -1553,4 +1563,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1559,6 +1570,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1573,6 +1584,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1581,9 +1592,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1601,4 +1613,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1611,4 +1624,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1623,6 +1637,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1638,11 +1652,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1663,8 +1677,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26494) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26495) @@ -171,6 +171,14 @@ FrictionVoidRatioEnum, FrontalForcingsBasinIcefrontAreaEnum, + FrontalForcingsAutoregressionInitialTimeEnum, + FrontalForcingsAutoregressionTimestepEnum, + FrontalForcingsAutoregressiveOrderEnum, + FrontalForcingsBeta0Enum, + FrontalForcingsBeta1Enum, + FrontalForcingsCovmatEnum, FrontalForcingsNumberofBasinsEnum, FrontalForcingsParamEnum, + FrontalForcingsPhiEnum, + FrontalForcingsRandomflagEnum, GrdModelEnum, GroundinglineFrictionInterpolationEnum, @@ -1013,4 +1021,5 @@ TemperaturePicardEnum, TemperatureSEMICEnum, + ThermalforcingValuesAutoregressionEnum, ThermalSpctemperatureEnum, ThicknessAbsGradientEnum, @@ -1261,4 +1270,5 @@ FrontalForcingsDefaultEnum, FrontalForcingsRignotEnum, + FrontalForcingsRignotAutoregressionEnum, FsetEnum, FullMeltOnPartiallyFloatingEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26495) @@ -179,6 +179,14 @@ case FrictionVoidRatioEnum : return "FrictionVoidRatio"; case FrontalForcingsBasinIcefrontAreaEnum : return "FrontalForcingsBasinIcefrontArea"; + case FrontalForcingsAutoregressionInitialTimeEnum : return "FrontalForcingsAutoregressionInitialTime"; + case FrontalForcingsAutoregressionTimestepEnum : return "FrontalForcingsAutoregressionTimestep"; + case FrontalForcingsAutoregressiveOrderEnum : return "FrontalForcingsAutoregressiveOrder"; + 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"; @@ -1018,4 +1026,5 @@ case TemperaturePicardEnum : return "TemperaturePicard"; case TemperatureSEMICEnum : return "TemperatureSEMIC"; + case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression"; case ThermalSpctemperatureEnum : return "ThermalSpctemperature"; case ThicknessAbsGradientEnum : return "ThicknessAbsGradient"; @@ -1264,4 +1273,5 @@ case FrontalForcingsDefaultEnum : return "FrontalForcingsDefault"; case FrontalForcingsRignotEnum : return "FrontalForcingsRignot"; + case FrontalForcingsRignotAutoregressionEnum : return "FrontalForcingsRignotAutoregression"; case FsetEnum : return "Fset"; case FullMeltOnPartiallyFloatingEnum : return "FullMeltOnPartiallyFloating"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26495) @@ -182,6 +182,14 @@ else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; else if (strcmp(name,"FrontalForcingsBasinIcefrontArea")==0) return FrontalForcingsBasinIcefrontAreaEnum; + else if (strcmp(name,"FrontalForcingsAutoregressionInitialTime")==0) return FrontalForcingsAutoregressionInitialTimeEnum; + else if (strcmp(name,"FrontalForcingsAutoregressionTimestep")==0) return FrontalForcingsAutoregressionTimestepEnum; + else if (strcmp(name,"FrontalForcingsAutoregressiveOrder")==0) return FrontalForcingsAutoregressiveOrderEnum; + 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; @@ -252,5 +260,8 @@ 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 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; @@ -260,8 +271,5 @@ else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum; else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; + else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; else if (strcmp(name,"LoveCoreMantleBoundary")==0) return LoveCoreMantleBoundaryEnum; @@ -375,5 +383,8 @@ 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 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; @@ -383,8 +394,5 @@ else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; + else if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum; else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum; @@ -498,5 +506,8 @@ 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 stage=5; + } + if(stage==5){ + if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum; else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum; else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum; @@ -506,8 +517,5 @@ else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; + else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum; else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum; @@ -621,5 +629,8 @@ else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; - else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"Converged")==0) return ConvergedEnum; @@ -629,8 +640,5 @@ else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum; else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"DamageF")==0) return DamageFEnum; + else if (strcmp(name,"DamageF")==0) return DamageFEnum; else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum; @@ -744,5 +752,8 @@ else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; - else if (strcmp(name,"Ice")==0) return IceEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"Input")==0) return InputEnum; @@ -752,8 +763,5 @@ else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; @@ -867,5 +875,8 @@ else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum; else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum; - else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum; else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum; else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum; @@ -875,8 +886,5 @@ else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum; else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; + else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; @@ -990,5 +998,8 @@ else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbValuesAutoregression")==0) return SmbValuesAutoregressionEnum; - else if (strcmp(name,"SmbV")==0) return SmbVEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SmbV")==0) return SmbVEnum; else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; @@ -998,8 +1009,5 @@ else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum; else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; + else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; else if (strcmp(name,"SmbZY")==0) return SmbZYEnum; else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum; @@ -1042,4 +1050,5 @@ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum; + else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum; else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum; else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; @@ -1112,5 +1121,8 @@ else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; - else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; @@ -1121,8 +1133,5 @@ else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum; else if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; + else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; @@ -1235,5 +1244,8 @@ else if (strcmp(name,"Contact")==0) return ContactEnum; else if (strcmp(name,"Contour")==0) return ContourEnum; - else if (strcmp(name,"Contours")==0) return ContoursEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum; @@ -1244,8 +1256,5 @@ else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; + else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; else if (strcmp(name,"DataSet")==0) return DataSetEnum; @@ -1294,4 +1303,5 @@ else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; + else if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum; @@ -1357,5 +1367,8 @@ else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; - else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; @@ -1367,8 +1380,5 @@ else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"LoveKi")==0) return LoveKiEnum; + else if (strcmp(name,"LoveKi")==0) return LoveKiEnum; else if (strcmp(name,"LoveKr")==0) return LoveKrEnum; else if (strcmp(name,"LoveLi")==0) return LoveLiEnum; @@ -1480,5 +1490,8 @@ else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; - else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; @@ -1490,8 +1503,5 @@ else if (strcmp(name,"SealevelchangePolarMotion")==0) return SealevelchangePolarMotionEnum; else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; + else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; else if (strcmp(name,"SealevelchangeAnalysis")==0) return SealevelchangeAnalysisEnum; else if (strcmp(name,"Seg")==0) return SegEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 26494) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 26495) @@ -274,4 +274,5 @@ case 1: return FrontalForcingsDefaultEnum; case 2: return FrontalForcingsRignotEnum; + case 55: return FrontalForcingsRignotAutoregressionEnum; default: _error_("Marshalled Frontalforcings code \""<',0); + md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'nan',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','frontalforcings.beta0','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); + md = checkfield(md,'fieldname','frontalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); + md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','frontalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); + 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 % }}} + function disp(self) % {{{ + disp(sprintf(' Frontalforcings parameters:')); + fielddisplay(self,'num_basins','number of different basins [unitless]'); + fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); + fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]'); + fielddisplay(self,'beta0','basin-specific intercept values [∘C]'); + fielddisplay(self,'beta1','basin-specific trend values [∘C yr^(-1)]'); + fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); + fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); + 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) % {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',55,'format','Integer'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_order','format','Integer'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_initialtime','format','Double','scale',yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_timestep','format','Double','scale',yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','basin_id','data',self.basin_id-1,'name','md.frontalforcings.basin_id','format','IntMat','mattype',2); %0-indexed + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta0','format','DoubleMat','name','md.frontalforcings.beta0'); + 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 +end Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26494) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26495) @@ -31,4 +31,20 @@ end end + +%SMB parameters +md.timestepping.start_time = 0; +md.timestepping.time_step = 1; +md.timestepping.final_time = 5; +md.smb = SMBautoregression(); +md.smb.num_basins = 3; %number of basins +md.smb.basin_id = idbasin; %prescribe basin ID number to elements +md.smb.beta0 = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr] +md.smb.beta1 = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2] +md.smb.ar_initialtime = md.timestepping.start_time; +md.smb.ar_order = 4; +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 md=solve(md,'Transient'); Index: /issm/trunk-jpl/test/NightlyRun/test542.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26494) +++ /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26495) @@ -9,5 +9,5 @@ %separate domain in 2 basins idbasin = zeros(md.mesh.numberofelements,1); -iid1 = find(md.mesh.y>=mean(md.mesh.y)); +iid1 = find(md.mesh.x<=-1.6e6); for ii=1:md.mesh.numberofelements for vertex=1:3 Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26495) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26495) @@ -0,0 +1,87 @@ +%Test Name: PigTranFrontalforcingsrignotautoregression +md=triangle(model(),'../Exp/Pig.exp',10000.); +md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); +md=parameterize(md,'../Par/Pig.par'); +md=setflowequation(md,'SSA','all'); +md.timestepping.start_time = 0; +md.timestepping.time_step = 1; +md.timestepping.final_time = 25; + +%separate domain in 2 basins +idbasin = zeros(md.mesh.numberofelements,1); +iid1 = find(md.mesh.x<=-1.6e6); +for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 + idbasin(ii) = 1; + end + end + if idbasin(ii)==0 %no vertex was found in basin 1 + idbasin(ii) = 2; + end +end + +%calving and frontalforcings parameters +md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5); +md.calving=calvingvonmises(); +md.frontalforcings=frontalforcingsrignotautoregression(); +md.frontalforcings.num_basins = 2; +md.frontalforcings.basin_id = idbasin; +md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1); +md.frontalforcings.beta0 = [0.05,0.01]; %intercept values of TF in basins [C] +md.frontalforcings.beta1 = [0.001,0.0001]; %trend values of TF in basins [C/yr] +md.frontalforcings.ar_initialtime = md.timestepping.start_time; %initial time in the AR model parameterization [yr] +md.frontalforcings.ar_order = 4; +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; + +md.transient.ismovingfront = 1; +md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); +md.levelset.migration_max = 1e10; + +md.transient.requested_outputs = {'default','CalvingMeltingrate'}; +md.cluster=generic('name',oshostname(),'np',2); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={... + 'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1', 'CalvingMeltingrate1',... + 'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2', 'CalvingMeltingrate2',... + 'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10', 'CalvingMeltingrate10',... + }; +field_tolerances={... + 1e-11,2e-11,2e-11,1e-12,2e-11,6e-12,9e-12,1e-12,1e-9,... + 2e-11,1e-11,1e-11,9e-12,2e-1,2e-11,2e-11,1e-11,1e-9,... + 2e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-9,1e-9,... + }; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).CalvingMeltingrate),... + (md.results.TransientSolution(1).MaskIceLevelset),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).MaskIceLevelset),... + (md.results.TransientSolution(2).CalvingMeltingrate),... + (md.results.TransientSolution(10).Vx),... + (md.results.TransientSolution(10).Vy),... + (md.results.TransientSolution(10).Vel),... + (md.results.TransientSolution(10).Pressure),... + (md.results.TransientSolution(10).Base),... + (md.results.TransientSolution(10).Surface),... + (md.results.TransientSolution(10).Thickness),... + (md.results.TransientSolution(10).MaskIceLevelset),... + (md.results.TransientSolution(10).CalvingMeltingrate),... + };