Changeset 26615

Timestamp:

11/12/21 06:25:41 (3 years ago)

Author:

vverjans

Message:

CHG: making stochasticforcing.m more generic, extending stochasticforcing capability to DefaultCalving (second attempt)

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.cpp (modified) (2 diffs)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (6 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• src/c/modules/StochasticForcingx/StochasticForcingx.cpp (modified) (2 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (1 diff)
• src/c/shared/Enum/Enum.vim (modified) (7 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (7 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (7 diffs)
• src/c/shared/Enum/Enumjl.vim (modified) (13 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (12 diffs)
• src/m/classes/stochasticforcing.m (modified) (6 diffs)
• src/m/classes/stochasticforcing.py (modified) (6 diffs)
• test/NightlyRun/test257.m (modified) (1 diff)
• test/NightlyRun/test257.py (modified) (1 diff)
• test/NightlyRun/test543.m (modified) (1 diff)
• test/NightlyRun/test543.py (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp ¶

@@ -83 +83 @@
 
         /*Get moving front parameters*/
-        int  calvinglaw;
-        iomodel->FindConstant(&calvinglaw,"md.calving.law");
-        switch(calvinglaw){
-                case DefaultCalvingEnum:
-                        iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
-                        break;
+        bool isstochastic;
+   int  calvinglaw;
+   iomodel->FindConstant(&calvinglaw,"md.calving.law");
+   iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
+   switch(calvinglaw){
+      case DefaultCalvingEnum:
+         iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
+         if(isstochastic){
+            iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
+            iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum);
+         }
+         break; 
                 case CalvingLevermannEnum:
                         iomodel->FetchDataToInput(inputs,elements,"md.calving.coeff",CalvinglevermannCoeffEnum);

issm/trunk-jpl/src/c/classes/Elements/Element.cpp ¶

@@ -99 +99 @@
    xDelete<IssmDouble>(phi_basin);      
 }/*}}}*/
-void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
-        
-        const int numvertices = this->GetNumberOfVertices();
-   int         basinid,M,N;
-   IssmDouble  beta0_basin,beta1_basin,noise_basin;
+void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/
+
+   const int numvertices = this->GetNumberOfVertices();
+   int         basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
+   IssmDouble  beta0_basin,beta1_basin,noiseterm;
    IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
    IssmDouble* valuesautoregression = NULL;
-
-   /*Get Basin ID and Basin coefficients*/
-   if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
-   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
+   Input*      noiseterm_input      = NULL;
+
+   /*Get field-specific enums*/
+   switch(enum_type){
+      case(SMBautoregressionEnum):
+         arenum_type    = SmbValuesAutoregressionEnum;
+         basinenum_type = SmbBasinsIdEnum;
+         noiseenum_type = SmbAutoregressionNoiseEnum;
+         outenum_type   = SmbMassBalanceEnum;
+         break;
+      case(FrontalForcingsRignotAutoregressionEnum):
+         arenum_type    = ThermalforcingValuesAutoregressionEnum;
+         basinenum_type = FrontalForcingsBasinIdEnum;
+         noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
+         outenum_type   = FrontalForcingsThermalForcingEnum;
+         break;
+   }
+
+   /*Get noise and autoregressive terms*/
+   this->GetInputValue(&basinid,basinenum_type);
+   if(isfieldstochastic){
+      noiseterm_input = this->GetInput(noiseenum_type);
+      Gauss* gauss = this->NewGauss();
+      noiseterm_input->GetInputValue(&noiseterm,gauss);
+      delete gauss;
+   }
+   else noiseterm = 0.0;
+   this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
+
+   /*Get basin coefficients*/
    for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
    beta0_basin   = beta0[basinid];
    beta1_basin   = beta1[basinid];
-   noise_basin   = noiseterms[basinid];
-   if(enum_type==SMBautoregressionEnum)                   this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
-   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
-
-   /*If not AR model timestep: take the old values of variable*/
+
+        /*If not AR model timestep: take the old values of variable*/
    if(isstepforar==false){
       for(int i=0;i<numvertices;i++) varlist[i]=valuesautoregression[i];
@@ -131 +154 @@
 
          /*Stochastic variable value*/
-         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
+         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
       }
-      /*Update autoregression TF values*/
+      /*Update autoregression values*/
       IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
       /*Assign newest values and shift older values*/
       for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
       for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
-      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);
+      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
       xDelete<IssmDouble>(temparray);
    }
 
    /*Add input to element*/
-   if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
-   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
+   this->AddInput(outenum_type,varlist,P1Enum);
 
    /*Cleanup*/

issm/trunk-jpl/src/c/classes/Elements/Element.h ¶

@@ -69 +69 @@
                 bool               AnyFSet(void);
                 void                                     AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type);
-      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
+      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type);
                 void               ComputeLambdaS(void);
                 void               ComputeNewDamage();

issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp ¶

@@ -80 +80 @@
    /*Load parameters*/
         bool isstochastic;
-   int M,N,Nphi,arorder,numbasins,my_rank;
+   bool istfstochastic = false;
+        int M,N,Nphi,arorder,numbasins,my_rank;
    femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
@@ -87 +88 @@
    IssmDouble* beta1      = NULL;
    IssmDouble* phi        = NULL;
-   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
 
         femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
@@ -94 +94 @@
    femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
 
-        /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
         femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
         if(isstochastic){
@@ -102 +101 @@
                 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
                 for(int i=0;i<numstochasticfields;i++){
-                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
-                                femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum);  _assert_(M==numbasins);
-                        }
+                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
                 }
+                xDelete<int>(stochasticfields);
         }
    /*Time elapsed with respect to AR model initial time*/
@@ -113 +111 @@
    for(Object* &object:femmodel->elements->objects){
       Element* element = xDynamicCast<Element*>(object);
-      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
    }
 
@@ -120 +118 @@
    xDelete<IssmDouble>(beta1);
    xDelete<IssmDouble>(phi);
-   xDelete<IssmDouble>(noiseterms);
 }/*}}}*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -500 +500 @@
 
         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<int>(num_fields);
-                for(int i=0;i<num_fields;i++){
-                        stochasticforcing_enums[i] = StringToEnumx(fields[i]);
-                        xDelete<char>(fields[i]);
-                }
-                xDelete<char*>(fields);
-                parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
-                xDelete<int>(stochasticforcing_enums);
-
+   parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
+   if(isstochasticforcing){
+      int num_fields,stochastic_dim;
+      char** fields;
+      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
+      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
+      iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
+      if(num_fields<1) _error_("no stochasticforcing fields found");
+      int* stochasticforcing_enums = xNew<int>(num_fields);
+      for(int i=0;i<num_fields;i++){
+         stochasticforcing_enums[i] = StringToEnumx(fields[i]);
+         xDelete<char>(fields[i]);
+      }
+      xDelete<char*>(fields);
+      parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
+      xDelete<int>(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<IssmDouble>(transparam);
-                iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
-                parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
-                xDelete<IssmDouble>(transparam);
-        }
+      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
+      parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
+      xDelete<IssmDouble>(transparam);
+      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
+      parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
+      xDelete<IssmDouble>(transparam);
+   }
 
         /*Deal with mass flux segments: {{{*/

issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp ¶

@@ -36 +36 @@
       if(randomflag) fixedseed=-1;
       else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
-      /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
+                /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
       IssmDouble* temparray = NULL;
       multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
@@ -46 +46 @@
         int i=0;
    for(int j=0;j<numfields;j++){
-      int enum_type;
+      int dimenum_type,noiseenum_type;
       IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
       for(int k=0;k<dimensions[j];k++){
          noisefield[k]=noiseterms[i+k];
       }
-      
-      if(fields[j]==SMBautoregressionEnum)                        enum_type = SmbAutoregressionNoiseEnum;
-                else if(fields[j]==FrontalForcingsRignotAutoregressionEnum) enum_type = ThermalforcingAutoregressionNoiseEnum;
-                else _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.\n"); 
-      femmodel->parameters->SetParam(noisefield,dimensions[j],enum_type);
-      i=i+dimensions[j];
+     
+                int dimensionid;
+
+                /*Deal with the autoregressive models*/
+                if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
+                        switch(fields[j]){
+                                case SMBautoregressionEnum:
+                                        dimenum_type   = SmbBasinsIdEnum;
+                                        noiseenum_type = SmbAutoregressionNoiseEnum;
+                                        break;
+                                case FrontalForcingsRignotAutoregressionEnum:
+                                        dimenum_type   = FrontalForcingsBasinIdEnum;
+                                        noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
+                                        break;
+                        }
+                        for(Object* &object:femmodel->elements->objects){
+            Element* element = xDynamicCast<Element*>(object);
+            int numvertices  = element->GetNumberOfVertices();
+            IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
+            element->GetInputValue(&dimensionid,dimenum_type);
+            for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
+            element->AddInput(noiseenum_type,noise_element,P0Enum);
+            xDelete<IssmDouble>(noise_element);
+                        }
+                }
+                else{
+                        switch(fields[j]){
+                                case SMBautoregressionEnum:
+                                case FrontalForcingsRignotAutoregressionEnum:
+                                        /*Already done above*/
+                                        break;
+                                case DefaultCalvingEnum:
+                                        /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
+                                        femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
+                                        for(Object* &object:femmodel->elements->objects){
+                                                Element* element = xDynamicCast<Element*>(object);
+                                                int numvertices  = element->GetNumberOfVertices();
+                                                IssmDouble baselinecalvingrate;
+                                                IssmDouble calvingrate_tot[numvertices];
+                                                Input* baselinecalvingrate_input  = NULL;
+                                                baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
+                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
+                                                Gauss* gauss = element->NewGauss();
+                                                for(int i=0;i<numvertices;i++){
+                                                        gauss->GaussVertex(i);
+                                                        baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
+                                                        calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
+                                                }
+                                                element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
+                                                delete gauss;
+                                        }
+                                        break;
+                                default:
+                                        _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
+                        }
+                }
+                i=i+dimensions[j];
       xDelete<IssmDouble>(noisefield);
    }

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp ¶

@@ -177 +177 @@
 void Smbautoregressionx(FemModel* femmodel){/*{{{*/
 
-        /*Get time parameters*/
-        IssmDouble time,dt,starttime,tstep_ar;
-        femmodel->parameters->FindParam(&time,TimeEnum);
-        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   /*Get time parameters*/
+   IssmDouble time,dt,starttime,tstep_ar;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
 
-        /*Initialize module at first time step*/
-        if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
-        /*Determine if this is a time step for the AR model*/
-        bool isstepforar = false;
-
-        #ifndef _HAVE_AD_
+   /*Initialize module at first time step*/
+   if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
+   /*Determine if this is a time step for the AR model*/
+   bool isstepforar = false;
+
+   #ifndef _HAVE_AD_
    if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
-        #else
-        _error_("not implemented yet");
-        #endif
-
-        /*Load parameters*/
-        bool isstochastic;
-        int M,N,Nphi,arorder,numbasins,my_rank;
-        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
-        femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
-        IssmDouble tinit_ar;
-        IssmDouble* beta0      = NULL; 
-        IssmDouble* beta1      = NULL;
-        IssmDouble* phi        = NULL;
-
-        femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
-        femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
-        femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
-        femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
-
-        /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
-   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
-        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
+   #else
+   _error_("not implemented yet");
+   #endif
+
+   /*Load parameters*/
+   bool isstochastic;
+   bool issmbstochastic = false;
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
+   IssmDouble tinit_ar;
+   IssmDouble* beta0      = NULL; 
+   IssmDouble* beta1      = NULL;
+   IssmDouble* phi        = NULL;
+
+   femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+
+   femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    if(isstochastic){
-                int  numstochasticfields;
+      int  numstochasticfields;
       int* stochasticfields;
       femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
       femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
       for(int i=0;i<numstochasticfields;i++){
-         if(stochasticfields[i]==SMBautoregressionEnum){
-                                femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum);  _assert_(M==numbasins);
-                        }
-                }
-                xDelete<int>(stochasticfields);
-        }
-        /*Time elapsed with respect to AR model initial time*/
-        IssmDouble telapsed_ar = time-tinit_ar; 
-
-        /*Loop over each element to compute SMB at vertices*/
-        for(Object* &object:femmodel->elements->objects){
-                Element* element = xDynamicCast<Element*>(object);
-                element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
-        }
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(beta0);
-        xDelete<IssmDouble>(beta1);
-        xDelete<IssmDouble>(phi);
-        xDelete<IssmDouble>(noiseterms);
+         if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true;
+      }
+      xDelete<int>(stochasticfields);
+   }
+   /*Time elapsed with respect to AR model initial time*/
+   IssmDouble telapsed_ar = time-tinit_ar;
+
+   /*Loop over each element to compute SMB at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
+   }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
 }/*}}}*/
 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/

issm/trunk-jpl/src/c/shared/Enum/Enum.vim ¶

@@ -398 +398 @@
 syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
 syn keyword cConstant StochasticForcingCovarianceEnum
+syn keyword cConstant StochasticForcingDefaultDimensionEnum
 syn keyword cConstant StochasticForcingDimensionsEnum
 syn keyword cConstant StochasticForcingFieldsEnum
@@ -428 +429 @@
 syn keyword cConstant SmbAdThreshEnum
 syn keyword cConstant SmbAutoregressionInitialTimeEnum
-syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbAutoregressionTimestepEnum
 syn keyword cConstant SmbAutoregressiveOrderEnum
@@ -504 +504 @@
 syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
 syn keyword cConstant StressbalanceShelfDampeningEnum
-syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalIsdrainicecolumnEnum
 syn keyword cConstant ThermalIsdynamicbasalspcEnum
@@ -602 +601 @@
 syn keyword cConstant BaseSlopeXEnum
 syn keyword cConstant BaseSlopeYEnum
+syn keyword cConstant BaselineCalvingCalvingrateEnum
 syn keyword cConstant BedEnum
 syn keyword cConstant BedGRDEnum
@@ -895 +895 @@
 syn keyword cConstant SmbAdiffiniEnum
 syn keyword cConstant SmbAiniEnum
+syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbBasinsIdEnum
 syn keyword cConstant SmbBMaxEnum
@@ -997 +998 @@
 syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
 syn keyword cConstant SolidearthExternalGeoidRateEnum
+syn keyword cConstant StochasticForcingDefaultIdEnum
 syn keyword cConstant StrainRateeffectiveEnum
 syn keyword cConstant StrainRateparallelEnum
@@ -1032 +1034 @@
 syn keyword cConstant TemperaturePicardEnum
 syn keyword cConstant TemperatureSEMICEnum
+syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalforcingValuesAutoregressionEnum
 syn keyword cConstant ThermalSpctemperatureEnum

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -392 +392 @@
         SolidearthSettingsOceanAreaScalingEnum,
         StochasticForcingCovarianceEnum,
+        StochasticForcingDefaultDimensionEnum,
         StochasticForcingDimensionsEnum,
         StochasticForcingFieldsEnum,
@@ -422 +423 @@
         SmbAdThreshEnum,
         SmbAutoregressionInitialTimeEnum,
-   SmbAutoregressionNoiseEnum,
         SmbAutoregressionTimestepEnum,
    SmbAutoregressiveOrderEnum,
@@ -498 +498 @@
         StressbalanceRiftPenaltyThresholdEnum,
         StressbalanceShelfDampeningEnum,
-   ThermalforcingAutoregressionNoiseEnum,
         ThermalIsdrainicecolumnEnum,
         ThermalIsdynamicbasalspcEnum,
@@ -598 +597 @@
         BaseSlopeXEnum,
         BaseSlopeYEnum,
+        BaselineCalvingCalvingrateEnum,
         BedEnum,
         BedGRDEnum,
@@ -891 +891 @@
         SmbAdiffiniEnum,
         SmbAiniEnum,
+   SmbAutoregressionNoiseEnum,
         SmbBasinsIdEnum,
         SmbBMaxEnum,
@@ -994 +995 @@
         SolidearthExternalDisplacementUpRateEnum,
         SolidearthExternalGeoidRateEnum,
+        StochasticForcingDefaultIdEnum,
         StrainRateeffectiveEnum,
         StrainRateparallelEnum,
@@ -1029 +1031 @@
         TemperaturePicardEnum,
         TemperatureSEMICEnum,
+   ThermalforcingAutoregressionNoiseEnum,
         ThermalforcingValuesAutoregressionEnum, 
         ThermalSpctemperatureEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -400 +400 @@
                 case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
                 case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
+                case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension";
                 case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
                 case StochasticForcingFieldsEnum : return "StochasticForcingFields";
@@ -430 +431 @@
                 case SmbAdThreshEnum : return "SmbAdThresh";
                 case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
-                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
                 case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
                 case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
@@ -506 +506 @@
                 case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
                 case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
-                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
                 case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
                 case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
@@ -604 +603 @@
                 case BaseSlopeXEnum : return "BaseSlopeX";
                 case BaseSlopeYEnum : return "BaseSlopeY";
+                case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
                 case BedEnum : return "Bed";
                 case BedGRDEnum : return "BedGRD";
@@ -897 +897 @@
                 case SmbAdiffiniEnum : return "SmbAdiffini";
                 case SmbAiniEnum : return "SmbAini";
+                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
                 case SmbBasinsIdEnum : return "SmbBasinsId";
                 case SmbBMaxEnum : return "SmbBMax";
@@ -999 +1000 @@
                 case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
                 case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
+                case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
                 case StrainRateeffectiveEnum : return "StrainRateeffective";
                 case StrainRateparallelEnum : return "StrainRateparallel";
@@ -1034 +1036 @@
                 case TemperaturePicardEnum : return "TemperaturePicard";
                 case TemperatureSEMICEnum : return "TemperatureSEMIC";
+                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
                 case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
                 case ThermalSpctemperatureEnum : return "ThermalSpctemperature";

issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim ¶

@@ -170 +170 @@
 syn keyword juliaConstC FrictionVoidRatioEnum
 syn keyword juliaConstC FrontalForcingsBasinIcefrontAreaEnum
+syn keyword juliaConstC FrontalForcingsAutoregressionInitialTimeEnum
+syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum
+syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum
+syn keyword juliaConstC FrontalForcingsBeta0Enum
+syn keyword juliaConstC FrontalForcingsBeta1Enum
 syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
 syn keyword juliaConstC FrontalForcingsParamEnum
+syn keyword juliaConstC FrontalForcingsPhiEnum
 syn keyword juliaConstC GrdModelEnum
 syn keyword juliaConstC GroundinglineFrictionInterpolationEnum
@@ -384 +390 @@
 syn keyword juliaConstC SolidearthSettingsGrdOceanEnum
 syn keyword juliaConstC SolidearthSettingsOceanAreaScalingEnum
+syn keyword juliaConstC StochasticForcingCovarianceEnum
+syn keyword juliaConstC StochasticForcingDefaultDimensionEnum
+syn keyword juliaConstC StochasticForcingDimensionsEnum
+syn keyword juliaConstC StochasticForcingFieldsEnum
+syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
+syn keyword juliaConstC StochasticForcingNumFieldsEnum
+syn keyword juliaConstC StochasticForcingRandomflagEnum
 syn keyword juliaConstC RotationalPolarMoiEnum
 syn keyword juliaConstC SolidearthSettingsReltolEnum
@@ -414 +427 @@
 syn keyword juliaConstC SmbBeta0Enum
 syn keyword juliaConstC SmbBeta1Enum
-syn keyword juliaConstC SmbCovmatEnum
 syn keyword juliaConstC SmbDesfacEnum
 syn keyword juliaConstC SmbDpermilEnum
@@ -424 +436 @@
 syn keyword juliaConstC SmbDtEnum
 syn keyword juliaConstC SmbEnum
+syn keyword juliaConstC SmbEIdxEnum
 syn keyword juliaConstC SmbFEnum
 syn keyword juliaConstC SmbInitDensityScalingEnum
@@ -432 +445 @@
 syn keyword juliaConstC SmbIsdelta18oEnum
 syn keyword juliaConstC SmbIsdensificationEnum
+syn keyword juliaConstC SmbIsdeltaLWupEnum
 syn keyword juliaConstC SmbIsfirnwarmingEnum
 syn keyword juliaConstC SmbIsgraingrowthEnum
@@ -447 +461 @@
 syn keyword juliaConstC SmbPfacEnum
 syn keyword juliaConstC SmbPhiEnum
-syn keyword juliaConstC SmbRandomflagEnum
 syn keyword juliaConstC SmbRdlEnum
 syn keyword juliaConstC SmbRequestedOutputsEnum
@@ -460 +473 @@
 syn keyword juliaConstC SmbT0dryEnum
 syn keyword juliaConstC SmbT0wetEnum
+syn keyword juliaConstC SmbTeThreshEnum
 syn keyword juliaConstC SmbTdiffEnum
 syn keyword juliaConstC SmbThermoDeltaTScalingEnum
@@ -580 +594 @@
 syn keyword juliaConstC BaseSlopeXEnum
 syn keyword juliaConstC BaseSlopeYEnum
+syn keyword juliaConstC BaselineCalvingCalvingrateEnum
 syn keyword juliaConstC BedEnum
 syn keyword juliaConstC BedGRDEnum
@@ -873 +888 @@
 syn keyword juliaConstC SmbAdiffiniEnum
 syn keyword juliaConstC SmbAiniEnum
+syn keyword juliaConstC SmbAutoregressionNoiseEnum
 syn keyword juliaConstC SmbBasinsIdEnum
 syn keyword juliaConstC SmbBMaxEnum
@@ -894 +910 @@
 syn keyword juliaConstC SmbDiniEnum
 syn keyword juliaConstC SmbDlwrfEnum
+syn keyword juliaConstC SmbDulwrfValueEnum
 syn keyword juliaConstC SmbDswrfEnum
 syn keyword juliaConstC SmbDswdiffrfEnum
@@ -974 +991 @@
 syn keyword juliaConstC SolidearthExternalDisplacementUpRateEnum
 syn keyword juliaConstC SolidearthExternalGeoidRateEnum
+syn keyword juliaConstC StochasticForcingDefaultIdEnum
 syn keyword juliaConstC StrainRateeffectiveEnum
 syn keyword juliaConstC StrainRateparallelEnum
@@ -1009 +1027 @@
 syn keyword juliaConstC TemperaturePicardEnum
 syn keyword juliaConstC TemperatureSEMICEnum
+syn keyword juliaConstC ThermalforcingAutoregressionNoiseEnum
+syn keyword juliaConstC ThermalforcingValuesAutoregressionEnum
 syn keyword juliaConstC ThermalSpctemperatureEnum
 syn keyword juliaConstC ThicknessAbsGradientEnum
@@ -1255 +1275 @@
 syn keyword juliaConstC FrontalForcingsDefaultEnum
 syn keyword juliaConstC FrontalForcingsRignotEnum
+syn keyword juliaConstC FrontalForcingsRignotAutoregressionEnum
 syn keyword juliaConstC FsetEnum
 syn keyword juliaConstC FullMeltOnPartiallyFloatingEnum

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -409 +409 @@
               else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
               else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
+              else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
               else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
               else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
@@ -439 +440 @@
               else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
               else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
-              else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
               else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
               else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
@@ -518 +518 @@
               else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
               else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
-              else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
               else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
               else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
@@ -616 +615 @@
               else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
               else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum;
+              else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
               else if (strcmp(name,"Bed")==0) return BedEnum;
               else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
@@ -918 +918 @@
               else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
               else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+              else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
               else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
               else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
@@ -997 +998 @@
               else if (strcmp(name,"SmbT")==0) return SmbTEnum;
               else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
-              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
+              if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
+              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
               else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
               else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
@@ -1023 +1024 @@
               else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
               else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
+              else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum;
               else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
               else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
@@ -1058 +1060 @@
               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
               else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum;
+              else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
               else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum;
               else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
@@ -1118 +1121 @@
               else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
               else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
-              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
-              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
-              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
+              if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
+              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
+              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
+              else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
               else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
               else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
@@ -1241 +1244 @@
               else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
               else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
-              else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
-              else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
-              else if (strcmp(name,"Channel")==0) return ChannelEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
+              if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
+              else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
+              else if (strcmp(name,"Channel")==0) return ChannelEnum;
+              else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
               else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
               else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
@@ -1364 +1367 @@
               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
               else if (strcmp(name,"Inputs")==0) return InputsEnum;
-              else if (strcmp(name,"Internal")==0) return InternalEnum;
-              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
-              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"J")==0) return JEnum;
+              if (strcmp(name,"Internal")==0) return InternalEnum;
+              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
+              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
+              else if (strcmp(name,"J")==0) return JEnum;
               else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
               else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
@@ -1487 +1490 @@
               else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
               else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
-              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
-              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
-              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
          else stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
+              if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
+              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
+              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
+              else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
               else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
               else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;

issm/trunk-jpl/src/m/classes/stochasticforcing.m ¶

@@ -8 +8 @@
                 isstochasticforcing  = 0;
                 fields               = NaN;
-                dimensions           = NaN;
+                defaultdimension     = 0;
+      default_id           = NaN;
                 covariance           = NaN;
                 randomflag           = 1;
@@ -34 +35 @@
 
                         num_fields = numel(self.fields);
-                        size_tot   = sum(self.dimensions);
+                        
+                        %Check that covariance matrix is positive definite
+                        try
+                                chol(self.covariance);
+                        catch
+                                error('md.stochasticforcing.covariance is not positive definite');
+                        end
+
+                        %Check that all fields agree with the corresponding md class and if any field needs the default params   
+         checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality
+                        for field=self.fields
+            %Checking agreement of classes
+            if(contains(field,'SMB'))
+               if~(isequal(class(md.smb),char(field)))
+                  error('md.smb does not agree with stochasticforcing field %s', char(field));
+               end
+            end
+            if(contains(field,'frontalforcings'))
+               if~(isequal(class(md.frontalforcings),char(field)))
+                  error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
+               end
+            end
+            %Checking for specific dimensions
+            if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
+               checkdefaults = true; %field with non-specific dimensionality
+            end
+         end
+
+                        %Retrieve sum of all the field dimensionalities
+                        size_tot   = self.defaultdimension*num_fields;
+                        indSMBar = -1; %about to check for index of SMBautoregression
+                        indTFar  = -1; %about to check for index of FrontalForcingsRignotAutoregression
+                        if any(contains(self.fields,'SMBautoregression'))
+            size_tot = size_tot-self.defaultdimension+md.smb.num_basins;
+                                indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021)
+         end
+         if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
+            size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins;
+                                indTFar  = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021)
+         end
+
+                        if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
+                                if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar))))~=0)
+                                        error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
+                                end
+                        end
 
                         md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
-                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021)
-                        md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field
+                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
                         md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
                         md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
-
-                        %Check that all fields agree with the corresponding md class
-                        for field=self.fields
-                                if(contains(field,'SMB'))
-                                        if~(isequal(class(md.smb),char(field)))
-                                                error('md.smb does not agree with stochasticforcing field %s', char(field));
-                                        end
-                                end
-                                if(contains(field,'frontalforcings'))
-                                        if~(isequal(class(md.frontalforcings),char(field)))
-                                                error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
-                                        end
-                                end
-                        end
+                        if(checkdefaults) %need to check the defaults
+            md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0);
+            md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]);
+         end
                 end % }}}
                 function disp(self) % {{{
@@ -60 +95 @@
                         fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
                         fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
-                        fielddisplay(self,'dimensions','dimensionality of each field');
+                        fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
+         fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
                         fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
                         fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
@@ -77 +113 @@
                                 return
                         else
+
+                                %Retrieve dimensionality of each field
+                                dimensions = self.defaultdimension*ones(1,num_fields);
+                                ind = 1;
+                                for field=self.fields
+                                        %Checking for specific dimensions
+                                        if(strcmp(field,'SMBautoregression'))
+                                                dimensions(ind) = md.smb.num_basins;
+                                        end
+                                        if(strcmp(field,'FrontalForcingsRignotAutoregression'))
+                                                dimensions(ind) = md.frontalforcings.num_basins;
+                                        end
+                                        ind = ind+1;
+                                end
+
                                 %Scaling covariance matrix (scale column-by-column and row-by-row)
-                                scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts
+                                scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
                                 tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
                                 for i=1:num_fields
                                         if any(strcmp(scaledfields,self.fields(i)))
-                                                inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
+                                                inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
                                                 for row=inds %scale rows corresponding to scaled field
                                                         tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
@@ -93 +144 @@
                                 WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
                                 WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
-                                WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat');
+                                WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat');
+            WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed
+                                WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
                                 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
                                 WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
@@ -105 +158 @@
 
    list = {...
-      'SMBautoregression',...
-      'FrontalForcingsRignotAutoregression'
+      'DefaultCalving',...
+                'FrontalForcingsRignotAutoregression',...
+      'SMBautoregression'
       };
 end % }}}

issm/trunk-jpl/src/m/classes/stochasticforcing.py ¶

@@ -15 +15 @@
     def __init__(self, *args):  # {{{
         self.isstochasticforcing = 0
-        self.fields = np.nan
-        self.dimensions = np.nan
-        self.covariance = np.nan
-        self.randomflag = 1
+        self.fields              = np.nan
+        self.defaultdimension    = 0
+        self.default_id          = np.nan
+        self.covariance          = np.nan
+        self.randomflag          = 1
 
         if len(args) == 0:
@@ -29 +30 @@
         s += '{}\n'.format(fielddisplay(self, 'isstochasticforcing', 'is stochasticity activated?'))
         s += '{}\n'.format(fielddisplay(self, 'fields', 'fields with stochasticity applied, ex: [\'SMBautoregression\'], or [\'FrontalForcingsRignotAutoregression\']'))
+        s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)'))
+        s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'))
         s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)'))
         s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
         s += 'Available fields:\n'
+        s += '   DefaultCalving\n'
         s += '   SMBautoregression\n'
         s += '   FrontalForcingsRignotAutoregression (thermal forcing)\n'
@@ -51 +55 @@
 
         num_fields  = numel(self.fields)
-        size_tot    = np.sum(self.dimensions)
 
-        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
-        md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) # VV check here 'cell' (19Oct2021)
-        md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
-        md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
-        md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
+        #Check that covariance matrix is positive definite
+        try:
+            np.linalg.cholesky(self.covariance);
+        except:
+            error('md.stochasticforcing.covariance is not positive definite');
 
         # Check that all fields agree with the corresponding md class
+        checkdefaults = False
         for field in self.fields:
             if (contains(field, 'SMB')):
@@ -67 +71 @@
                 if not (type(md.frontalforcings) == field):
                     error('md.frontalforcings does not agree with stochasticforcing field {}'.format(field))
+            #Checking for specific dimensions
+            if (field!='SMBautoregression' or field!='FrontalForcingsRignotAutoregression'):
+                checkdefaults = True #field with non-specific dimensionality
+
+        #Retrieve sum of all the field dimensionalities
+        size_tot   = self.defaultdimension*num_fields;
+        indSMBar = -1; #about to check for index of SMBautoregression
+        indTFar  = -1; #about to check for index of FrontalForcingsRignotAutoregression
+        if ('SMBautoregression' in self.fields):
+            size_tot = size_tot-self.defaultdimension+md.smb.num_basins
+            indSMBar = np.where(self.fields=='SMBautoregression')[0][0]
+        if ('FrontalForcingsRignotAutoregression' in self.fields):
+            size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins
+            indSMBar = np.where(self.fields=='FrontalForcingsRignotAutoregression')[0][0]
+        if (indSMBar!=-1 and indTFar!=-1):
+            if((md.smb.ar_timestep!=md.frontalforcings.ar_timestep) and np.any(self.covariance[np.sum(self.dimensions[0:indSMBar]).astype(int):np.sum(self.dimensions[0:indSMBar+1]).astype(int),np.sum(self.dimensions[0:indTFar]).astype(int):np.sum(self.dimensions[0:indTFar+1]).astype(int)]!=0)):
+                error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance')
+
+        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
+        md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings())
+        #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
+        md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
+        md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
+        if (checkdefaults):
+            md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
+            md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements,1])
         return md
     # }}}
@@ -83 +113 @@
             return md
         else:
+            #Retrieve dimensionality of each field
+            dimensions = self.defaultdimension*np.ones(num_fields)
+            for ind,field in enumerate(self.fields):
+                #Checking for specific dimensions
+                if (field=='SMBautoregression'):
+                    dimensions[ind] = md.smb.num_basins
+                if (field=='FrontalForcingsRignotAutoregression'):
+                    dimensions[ind] = md.frontalforcings.num_basins
+
             # Scaling covariance matrix (scale column-by-column and row-by-row)
-            scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts
-            tempcovariance = self.covariance
+            scaledfields = ['DefaultCalving','SMBautoregression'] # list of fields that need scaling * 1/yts
+            tempcovariance = np.copy(self.covariance)
             for i in range(num_fields):
                 if self.fields[i] in scaledfields:
                     inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1])))
                     for row in inds: # scale rows corresponding to scaled field
-                        tempcovariance[row, :] = 1 / yts * self.covariance[row, :]
+                        tempcovariance[row, :] = 1 / yts * tempcovariance[row, :]
                     for col in inds: # scale columns corresponding to scaled field
-                        tempcovariance[:, col] = 1 / yts * self.covariance[:, col]
+                        tempcovariance[:, col] = 1 / yts * tempcovariance[:, col]
             WriteData(fid, prefix, 'data', num_fields, 'name', 'md.stochasticforcing.num_fields', 'format', 'Integer')
             WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray')
-            WriteData(fid, prefix, 'object', self, 'fieldname','dimensions', 'format', 'IntMat')
+            WriteData(fid, prefix, 'data', 'dimensions', 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat')
+            WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed!
+            WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer') 
             WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
             WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
@@ -104 +145 @@
     """
     return [
-        'SMBautoregression',
-        'FrontalForcingsRignotAutoregression'
+        'DefaultCalving',
+        'FrontalForcingsRignotAutoregression',
+        'SMBautoregression'
     ]

issm/trunk-jpl/test/NightlyRun/test257.m ¶

@@ -49 +49 @@
 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

issm/trunk-jpl/test/NightlyRun/test257.py ¶

@@ -56 +56 @@
 md.stochasticforcing.isstochasticforcing = 1
 md.stochasticforcing.fields = ['SMBautoregression']
-md.stochasticforcing.dimensions = [md.smb.num_basins] # dimension of each field
 md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields
 md.stochasticforcing.randomflag = 0 # fixed random seeds

issm/trunk-jpl/test/NightlyRun/test543.m ¶

@@ -39 +39 @@
 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

issm/trunk-jpl/test/NightlyRun/test543.py ¶

@@ -49 +49 @@
 md.stochasticforcing.isstochasticforcing = 1
 md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression']
-md.stochasticforcing.dimensions = [md.frontalforcings.num_basins] # dimension of each field
 md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields
 md.stochasticforcing.randomflag = 0 # determines true/false randomness