Changeset 26477


Ignore:
Timestamp:
10/13/21 07:11:29 (3 years ago)
Author:
vverjans
Message:

autoregression SMB module, Random utilities, Cholesky decomposition

Location:
issm/trunk-jpl/src
Files:
4 added
13 edited

Legend:

Unmodified
Added
Removed

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

    r26476 r26477  
    187187        ./shared/Elements/DrainageFunctionWaterfraction.cpp \
    188188        ./shared/Elements/EstarComponents.cpp \
     189        ./shared/Random/random.cpp \
    189190        ./shared/String/DescriptorIndex.cpp \
    190191        ./toolkits/issm/IssmToolkitUtils.cpp \

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

    r26468 r26477  
    179179                        iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
    180180                        break;
     181                case SMBautoregressionEnum:
     182         iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
     183         break;
    181184                case SMBhenningEnum:
    182185                        iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
    … …  
    375378                        /*Nothing to add to parameters*/
    376379                        break;
     380                case SMBautoregressionEnum:
     381         /*Nothing to add to parameters*/
     382         break;
    377383                case SMBhenningEnum:
    378384                        /*Nothing to add to parameters*/
    … …  
    467473                        SmbGradientsElax(femmodel);
    468474                        break;
     475                case SMBautoregressionEnum:
     476         if(VerboseSolution())_printf0_("    call smb autoregression module\n");
     477         Smbautoregressionx(femmodel);
     478         break;
    469479                case SMBhenningEnum:
    470480                        if(VerboseSolution())_printf0_("  call smb Henning module\n");

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

    r26468 r26477  
    35743574}
    35753575/*}}}*/
     3576void       Element::SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin){/*{{{*/
     3577        const int numvertices = this->GetNumberOfVertices();
     3578        int         basinid;
     3579        IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars
     3580        IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
     3581   IssmDouble* smbspin                 = xNew<IssmDouble>(numvertices*arorder);
     3582        this->GetInputValue(&basinid,SmbBasinsIdEnum);
     3583        for(int ii{0};ii<arorder;ii++){phi_basin[ii] = phi[basinid*arorder+ii];}
     3584        beta0_basin   = beta0[basinid];
     3585        beta1_basin   = beta1[basinid];
     3586        for(int jj{0};jj<nspin;jj++){   
     3587                tspin = starttime-((nspin-jj)*tstep_ar); //current time in spin-up loop
     3588      noisespin_basin = noisespin[jj*numbasins+basinid];
     3589      IssmDouble* oldsmbspin = xNewZeroInit<IssmDouble>(numvertices*arorder);
     3590      for(int ii{0};ii<numvertices*arorder;ii++){oldsmbspin[ii]=smbspin[ii];} //copy smbspin in oldsmbspin
     3591      for(int v=0;v<numvertices;v++){
     3592         double autoregressionterm{0.0};
     3593         for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];}
     3594         smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
     3595      }
     3596      for(int ii{0};ii<(arorder-1)*numvertices;ii++){smbspin[ii+numvertices]=oldsmbspin[ii];} //correct older values in smbspin
     3597      xDelete<IssmDouble>(oldsmbspin); //cleanup
     3598        }
     3599        this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);//save spin-up autoregression values
     3600   xDelete<IssmDouble>(smbspin); //cleanup
     3601   xDelete<IssmDouble>(phi_basin); //cleanup
     3602}/*}}}*/
     3603void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
     3604        const int numvertices = this->GetNumberOfVertices();
     3605        int         basinid,M,N;
     3606        IssmDouble  beta0_basin,beta1_basin,noise_basin; //initialize scalars
     3607        IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
     3608   IssmDouble* smblist                 = xNew<IssmDouble>(numvertices);
     3609        IssmDouble* smbvaluesautoregression = NULL; //array for past SMB values that we are about to retrieve
     3610        this->GetInputValue(&basinid,SmbBasinsIdEnum);
     3611        for(int ii{0};ii<arorder;ii++){phi_basin[ii] = phi[basinid*arorder+ii];}
     3612        beta0_basin   = beta0[basinid];
     3613        beta1_basin   = beta1[basinid];
     3614        noise_basin   = noiseterms[basinid]; //note that noiseterms is computed at every timestep
     3615        this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); //get array of past SMB values to compute AR model
     3616        /*If not AR model timestep: take the old SMB values*/
     3617        if(isstepforar==false){
     3618                //VV do something with the lapse rate here if needed (12Oct2021)
     3619      for(int ii{0};ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
     3620        }
     3621        /*If AR model timestep: compute SMB values on vertices from AR*/
     3622        else{
     3623                for(int v=0;v<numvertices;v++){
     3624                        double autoregressionterm{0.0};
     3625                        for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];} //compute autoregressive term
     3626                        smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; //stochastic SMB value
     3627                }
     3628                /*Update autoregression smb values*/
     3629                IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
     3630                for(int ii{0};ii<numvertices;ii++){temparray[ii] = smblist[ii];} //first store newly computed smb values
     3631                for(int ii{0};ii<(arorder-1)*numvertices;ii++){temparray[ii+numvertices] = smbvaluesautoregression[ii];} //second shift older smb values
     3632                this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values
     3633                xDelete<IssmDouble>(temparray); //cleanup
     3634        }
     3635        /*Add input to element*/
     3636        this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
     3637        /*Cleanup*/
     3638        xDelete<IssmDouble>(phi_basin);
     3639        xDelete<IssmDouble>(smblist);
     3640        xDelete<IssmDouble>(smbvaluesautoregression);
     3641}/*}}}*/
    35763642void       Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/
    35773643

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

    r26469 r26477  
    176176                void               SmbSemic();
    177177                int                Sid();
     178                void               SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin);
     179                void               Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms);
    178180                void               SmbGemb(IssmDouble timeinputs, int count);
    179181                void               StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

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

    r26468 r26477  
    396396                        /*Nothing to add*/
    397397                        break;
     398                case SMBautoregressionEnum:
     399         /*Add parameters that are not in standard nbvertices format*/
     400         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum));
     401         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbAutoregressiveOrderEnum));
     402         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum));
     403         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum));
     404         iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0");
     405         parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N));
     406         xDelete<IssmDouble>(transparam);
     407         iomodel->FetchData(&transparam,&M,&N,"md.smb.beta1");
     408         parameters->AddObject(new DoubleVecParam(SmbBeta1Enum,transparam,N));
     409         xDelete<IssmDouble>(transparam);
     410         iomodel->FetchData(&transparam,&M,&N,"md.smb.phi");
     411         parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
     412         xDelete<IssmDouble>(transparam);
     413         iomodel->FetchData(&transparam,&M,&N,"md.smb.covmat");
     414         parameters->AddObject(new DoubleMatParam(SmbCovmatEnum,transparam,M,N));
     415         xDelete<IssmDouble>(transparam);
     416         break;
    398417                case SMBgembEnum:
    399418                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));

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

    r26208 r26477  
    88#include "../modules.h"
    99#include "../../classes/Inputs/TransientInput.h"
     10#include "../../shared/Random/random.h"
    1011
    1112void SmbForcingx(FemModel* femmodel){/*{{{*/
    … …  
    144145        }
    145146
     147}/*}}}*/
     148void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
     149        /*Initialization step of Smbautoregressionx*/
     150        int M,N,Nphi,arorder,numbasins;
     151        IssmDouble starttime,tstep_ar,tinit_ar;
     152        IssmDouble* beta0    = xNew<IssmDouble>(numbasins);
     153        IssmDouble* beta1    = xNew<IssmDouble>(numbasins);
     154        IssmDouble* phi      = xNew<IssmDouble>(numbasins*arorder);
     155        IssmDouble* covmat   = xNew<IssmDouble>(numbasins*numbasins);
     156        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     157   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
     158        femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
     159        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
     160        femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
     161        femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
     162        femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
     163        femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
     164        femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
     165        /*AR model spin-up*/
     166        int nspin{2*arorder+5}; //number of spin-up steps to be executed
     167        IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); //sample of basin-specific noise at each spinup step
     168        for(int ii{0};ii<nspin;ii++){
     169                IssmDouble* temparray = xNew<IssmDouble>(numbasins);
     170                multivariateNormal(&temparray,numbasins,0.0,covmat);
     171                for(int jj{0};jj<numbasins;jj++){noisespin[ii*numbasins+jj]=temparray[jj];}
     172                xDelete<IssmDouble>(temparray);
     173        }
     174        /*Initialize SmbValuesAutoregressionEnum*/
     175        for(Object* &object:femmodel->elements->objects){
     176      Element* element      = xDynamicCast<Element*>(object); //generate element object
     177                element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
     178        }
     179        /*Cleanup*/
     180        xDelete<IssmDouble>(beta0);
     181        xDelete<IssmDouble>(beta1);
     182        xDelete<IssmDouble>(phi);
     183        xDelete<IssmDouble>(noisespin);
     184        xDelete<IssmDouble>(covmat);
     185}/*}}}*/
     186void Smbautoregressionx(FemModel* femmodel){/*{{{*/
     187        /*Get time parameters*/
     188        IssmDouble time,dt,starttime,tstep_ar;
     189        femmodel->parameters->FindParam(&time,TimeEnum);
     190        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     191        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     192   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
     193        /*Initialize module at first time step*/
     194        if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
     195        /*Determine if this is a time step for the AR model*/
     196        bool isstepforar{false};
     197   if((std::fmod(time,tstep_ar)<std::fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
     198        /*Load parameters*/
     199        int M,N,Nphi,arorder,numbasins;
     200        femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
     201        femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
     202        IssmDouble tinit_ar;
     203        IssmDouble* beta0                = xNew<IssmDouble>(numbasins);
     204        IssmDouble* beta1                = xNew<IssmDouble>(numbasins);
     205        IssmDouble* phi                  = xNew<IssmDouble>(numbasins*arorder);
     206        IssmDouble* noiseterms           = xNew<IssmDouble>(numbasins);
     207        IssmDouble* covmat               = xNew<IssmDouble>(numbasins*numbasins);
     208        femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
     209        femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
     210        femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
     211        femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
     212        femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
     213        IssmDouble telapsed_ar = time-tinit_ar; //time elapsed with respect to AR model initial time
     214        /*Before looping through elements: compute noise term specific to each basin from covmat*/
     215        multivariateNormal(&noiseterms,numbasins,0.0,covmat);
     216        /*Loop over each element to compute SMB at vertices*/
     217        for(Object* &object:femmodel->elements->objects){
     218                Element* element      = xDynamicCast<Element*>(object); //generate element object
     219                element->Smbautoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms);
     220        }
     221        /*Cleanup*/
     222        xDelete<IssmDouble>(beta0);
     223        xDelete<IssmDouble>(beta1);
     224        xDelete<IssmDouble>(phi);
     225        xDelete<IssmDouble>(noiseterms);
     226        xDelete<IssmDouble>(covmat);
    146227}/*}}}*/
    147228void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/

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

    r26271 r26477  
    1313void SmbGradientsx(FemModel* femmodel);
    1414void SmbGradientsElax(FemModel* femmodel);
     15void SmbautoregressionInitx(FemModel* femmodel);
     16void Smbautoregressionx(FemModel* femmodel);
    1517void Delta18oParameterizationx(FemModel* femmodel);
    1618void MungsmtpParameterizationx(FemModel* femmodel);

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

    r26457 r26477  
    409409        SmbAccurefEnum,
    410410        SmbAdThreshEnum,
     411        SmbAutoregressionInitialTimeEnum,
     412   SmbAutoregressionTimestepEnum,
     413   SmbAutoregressiveOrderEnum,
    411414        SmbAveragingEnum,
     415        SmbBeta0Enum,
     416   SmbBeta1Enum,
     417   SmbCovmatEnum,
    412418        SmbDesfacEnum,
    413419        SmbDpermilEnum,
    … …  
    438444        SmbIsturbulentfluxEnum,
    439445        SmbKEnum,
     446        SmbNumBasinsEnum,
    440447        SmbNumRequestedOutputsEnum,
    441448        SmbPfacEnum,
     449        SmbPhiEnum,
    442450        SmbRdlEnum,
    443451        SmbRequestedOutputsEnum,
    … …  
    863871        SmbAdiffiniEnum,
    864872        SmbAiniEnum,
     873        SmbBasinsIdEnum,
    865874        SmbBMaxEnum,
    866875        SmbBMinEnum,
    … …  
    949958        SmbTmeanEnum,
    950959        SmbTzEnum,
     960        SmbValuesAutoregressionEnum,
    951961        SmbVEnum,
    952962        SmbVmeanEnum,
    … …  
    14141424        SamplingSolutionEnum,
    14151425        SIAApproximationEnum,
     1426        SMBautoregressionEnum,
    14161427        SMBcomponentsEnum,
    14171428        SMBd18opddEnum,

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

    r26445 r26477  
    417417                case SmbAccurefEnum : return "SmbAccuref";
    418418                case SmbAdThreshEnum : return "SmbAdThresh";
     419                case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
     420                case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
     421                case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
    419422                case SmbAveragingEnum : return "SmbAveraging";
     423                case SmbBeta0Enum : return "SmbBeta0";
     424                case SmbBeta1Enum : return "SmbBeta1";
     425                case SmbCovmatEnum : return "SmbCovmat";
    420426                case SmbDesfacEnum : return "SmbDesfac";
    421427                case SmbDpermilEnum : return "SmbDpermil";
    … …  
    446452                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
    447453                case SmbKEnum : return "SmbK";
     454                case SmbNumBasinsEnum : return "SmbNumBasins";
    448455                case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs";
    449456                case SmbPfacEnum : return "SmbPfac";
     457                case SmbPhiEnum : return "SmbPhi";
    450458                case SmbRdlEnum : return "SmbRdl";
    451459                case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
    … …  
    542550                case AdjointpEnum : return "Adjointp";
    543551                case AdjointxEnum : return "Adjointx";
     552                case AdjointxBaseEnum : return "AdjointxBase";
     553                case AdjointxShearEnum : return "AdjointxShear";
    544554                case AdjointyEnum : return "Adjointy";
     555                case AdjointyBaseEnum : return "AdjointyBase";
     556                case AdjointyShearEnum : return "AdjointyShear";
    545557                case AdjointzEnum : return "Adjointz";
    546558                case AirEnum : return "Air";
    … …  
    865877                case SmbAdiffiniEnum : return "SmbAdiffini";
    866878                case SmbAiniEnum : return "SmbAini";
     879                case SmbBasinsIdEnum : return "SmbBasinsId";
    867880                case SmbBMaxEnum : return "SmbBMax";
    868881                case SmbBMinEnum : return "SmbBMin";
    … …  
    950963                case SmbTmeanEnum : return "SmbTmean";
    951964                case SmbTzEnum : return "SmbTz";
     965                case SmbValuesAutoregressionEnum : return "SmbValuesAutoregression";
    952966                case SmbVEnum : return "SmbV";
    953967                case SmbVmeanEnum : return "SmbVmean";
    … …  
    14131427                case SamplingSolutionEnum : return "SamplingSolution";
    14141428                case SIAApproximationEnum : return "SIAApproximation";
     1429                case SMBautoregressionEnum : return "SMBautoregression";
    14151430                case SMBcomponentsEnum : return "SMBcomponents";
    14161431                case SMBd18opddEnum : return "SMBd18opdd";

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

    r26445 r26477  
    426426              else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
    427427              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
     428              else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
     429              else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
     430              else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
    428431              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
     432              else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum;
     433              else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum;
     434              else if (strcmp(name,"SmbCovmat")==0) return SmbCovmatEnum;
    429435              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    430436              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    … …  
    455461              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
    456462              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
     463              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    457464              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    458465              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
     466              else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
    459467              else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
    460468              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    … …  
    498506              else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
    499507              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    500               else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    501512              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    502513              else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
    … …  
    506517              else if (strcmp(name,"Time")==0) return TimeEnum;
    507518              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"TimesteppingCouplingTime")==0) return TimesteppingCouplingTimeEnum;
     519              else if (strcmp(name,"TimesteppingCouplingTime")==0) return TimesteppingCouplingTimeEnum;
    512520              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    513521              else if (strcmp(name,"TimesteppingInterpForcing")==0) return TimesteppingInterpForcingEnum;
    … …  
    554562              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    555563              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
     564              else if (strcmp(name,"AdjointxBase")==0) return AdjointxBaseEnum;
     565              else if (strcmp(name,"AdjointxShear")==0) return AdjointxShearEnum;
    556566              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
     567              else if (strcmp(name,"AdjointyBase")==0) return AdjointyBaseEnum;
     568              else if (strcmp(name,"AdjointyShear")==0) return AdjointyShearEnum;
    557569              else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
    558570              else if (strcmp(name,"Air")==0) return AirEnum;
    … …  
    617629              else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum;
    618630              else if (strcmp(name,"DamageF")==0) return DamageFEnum;
    619               else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
    620635              else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
    621636              else if (strcmp(name,"DeltaIceThickness")==0) return DeltaIceThicknessEnum;
    … …  
    629644              else if (strcmp(name,"Str")==0) return StrEnum;
    630645              else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
     646              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
    635647              else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
    636648              else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
    … …  
    740752              else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum;
    741753              else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
    742               else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
    743758              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
    744759              else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum;
    … …  
    752767              else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
    753768              else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
     769              else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
    758770              else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
    759771              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
    … …  
    863875              else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
    864876              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    865               else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    866881              else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
    867882              else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum;
    … …  
    875890              else if (strcmp(name,"SmbAccumulatedPrecipitation")==0) return SmbAccumulatedPrecipitationEnum;
    876891              else if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
     892              else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
    881893              else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum;
    882894              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    … …  
    886898              else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
    887899              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
     900              else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
    888901              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    889902              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    … …  
    971984              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
    972985              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
     986              else if (strcmp(name,"SmbValuesAutoregression")==0) return SmbValuesAutoregressionEnum;
    973987              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
    974988              else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
    … …  
    984998              else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum;
    985999              else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
    986               else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
    9871004              else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    9881005              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    … …  
    9981015              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    9991016              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     1017              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    10041018              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    10051019              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    … …  
    11071121              else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    11081122              else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
    1109               else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
    11101127              else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
    11111128              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    … …  
    11211138              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    11221139              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
     1140              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    11271141              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    11281142              else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
    … …  
    12301244              else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
    12311245              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    1232               else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    12331250              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    12341251              else if (strcmp(name,"Dense")==0) return DenseEnum;
    … …  
    12441261              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    12451262              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
     1263              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    12501264              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    12511265              else if (strcmp(name,"Element")==0) return ElementEnum;
    … …  
    13531367              else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
    13541368              else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    1355               else if (strcmp(name,"MINI")==0) return MINIEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"MINI")==0) return MINIEnum;
    13561373              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    13571374              else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
    … …  
    13671384              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    13681385              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
     1386              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    13731387              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    13741388              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    … …  
    14461460              else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
    14471461              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
     1462              else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum;
    14481463              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    14491464              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    … …  
    14751490              else if (strcmp(name,"SegInput")==0) return SegInputEnum;
    14761491              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    1477               else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    14781496              else if (strcmp(name,"Separate")==0) return SeparateEnum;
    14791497              else if (strcmp(name,"Seq")==0) return SeqEnum;
    … …  
    14901508              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    14911509              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
     1510              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    14961511              else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
    14971512              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r24754 r26477  
    673673
    674674} /*}}}*/
     675void CholeskyRealPositiveDefinite(IssmDouble* Lchol, IssmDouble* A, int ndim) { /*{{{*/
     676   /*CholeskyRealPositiveDefinite   computes lower triangular matrix of the Cholesky decomposition of A
     677   Follows the Cholesky–Banachiewicz algorithm
     678   Lchol should point to an IssmDouble* of same dimensions as A*/
     679   for(int ii{0};ii<(ndim*ndim);ii++){Lchol[ii]=0;}; //ensure zero-initialization
     680   for(int ii{0};ii<ndim;ii++){
     681                for(int jj{0};jj<=ii;jj++){
     682                        double sum{0.0};
     683                        for(int kk{0};kk<jj;kk++) {
     684                                sum += Lchol[ii*ndim+kk]*Lchol[jj*ndim+kk];
     685                        }
     686                        if(ii==jj){Lchol[ii*ndim+jj] = sqrt(A[ii*ndim+jj]-sum);}
     687                        else{Lchol[ii*ndim+jj] = 1/Lchol[jj*ndim+jj] * (A[ii*ndim+jj]-sum);}
     688                }
     689        }
     690} /*}}}*/

  • issm/trunk-jpl/src/c/shared/Matrix/matrix.h

    r26469 r26477  
    3131void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale);
    3232void cellecho(int numcells, int m, ...);
     33void CholeskyRealPositiveDefinite(IssmDouble* Lchol, IssmDouble* A, int ndim);
    3334#endif //ifndef _MATRIXUTILS_H_

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r26047 r26477  
    241241                case 11: return SMBgradientscomponentsEnum;
    242242                case 12: return SMBsemicEnum;   
     243                case 55: return SMBautoregressionEnum;
    243244                default: _error_("Marshalled SMB code \""<<enum_in<<"\" not supported yet");
    244245        }
Note: See TracChangeset for help on using the changeset viewer.