Changeset 26526


Ignore:
Timestamp:
11/02/21 07:18:54 (3 years ago)
Author:
vverjans
Message:

CHG: matlab stochasticforcing class, update of SMBautoregression and FrontalForcingsRignotAutoregression, no noise in AutoregressionInit

Location:
issm/trunk-jpl
Files:
2 added
19 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r26477 r26526  
    247247        ./modules/ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.cpp \
    248248        ./modules/Solverx/Solverx.cpp \
     249        ./modules/StochasticForcingx/StochasticForcingx.cpp \
    249250        ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \
    250251        ./cores/ProcessArguments.cpp \

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

    r26495 r26526  
    147147                case FrontalForcingsRignotAutoregressionEnum:
    148148                        /*Retrieve autoregressive parameters*/
    149                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.randomflag",FrontalForcingsRandomflagEnum));
    150149         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsAutoregressiveOrderEnum));
    151150         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_initialtime",FrontalForcingsAutoregressionInitialTimeEnum));
    … …  
    159158         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.phi");
    160159         parameters->AddObject(new DoubleMatParam(FrontalForcingsPhiEnum,transparam,M,N));
    161          xDelete<IssmDouble>(transparam);
    162          iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.covmat");
    163          parameters->AddObject(new DoubleMatParam(FrontalForcingsCovmatEnum,transparam,M,N));
    164160         xDelete<IssmDouble>(transparam);
    165161                        /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/

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

    r26495 r26526  
    6060        return false;
    6161}/*}}}*/
    62 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){/*{{{*/
    63 
    64         const int numvertices = this->GetNumberOfVertices();
     62void       Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type){/*{{{*/
     63
     64const int numvertices = this->GetNumberOfVertices();
    6565   int         basinid;
    66    IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; 
     66   IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;
    6767   IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
    6868   IssmDouble* varspin     = xNewZeroInit<IssmDouble>(numvertices*arorder);
    6969
    7070   /*Get Basin ID and Basin coefficients*/
    71         if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
     71   if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
    7272   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
    7373   for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
    … …  
    7575   beta1_basin   = beta1[basinid];
    7676
    77    /*Loop over number of spin-up steps for all vertices*/   
    78    for(int j=0;j<nspin;j++){ 
    79       tspin = starttime-((nspin-j)*tstep_ar);
    80       noisespin_basin = noisespin[j*numbasins+basinid];
     77   /*Loop over number of spin-up steps for all vertices*/
     78   for(int j=0;j<nspin;j++){
     79      tspin = starttime-((nspin-j)*tstep_ar);
    8180      IssmDouble* oldvarspin = xNew<IssmDouble>(numvertices*arorder);
    8281      for(int i=0;i<numvertices*arorder;i++) oldvarspin[i]=varspin[i];
    … …  
    8584         IssmDouble autoregressionterm = 0.;
    8685         for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*varspin[v+i*numvertices];
    87          varspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin;
     86         varspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm;
    8887      }
    8988
    90       /*Adjustt older values in varspin*/
    91       for(int i=0;i<(arorder-1)*numvertices;i++) varspin[i+numvertices]=oldvarspin[i]; 
    92 
    93       xDelete<IssmDouble>(oldvarspin); 
     89      /*Adjust older values in varspin*/
     90      for(int i=0;i<(arorder-1)*numvertices;i++) varspin[i+numvertices]=oldvarspin[i];
     91
     92      xDelete<IssmDouble>(oldvarspin);
    9493   }
    95         if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
     94   if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
    9695   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
    9796
    98  
    9997   /*Cleanup and return*/
    100    xDelete<IssmDouble>(varspin); 
    101    xDelete<IssmDouble>(phi_basin);
     98   xDelete<IssmDouble>(varspin);
     99   xDelete<IssmDouble>(phi_basin);     
    102100}/*}}}*/
    103101void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
    104    
    105    const int numvertices = this->GetNumberOfVertices();
     102       
     103        const int numvertices = this->GetNumberOfVertices();
    106104   int         basinid,M,N;
    107    IssmDouble  beta0_basin,beta1_basin,noise_basin; 
     105   IssmDouble  beta0_basin,beta1_basin,noise_basin;
    108106   IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    109107   IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
    … …  
    113111   if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
    114112   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
    115    for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii]; 
     113   for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
    116114   beta0_basin   = beta0[basinid];
    117115   beta1_basin   = beta1[basinid];
    … …  
    138136      IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
    139137      /*Assign newest values and shift older values*/
    140       for(int i=0;i<numvertices;i++) temparray[i] = varlist[i]; 
     138      for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
    141139      for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
    142       if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
    143       if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
     140      if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
     141      if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
    144142      xDelete<IssmDouble>(temparray);
    145143   }
    146    
    147         //if(this->lid==55){_printf_("varlist: "<<varlist[0]<<"  "<<varlist[1]<<"  "<<varlist[2]<<"  "<<"noise_basin: "<<noise_basin<<"  "<<'\n');}
    148 
    149         /*Add input to element*/
     144
     145   /*Add input to element*/
    150146   if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
    151147   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r26432 r26526  
    127127        /*parameters: */
    128128        bool isstressbalance,ismasstransport,isoceantransport,issmb,isthermal,isgroundingline,isesa,issampling;;
    129         bool isslc,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,save_results;
     129        bool isslc,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,isstochasticforcing,save_results;
    130130        int  step,sb_coupling_frequency;
    131131        int  domaintype,numoutputs;
    … …  
    150150        femmodel->parameters->FindParam(&issampling,TransientIssamplingEnum);
    151151        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
     152        femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
    152153
    153154        #if defined(_HAVE_OCEAN_ )
    154155        if(isoceancoupling) OceanExchangeDatax(femmodel,false);
    155156        #endif
     157
     158        if(isstochasticforcing){
     159                StochasticForcingx(femmodel);
     160        }
    156161
    157162        if(isthermal && domaintype==Domain3DEnum){

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

    r26495 r26526  
    3939   IssmDouble* beta1    = NULL;
    4040   IssmDouble* phi      = NULL;
    41    IssmDouble* covmat   = NULL;
    4241   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    4342   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
    … …  
    4645   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
    4746   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    48    femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
    4947
    50    /*AR model spin-up*/
    51    int nspin{2*arorder+5};
    52    IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin);
    53    my_rank=IssmComm::GetRank();
    54    if(my_rank==0){
    55       bool randomflag{};
    56       femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum);
    57       int fixedseed;
    58       for(int i=0;i<nspin;i++){
    59          IssmDouble* temparray = NULL;
    60          /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
    61          if(randomflag) fixedseed=-1;
    62          else fixedseed = i;
    63          multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
    64          for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
    65          xDelete<IssmDouble>(temparray);
    66       }
    67    }
    68    ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    69 
    70    /*Initialize ThermalforcingValuesAutoregressionEnum*/
     48   /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum*/
     49        int nspin{2*arorder+5};
    7150   for(Object* &object:femmodel->elements->objects){
    7251      Element* element      = xDynamicCast<Element*>(object); //generate element object
    73       element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,FrontalForcingsRignotAutoregressionEnum);
     52      element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum);
    7453   }
    7554
    … …  
    7857   xDelete<IssmDouble>(beta1);
    7958   xDelete<IssmDouble>(phi);
    80    xDelete<IssmDouble>(noisespin);
    81    xDelete<IssmDouble>(covmat);
    8259}/*}}}*/
    8360void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/
    … …  
    10279
    10380   /*Load parameters*/
     81        bool isstochastic;
    10482   int M,N,Nphi,arorder,numbasins,my_rank;
    10583   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    … …  
    10987   IssmDouble* beta1      = NULL;
    11088   IssmDouble* phi        = NULL;
    111    IssmDouble* covmat     = NULL;
    112    IssmDouble* noiseterms = xNew<IssmDouble>(numbasins);
    113    femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
     89   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
     90
     91        femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
    11492   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
    11593   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
    11694   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    117    femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
    11895
     96        /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
     97        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
     98        if(isstochastic){
     99                int  numstochasticfields;
     100                int* stochasticfields;
     101                femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
     102                femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
     103                for(int i=0;i<numstochasticfields;i++){
     104                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
     105                                femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum);  _assert_(M==numbasins);
     106                        }
     107                }
     108        }
    119109   /*Time elapsed with respect to AR model initial time*/
    120110   IssmDouble telapsed_ar = time-tinit_ar;
    121111
    122    /*Before looping through elements: compute noise term specific to each basin from covmat*/
    123    my_rank=IssmComm::GetRank();
    124    if(my_rank==0){
    125       bool randomflag{};
    126       femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum);
    127       /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
    128       int fixedseed;
    129       if(randomflag) fixedseed=-1;
    130       else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
    131                 /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
    132       IssmDouble* temparray = NULL;
    133       multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
    134       for(int i=0;i<numbasins;i++) noiseterms[i]=temparray[i];
    135       xDelete<IssmDouble>(temparray);
    136    }
    137    ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    138 
    139    /*Loop over each element to compute SMB at vertices*/
     112   /*Loop over each element to compute Thermal Forcing at vertices*/
    140113   for(Object* &object:femmodel->elements->objects){
    141114      Element* element = xDynamicCast<Element*>(object);
    … …  
    148121   xDelete<IssmDouble>(phi);
    149122   xDelete<IssmDouble>(noiseterms);
    150    xDelete<IssmDouble>(covmat);
    151123}/*}}}*/

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

    r26483 r26526  
    9090                parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));
    9191                parameters->AddObject(iomodel->CopyConstantObject("md.transient.issampling",TransientIssamplingEnum));
     92                parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.isstochasticforcing",StochasticForcingIsStochasticForcingEnum));
    9293
    9394                /*For stress balance only*/
    … …  
    399400         /*Add parameters that are not in standard nbvertices format*/
    400401         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum));
    401          parameters->AddObject(iomodel->CopyConstantObject("md.smb.randomflag",SmbRandomflagEnum));
    402402         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbAutoregressiveOrderEnum));
    403403         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum));
    … …  
    412412         parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
    413413         xDelete<IssmDouble>(transparam);
    414          iomodel->FetchData(&transparam,&M,&N,"md.smb.covmat");
    415          parameters->AddObject(new DoubleMatParam(SmbCovmatEnum,transparam,M,N));
    416          xDelete<IssmDouble>(transparam);
    417414         break;
    418415                case SMBgembEnum:
    … …  
    500497                if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
    501498                iomodel->DeleteData(&requestedoutputs,numoutputs,"md.damage.requested_outputs");
     499        }
     500
     501        bool isstochasticforcing;
     502        parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
     503        if(isstochasticforcing){
     504                int num_fields;
     505                char** fields;
     506                parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
     507                iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
     508                if(num_fields<1) _error_("no stochasticforcing fields found");
     509                int* stochasticforcing_enums = xNew<int>(num_fields);
     510                for(int i=0;i<num_fields;i++){
     511                        stochasticforcing_enums[i] = StringToEnumx(fields[i]);
     512                        xDelete<char>(fields[i]);
     513                }
     514                xDelete<char*>(fields);
     515                parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
     516                xDelete<int>(stochasticforcing_enums);
     517
     518      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum));
     519                iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
     520                parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
     521                xDelete<IssmDouble>(transparam);
     522                iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
     523                parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
     524                xDelete<IssmDouble>(transparam);
    502525        }
    503526

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

    r26495 r26526  
    157157        IssmDouble* beta1    = NULL;
    158158        IssmDouble* phi      = NULL;
    159         IssmDouble* covmat   = NULL;
    160159        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    161160   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
    … …  
    164163        femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
    165164        femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    166         femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
    167165       
    168         /*AR model spin-up*/
    169         int nspin{2*arorder+5};
    170         IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin);
    171         my_rank=IssmComm::GetRank();
    172         if(my_rank==0){
    173                 bool randomflag{};
    174                 femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
    175                 int fixedseed;
    176                 for(int i=0;i<nspin;i++){
    177                         //IssmDouble* temparray = xNew<IssmDouble>(numbasins);
    178                         IssmDouble* temparray = NULL;
    179                         /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
    180                         if(randomflag) fixedseed=-1;
    181                         else fixedseed = i;
    182                         multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
    183                         for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
    184                         xDelete<IssmDouble>(temparray);
    185                 }
    186         }
    187         ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    188        
    189         /*Initialize SmbValuesAutoregressionEnum*/
     166        /*AR model spin-up with 0 noise to initialize SmbValuesAutoregressionEnum*/
     167        int nspin{2*arorder+5};
    190168        for(Object* &object:femmodel->elements->objects){
    191169      Element* element      = xDynamicCast<Element*>(object); //generate element object
    192                 element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,SMBautoregressionEnum);
    193         }
    194        
     170                element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,SMBautoregressionEnum);
     171        }
    195172        /*Cleanup*/
    196173        xDelete<IssmDouble>(beta0);
    197174        xDelete<IssmDouble>(beta1);
    198175        xDelete<IssmDouble>(phi);
    199         xDelete<IssmDouble>(noisespin);
    200         xDelete<IssmDouble>(covmat);
    201176}/*}}}*/
    202177void Smbautoregressionx(FemModel* femmodel){/*{{{*/
    … …  
    221196
    222197        /*Load parameters*/
     198        bool isstochastic;
    223199        int M,N,Nphi,arorder,numbasins,my_rank;
    224200        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
    … …  
    228204        IssmDouble* beta1      = NULL;
    229205        IssmDouble* phi        = NULL;
    230         IssmDouble* covmat     = NULL;
    231         IssmDouble* noiseterms = xNew<IssmDouble>(numbasins);
     206
    232207        femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
    233208        femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
    234209        femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
    235210        femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    236         femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
    237 
     211
     212        /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
     213   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
     214        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
     215   if(isstochastic){
     216                int  numstochasticfields;
     217      int* stochasticfields;
     218      femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
     219      femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
     220      for(int i=0;i<numstochasticfields;i++){
     221         if(stochasticfields[i]==SMBautoregressionEnum){
     222                                femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum);  _assert_(M==numbasins);
     223                        }
     224                }
     225                xDelete<int>(stochasticfields);
     226        }
    238227        /*Time elapsed with respect to AR model initial time*/
    239228        IssmDouble telapsed_ar = time-tinit_ar;
    240 
    241         /*Before looping through elements: compute noise term specific to each basin from covmat*/
    242         my_rank=IssmComm::GetRank();
    243         if(my_rank==0){
    244                 bool randomflag{};
    245                 femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum);
    246                 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
    247                 int fixedseed;
    248                 if(randomflag) fixedseed=-1;
    249                 else fixedseed =  reCast<int,IssmDouble>((time-starttime)/dt);
    250                 /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/   
    251                 IssmDouble* temparray = NULL;
    252                 multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
    253                 for(int i=0;i<numbasins;i++) noiseterms[i]=temparray[i];
    254                 xDelete<IssmDouble>(temparray);
    255         }
    256         ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    257229
    258230        /*Loop over each element to compute SMB at vertices*/
    … …  
    267239        xDelete<IssmDouble>(phi);
    268240        xDelete<IssmDouble>(noiseterms);
    269         xDelete<IssmDouble>(covmat);
    270241}/*}}}*/
    271242void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/

  • issm/trunk-jpl/src/c/modules/modules.h

    r26099 r26526  
    9191#include "./Scotchx/Scotchx.h"
    9292#include "./Shp2Kmlx/Shp2Kmlx.h"
     93#include "./StochasticForcingx/StochasticForcingx.h"
    9394#include "./SurfaceMassBalancex/SurfaceMassBalancex.h"
    9495#include "./Solverx/Solverx.h"

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

    r26495 r26526  
    182182syn keyword cConstant FrontalForcingsBeta0Enum
    183183syn keyword cConstant FrontalForcingsBeta1Enum
    184 syn keyword cConstant FrontalForcingsCovmatEnum
    185184syn keyword cConstant FrontalForcingsNumberofBasinsEnum
    186185syn keyword cConstant FrontalForcingsParamEnum
    187186syn keyword cConstant FrontalForcingsPhiEnum
    188 syn keyword cConstant FrontalForcingsRandomflagEnum
    189187syn keyword cConstant GrdModelEnum
    190188syn keyword cConstant GroundinglineFrictionInterpolationEnum
    … …  
    399397syn keyword cConstant SolidearthSettingsGrdOceanEnum
    400398syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
     399syn keyword cConstant StochasticForcingCovarianceEnum
     400syn keyword cConstant StochasticForcingDimensionsEnum
     401syn keyword cConstant StochasticForcingFieldsEnum
     402syn keyword cConstant StochasticForcingIsStochasticForcingEnum
     403syn keyword cConstant StochasticForcingNumFieldsEnum
     404syn keyword cConstant StochasticForcingRandomflagEnum
    401405syn keyword cConstant RotationalPolarMoiEnum
    402406syn keyword cConstant SolidearthSettingsReltolEnum
    … …  
    424428syn keyword cConstant SmbAdThreshEnum
    425429syn keyword cConstant SmbAutoregressionInitialTimeEnum
     430syn keyword cConstant SmbAutoregressionNoiseEnum
    426431syn keyword cConstant SmbAutoregressionTimestepEnum
    427432syn keyword cConstant SmbAutoregressiveOrderEnum
    … …  
    429434syn keyword cConstant SmbBeta0Enum
    430435syn keyword cConstant SmbBeta1Enum
    431 syn keyword cConstant SmbCovmatEnum
    432436syn keyword cConstant SmbDesfacEnum
    433437syn keyword cConstant SmbDpermilEnum
    … …  
    462466syn keyword cConstant SmbPfacEnum
    463467syn keyword cConstant SmbPhiEnum
    464 syn keyword cConstant SmbRandomflagEnum
    465468syn keyword cConstant SmbRdlEnum
    466469syn keyword cConstant SmbRequestedOutputsEnum
    … …  
    498501syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
    499502syn keyword cConstant StressbalanceShelfDampeningEnum
     503syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
    500504syn keyword cConstant ThermalIsdrainicecolumnEnum
    501505syn keyword cConstant ThermalIsdynamicbasalspcEnum

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

    r26495 r26526  
    176176        FrontalForcingsBeta0Enum,
    177177   FrontalForcingsBeta1Enum,
    178    FrontalForcingsCovmatEnum,
    179178        FrontalForcingsNumberofBasinsEnum,
    180179        FrontalForcingsParamEnum,
    181180   FrontalForcingsPhiEnum,
    182         FrontalForcingsRandomflagEnum,
    183181        GrdModelEnum,
    184182        GroundinglineFrictionInterpolationEnum,
    … …  
    393391        SolidearthSettingsGrdOceanEnum,
    394392        SolidearthSettingsOceanAreaScalingEnum,
     393        StochasticForcingCovarianceEnum,
     394        StochasticForcingDimensionsEnum,
     395        StochasticForcingFieldsEnum,
     396        StochasticForcingIsStochasticForcingEnum,
     397        StochasticForcingNumFieldsEnum,
     398        StochasticForcingRandomflagEnum,
    395399        RotationalPolarMoiEnum,
    396400        SolidearthSettingsReltolEnum,
    … …  
    418422        SmbAdThreshEnum,
    419423        SmbAutoregressionInitialTimeEnum,
    420    SmbAutoregressionTimestepEnum,
     424   SmbAutoregressionNoiseEnum,
     425        SmbAutoregressionTimestepEnum,
    421426   SmbAutoregressiveOrderEnum,
    422427        SmbAveragingEnum,
    423428        SmbBeta0Enum,
    424429   SmbBeta1Enum,
    425    SmbCovmatEnum,
    426430        SmbDesfacEnum,
    427431        SmbDpermilEnum,
    … …  
    456460        SmbPfacEnum,
    457461        SmbPhiEnum,
    458         SmbRandomflagEnum,
    459462        SmbRdlEnum,
    460463        SmbRequestedOutputsEnum,
    … …  
    492495        StressbalanceRiftPenaltyThresholdEnum,
    493496        StressbalanceShelfDampeningEnum,
     497   ThermalforcingAutoregressionNoiseEnum,
    494498        ThermalIsdrainicecolumnEnum,
    495499        ThermalIsdynamicbasalspcEnum,

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

    r26495 r26526  
    184184                case FrontalForcingsBeta0Enum : return "FrontalForcingsBeta0";
    185185                case FrontalForcingsBeta1Enum : return "FrontalForcingsBeta1";
    186                 case FrontalForcingsCovmatEnum : return "FrontalForcingsCovmat";
    187186                case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
    188187                case FrontalForcingsParamEnum : return "FrontalForcingsParam";
    189188                case FrontalForcingsPhiEnum : return "FrontalForcingsPhi";
    190                 case FrontalForcingsRandomflagEnum : return "FrontalForcingsRandomflag";
    191189                case GrdModelEnum : return "GrdModel";
    192190                case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation";
    … …  
    401399                case SolidearthSettingsGrdOceanEnum : return "SolidearthSettingsGrdOcean";
    402400                case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
     401                case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
     402                case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
     403                case StochasticForcingFieldsEnum : return "StochasticForcingFields";
     404                case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
     405                case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
     406                case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
    403407                case RotationalPolarMoiEnum : return "RotationalPolarMoi";
    404408                case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";
    … …  
    426430                case SmbAdThreshEnum : return "SmbAdThresh";
    427431                case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
     432                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
    428433                case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
    429434                case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
    … …  
    431436                case SmbBeta0Enum : return "SmbBeta0";
    432437                case SmbBeta1Enum : return "SmbBeta1";
    433                 case SmbCovmatEnum : return "SmbCovmat";
    434438                case SmbDesfacEnum : return "SmbDesfac";
    435439                case SmbDpermilEnum : return "SmbDpermil";
    … …  
    464468                case SmbPfacEnum : return "SmbPfac";
    465469                case SmbPhiEnum : return "SmbPhi";
    466                 case SmbRandomflagEnum : return "SmbRandomflag";
    467470                case SmbRdlEnum : return "SmbRdl";
    468471                case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
    … …  
    500503                case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
    501504                case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
     505                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
    502506                case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
    503507                case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";

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

    r26495 r26526  
    187187              else if (strcmp(name,"FrontalForcingsBeta0")==0) return FrontalForcingsBeta0Enum;
    188188              else if (strcmp(name,"FrontalForcingsBeta1")==0) return FrontalForcingsBeta1Enum;
    189               else if (strcmp(name,"FrontalForcingsCovmat")==0) return FrontalForcingsCovmatEnum;
    190189              else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum;
    191190              else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum;
    192191              else if (strcmp(name,"FrontalForcingsPhi")==0) return FrontalForcingsPhiEnum;
    193               else if (strcmp(name,"FrontalForcingsRandomflag")==0) return FrontalForcingsRandomflagEnum;
    194192              else if (strcmp(name,"GrdModel")==0) return GrdModelEnum;
    195193              else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum;
    … …  
    260258              else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
    261259              else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
     260              else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
     261              else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
    266               else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    267               else if (strcmp(name,"InversionType")==0) return InversionTypeEnum;
     265              if (strcmp(name,"InversionType")==0) return InversionTypeEnum;
    268266              else if (strcmp(name,"Ivins")==0) return IvinsEnum;
    269267              else if (strcmp(name,"IsSlcCoupling")==0) return IsSlcCouplingEnum;
    … …  
    383381              else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
    384382              else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
     383              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
     384              else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
    389               else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
    390               else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
     388              if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
    391389              else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
    392390              else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum;
    … …  
    410408              else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
    411409              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
     410              else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
     411              else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
     412              else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
     413              else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
     414              else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
     415              else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
    412416              else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum;
    413417              else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    … …  
    435439              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
    436440              else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
     441              else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
    437442              else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
    438443              else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
    … …  
    440445              else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum;
    441446              else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum;
    442               else if (strcmp(name,"SmbCovmat")==0) return SmbCovmatEnum;
    443447              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    444448              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    … …  
    473477              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    474478              else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
    475               else if (strcmp(name,"SmbRandomflag")==0) return SmbRandomflagEnum;
    476479              else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
    477480              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    … …  
    503506              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    504507              else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
    505               else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
    506               else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
    507               else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
     511              if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
     512              else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
     513              else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
     514              else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
    512515              else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
    513516              else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
     517              else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
    514518              else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
    515519              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    … …  
    625629              else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
    626630              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    627               else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
    628635              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    629636              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    630637              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
     638              else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
    635639              else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum;
    636640              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    … …  
    748752              else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
    749753              else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum;
    750               else if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"HydrologyTwsAnalysis")==0) return HydrologyTwsAnalysisEnum;
    751758              else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum;
    752759              else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    753760              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Ice")==0) return IceEnum;
     761              else if (strcmp(name,"Ice")==0) return IceEnum;
    758762              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    759763              else if (strcmp(name,"Input")==0) return InputEnum;
    … …  
    871875              else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    872876              else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
    873               else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
    874881              else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
    875882              else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
    876883              else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
     884              else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
    881885              else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum;
    882886              else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum;
    … …  
    994998              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    995999              else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
    996               else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    9971004              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
    9981005              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    9991006              else if (strcmp(name,"SmbValuesAutoregression")==0) return SmbValuesAutoregressionEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmbV")==0) return SmbVEnum;
     1007              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
    10041008              else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
    10051009              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    … …  
    11171121              else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
    11181122              else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
    1119               else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
    11201127              else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
    11211128              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
    11221129              else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
     1130              else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
    11271131              else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
    11281132              else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
    … …  
    12401244              else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
    12411245              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    1242               else if (strcmp(name,"Colinear")==0) return ColinearEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Colinear")==0) return ColinearEnum;
    12431250              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    12441251              else if (strcmp(name,"Contact")==0) return ContactEnum;
    12451252              else if (strcmp(name,"Contour")==0) return ContourEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Contours")==0) return ContoursEnum;
     1253              else if (strcmp(name,"Contours")==0) return ContoursEnum;
    12501254              else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    12511255              else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum;
    … …  
    13631367              else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
    13641368              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    1365               else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    13661373              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    13671374              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    13681375              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
     1376              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    13731377              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    13741378              else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
    … …  
    14861490              else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
    14871491              else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
    1488               else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    14891496              else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
    14901497              else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
    14911498              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
     1499              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    14961500              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    14971501              else if (strcmp(name,"Scaled")==0) return ScaledEnum;

  • issm/trunk-jpl/src/m/classes/SMBautoregression.m

    r26486 r26526  
    1313                ar_timestep       = 0;
    1414                phi               = NaN;
    15                 covmat            = NaN;
    1615                basin_id          = NaN;
    17                 randomflag        = 1;
    1816                steps_per_step    = 1;
    1917                averaging         = 0;
    … …  
    5755                                disp('      smb.phi (lag coefficients) not specified: order of autoregressive model set to 0');
    5856                        end
    59                         if isnan(self.covmat)
    60                                 self.covmat = 1e-21*eye(self.num_basins); %no stochasticity and no covariance
    61                                 disp('      smb.covmat not specified: stochasticity set to 0');
    62                         end
    6357                end % }}}
    6458                function self = setdefaultparameters(self) % {{{
    6559                        self.ar_order    = 0.0; %autoregression model of order 0
    66                         self.randomflag  = 1;
    6760                end % }}}
    6861                function md = checkconsistency(self,md,solution,analyses)
    … …  
    7770                                md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
    7871                                md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
    79                                 md = checkfield(md,'fieldname','smb.covmat','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.num_basins]);
    80                                 md = checkfield(md,'fieldname','smb.randomflag','numel',[1],'values',[0 1]);
    8172                        end
    8273                        md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
    … …  
    9485                        fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    9586                        fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
    96                         fielddisplay(self,'covmat','inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]');
    97                         fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    9887                        fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
    9988                        fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
    … …  
    117106                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts);
    118107                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts);
    119                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',yts);
    120                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','randomflag','format','Boolean');
    121108                        WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
    122109                        WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');

  • issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m

    r26495 r26526  
    1313      ar_timestep          = 0;
    1414      phi                  = NaN;
    15       covmat               = NaN;
    1615      basin_id             = NaN;
    17       randomflag           = 1;
    1816                subglacial_discharge = NaN;
    1917        end
    … …  
    4644                        subglacial_discharge = NaN;
    4745                        self.ar_order        = 0.0; %autoregression model of order 0
    48          self.randomflag      = 1;
    4946
    5047                end % }}}
    … …  
    6259         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
    6360                        md = checkfield(md,'fieldname','frontalforcings.phi','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]);
    64          md = checkfield(md,'fieldname','frontalforcings.covmat','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.num_basins]);
    65          md = checkfield(md,'fieldname','frontalforcings.randomflag','numel',[1],'values',[0 1]);
    6661
    6762                end % }}}
    … …  
    7772         fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    7873         fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
    79          fielddisplay(self,'covmat','inter-basin covariance matrix for multivariate normal noise at each time step [∘C^2]');
    80          fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    8174                end % }}}
    8275                function marshall(self,prefix,md,fid) % {{{
    … …  
    9285         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta1','format','DoubleMat','name','md.frontalforcings.beta1','scale',1./yts,'yts',md.constants.yts);
    9386         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',md.constants.yts);
    94          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','covmat','format','DoubleMat','name','md.frontalforcings.covmat');
    95          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','randomflag','format','Boolean');
    9687                end % }}}
    9788        end

  • issm/trunk-jpl/src/m/classes/model.m

    r26364 r26526  
    5555                miscellaneous    = 0;
    5656                private          = 0;
     57                stochasticforcing= 0;
    5758
    5859                %}}}
    … …  
    191192                        %2021 February 17
    192193                        if isa(md.sampling,'double'); md.sampling=sampling(); end
     194                        %VV
     195                        if ~isa(md.stochasticforcing,'stochasticforcing'); md.stochasticforcing=stochasticforcing(); end
    193196                end% }}}
    194197        end
    … …  
    250253                        disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
    251254                        disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
     255                        disp(sprintf('%19s: %-22s -- %s','stochasticforcing',['[1x1 ' class(self.stochasticforcing) ']'],'stochasticity applied to model forcings'));
    252256                end % }}}
    253257                function md = setdefaultparameters(md,planet) % {{{
    … …  
    297301                        md.miscellaneous    = miscellaneous();
    298302                        md.private          = private();
     303                        md.stochasticforcing= stochasticforcing();
    299304                end
    300305                %}}}

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

    r26495 r26526  
    4545md.smb.ar_timestep    = 2.0; %timestep of the autoregressive model [yr]
    4646md.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]];
    47 md.smb.covmat         = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]];
    48 md.smb.randomflag     = 0; %fixed random seeds
     47
     48%Stochastic forcing
     49md.stochasticforcing.isstochasticforcing = 1;
     50md.stochasticforcing.fields              = [{'SMBautoregression'}];
     51md.stochasticforcing.dimensions          = [md.smb.num_basins]; %dimension of each field
     52md.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
     53md.stochasticforcing.randomflag          = 0; %fixed random seeds
    4954
    5055md=solve(md,'Transient');

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

    r26495 r26526  
    3535md.frontalforcings.ar_timestep          = 2; %timestep of the autoregressive model [yr]
    3636md.frontalforcings.phi                  = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters
    37 md.frontalforcings.covmat               = 1e-4*[[1.5,0.5];[0.5,0.4]];
    38 md.frontalforcings.randomflag           = 0;
     37
     38%Stochastic forcing
     39md.stochasticforcing.isstochasticforcing = 1;
     40md.stochasticforcing.fields              = [{'FrontalForcingsRignotAutoregression'}];
     41md.stochasticforcing.dimensions          = [md.frontalforcings.num_basins]; %dimension of each field
     42md.stochasticforcing.covariance          = 1e-4*[[1.5,0.5];[0.5,0.4]];; %global covariance among- and between-fields
     43md.stochasticforcing.randomflag          = 0; %determines true/false randomness
    3944
    4045md.transient.ismovingfront = 1;
Note: See TracChangeset for help on using the changeset viewer.