Changeset 27780


Ignore:
Timestamp:
06/05/23 07:35:05 (22 months ago)
Author:
vverjans
Message:

NEW allowing for time-varying stochasticforcing covariance matrix

Location:
issm/trunk-jpl
Files:
3 added
9 edited

Legend:

Unmodified
Added
Removed

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

    r27498 r27780  
    575575        parameters->AddObject(new BoolParam(StochasticForcingIsWaterPressureEnum,false));
    576576   if(isstochasticforcing){
    577       int num_fields,stochastic_dim;
     577      int num_fields,num_tcov,stochastic_dim;
    578578      char** fields;
    579579      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
    580580      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
    581581      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.stochastictimestep",StochasticForcingTimestepEnum));
     582      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_timescovariance",StochasticForcingNumTimesCovarianceEnum));
    582583      iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
    583584      if(num_fields<1) _error_("no stochasticforcing fields found");
    … …  
    594595      parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
    595596      xDelete<IssmDouble>(transparam);
     597      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.timecovariance");
     598      parameters->AddObject(new DoubleVecParam(StochasticForcingTimeCovarianceEnum,transparam,N));
     599      xDelete<IssmDouble>(transparam);
    596600      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
    597601      parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));

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

    r27462 r27780  
    1414   /*Retrieve parameters*/
    1515   bool randomflag;
    16    int M,N,numfields,my_rank;
    17    int* fields            = NULL;
    18    int* dimensions        = NULL;
    19    IssmDouble* covariance = NULL;
     16   int M,N,numfields,numtcov,my_rank;
     17   int* fields                = NULL;
     18   int* dimensions            = NULL;
     19   IssmDouble* timecovariance = NULL;
     20   IssmDouble* covariance     = NULL;
    2021   femmodel->parameters->FindParam(&randomflag,StochasticForcingRandomflagEnum);
    2122   femmodel->parameters->FindParam(&numfields,StochasticForcingNumFieldsEnum);
     23   femmodel->parameters->FindParam(&numtcov,StochasticForcingNumTimesCovarianceEnum);
    2224   femmodel->parameters->FindParam(&fields,&N,StochasticForcingFieldsEnum);    _assert_(N==numfields);
    2325   femmodel->parameters->FindParam(&dimensions,&N,StochasticForcingDimensionsEnum);    _assert_(N==numfields);
     26   femmodel->parameters->FindParam(&timecovariance,&N,StochasticForcingTimeCovarianceEnum);    _assert_(N==numtcov);
    2427   int dimtot=0;
    2528   for(int i=0;i<numfields;i++) dimtot = dimtot+dimensions[i];
    26    femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==dimtot); _assert_(N==dimtot);
     29   femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==numtcov); _assert_(N==dimtot*dimtot);
    2730
    2831        /*Check if this is a timestep for new noiseterms computation*/
    … …  
    4548
    4649   /*Compute noise terms*/
    47         IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
     50        IssmDouble* timestepcovariance = xNew<IssmDouble>(dimtot*dimtot);
     51        IssmDouble* noiseterms         = xNew<IssmDouble>(dimtot);
    4852   if(isstepforstoch){
     53                /*Find covariance to be applied at current time step*/
     54                int itime;
     55                if(numtcov>1){
     56                        for(int i=0;i<numtcov;i++){
     57                                if(time>=timecovariance[i]) itime=i;
     58                        }
     59                }
     60                else itime=0;
     61                for(int i=0;i<dimtot*dimtot;i++) timestepcovariance[i] = covariance[itime*dimtot*dimtot+i];
    4962                my_rank=IssmComm::GetRank();
    5063        if(my_rank==0){
    … …  
    5568                        /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
    5669           IssmDouble* temparray = NULL;
    57            multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
     70           multivariateNormal(&temparray,dimtot,0.0,timestepcovariance,fixedseed);
    5871           for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
    5972                        xDelete<IssmDouble>(temparray);
    … …  
    233246   xDelete<int>(dimensions);
    234247   xDelete<IssmDouble>(covariance);
     248   xDelete<IssmDouble>(timecovariance);
     249   xDelete<IssmDouble>(timestepcovariance);
    235250   xDelete<IssmDouble>(noiseterms);
    236251}/*}}}*/

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

    r27758 r27780  
    507507syn keyword cConstant StochasticForcingNoisetermsEnum
    508508syn keyword cConstant StochasticForcingNumFieldsEnum
     509syn keyword cConstant StochasticForcingNumTimesCovarianceEnum
    509510syn keyword cConstant StochasticForcingRandomflagEnum
     511syn keyword cConstant StochasticForcingTimeCovarianceEnum
    510512syn keyword cConstant StochasticForcingTimestepEnum
    511513syn keyword cConstant SolidearthSettingsReltolEnum
    … …  
    17751777syn keyword cType Cfsurfacesquaretransient
    17761778syn keyword cType Channel
     1779syn keyword cType classes
    17771780syn keyword cType Constraint
    17781781syn keyword cType Constraints
    … …  
    17821785syn keyword cType ControlParam
    17831786syn keyword cType Covertree
     1787syn keyword cType DatasetInput
    17841788syn keyword cType DataSetParam
    1785 syn keyword cType DatasetInput
    17861789syn keyword cType Definition
    17871790syn keyword cType DependentObject
    … …  
    17961799syn keyword cType ElementInput
    17971800syn keyword cType ElementMatrix
     1801syn keyword cType Elements
    17981802syn keyword cType ElementVector
    1799 syn keyword cType Elements
    18001803syn keyword cType ExponentialVariogram
    18011804syn keyword cType ExternalResult
    … …  
    18041807syn keyword cType Friction
    18051808syn keyword cType Gauss
     1809syn keyword cType GaussianVariogram
     1810syn keyword cType gaussobjects
    18061811syn keyword cType GaussPenta
    18071812syn keyword cType GaussSeg
    18081813syn keyword cType GaussTetra
    18091814syn keyword cType GaussTria
    1810 syn keyword cType GaussianVariogram
    18111815syn keyword cType GenericExternalResult
    18121816syn keyword cType GenericOption
    … …  
    18251829syn keyword cType IssmDirectApplicInterface
    18261830syn keyword cType IssmParallelDirectApplicInterface
     1831syn keyword cType krigingobjects
    18271832syn keyword cType Load
    18281833syn keyword cType Loads
    … …  
    18351840syn keyword cType Matice
    18361841syn keyword cType Matlitho
     1842syn keyword cType matrixobjects
    18371843syn keyword cType MatrixParam
    18381844syn keyword cType Misfit
    … …  
    18471853syn keyword cType Observations
    18481854syn keyword cType Option
     1855syn keyword cType Options
    18491856syn keyword cType OptionUtilities
    1850 syn keyword cType Options
    18511857syn keyword cType Param
    18521858syn keyword cType Parameters
    … …  
    18621868syn keyword cType Regionaloutput
    18631869syn keyword cType Results
     1870syn keyword cType Riftfront
    18641871syn keyword cType RiftStruct
    1865 syn keyword cType Riftfront
    18661872syn keyword cType SealevelGeometry
    18671873syn keyword cType Seg
    18681874syn keyword cType SegInput
     1875syn keyword cType Segment
    18691876syn keyword cType SegRef
    1870 syn keyword cType Segment
    18711877syn keyword cType SpcDynamic
    18721878syn keyword cType SpcStatic
    … …  
    18871893syn keyword cType Vertex
    18881894syn keyword cType Vertices
    1889 syn keyword cType classes
    1890 syn keyword cType gaussobjects
    1891 syn keyword cType krigingobjects
    1892 syn keyword cType matrixobjects
    18931895syn keyword cType AdjointBalancethickness2Analysis
    18941896syn keyword cType AdjointBalancethicknessAnalysis

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

    r27758 r27780  
    501501        StochasticForcingNoisetermsEnum,
    502502        StochasticForcingNumFieldsEnum,
     503        StochasticForcingNumTimesCovarianceEnum,
    503504        StochasticForcingRandomflagEnum,
     505        StochasticForcingTimeCovarianceEnum,
    504506        StochasticForcingTimestepEnum,
    505507        SolidearthSettingsReltolEnum,

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

    r27758 r27780  
    509509                case StochasticForcingNoisetermsEnum : return "StochasticForcingNoiseterms";
    510510                case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
     511                case StochasticForcingNumTimesCovarianceEnum : return "StochasticForcingNumTimesCovariance";
    511512                case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
     513                case StochasticForcingTimeCovarianceEnum : return "StochasticForcingTimeCovariance";
    512514                case StochasticForcingTimestepEnum : return "StochasticForcingTimestep";
    513515                case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";

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

    r27758 r27780  
    500500syn keyword juliaConstC StochasticForcingNoisetermsEnum
    501501syn keyword juliaConstC StochasticForcingNumFieldsEnum
     502syn keyword juliaConstC StochasticForcingNumTimesCovarianceEnum
    502503syn keyword juliaConstC StochasticForcingRandomflagEnum
     504syn keyword juliaConstC StochasticForcingTimeCovarianceEnum
    503505syn keyword juliaConstC StochasticForcingTimestepEnum
    504506syn keyword juliaConstC SolidearthSettingsReltolEnum

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

    r27758 r27780  
    521521              else if (strcmp(name,"StochasticForcingNoiseterms")==0) return StochasticForcingNoisetermsEnum;
    522522              else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
     523              else if (strcmp(name,"StochasticForcingNumTimesCovariance")==0) return StochasticForcingNumTimesCovarianceEnum;
    523524              else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
     525              else if (strcmp(name,"StochasticForcingTimeCovariance")==0) return StochasticForcingTimeCovarianceEnum;
    524526              else if (strcmp(name,"StochasticForcingTimestep")==0) return StochasticForcingTimestepEnum;
    525527              else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    … …  
    627629              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    628630              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
    629               else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    630               else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
     634              if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
     635              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
     636              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    635637              else if (strcmp(name,"Step")==0) return StepEnum;
    636638              else if (strcmp(name,"Steps")==0) return StepsEnum;
    … …  
    750752              else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
    751753              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    752               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    753               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
     757              if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     758              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
     759              else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
    758760              else if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum;
    759761              else if (strcmp(name,"BasalStress")==0) return BasalStressEnum;
    … …  
    873875              else if (strcmp(name,"UGia")==0) return UGiaEnum;
    874876              else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum;
    875               else if (strcmp(name,"Gradient")==0) return GradientEnum;
    876               else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum;
     880              if (strcmp(name,"Gradient")==0) return GradientEnum;
     881              else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum;
     882              else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum;
    881883              else if (strcmp(name,"HydraulicPotentialOld")==0) return HydraulicPotentialOldEnum;
    882884              else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
    … …  
    996998              else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
    997999              else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
    998               else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
    999               else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
     1003              if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
     1004              else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
     1005              else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
    10041006              else if (strcmp(name,"SealevelBarystaticBpArea")==0) return SealevelBarystaticBpAreaEnum;
    10051007              else if (strcmp(name,"SealevelBarystaticBpLoad")==0) return SealevelBarystaticBpLoadEnum;
    … …  
    11191121              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    11201122              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    1121               else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    1122               else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum;
     1126              if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
     1127              else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum;
     1128              else if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum;
    11271129              else if (strcmp(name,"SmbMassBalanceSemic")==0) return SmbMassBalanceSemicEnum;
    11281130              else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
    … …  
    12421244              else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
    12431245              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    1244               else if (strcmp(name,"Vel")==0) return VelEnum;
    1245               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"VxBase")==0) return VxBaseEnum;
     1249              if (strcmp(name,"Vel")==0) return VelEnum;
     1250              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     1251              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    12501252              else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum;
    12511253              else if (strcmp(name,"Vx")==0) return VxEnum;
    … …  
    13651367              else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum;
    13661368              else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum;
    1367               else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
    1368               else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
     1372              if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
     1373              else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
     1374              else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
    13731375              else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
    13741376              else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
    … …  
    14881490              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    14891491              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    1490               else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    1491               else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
     1495              if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
     1496              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
     1497              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    14961498              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
    14971499              else if (strcmp(name,"FemModel")==0) return FemModelEnum;
    … …  
    16111613              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    16121614              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    1613               else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    1614               else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
     1618              if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
     1619              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
     1620              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    16191621              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    16201622              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    … …  
    17341736              else if (strcmp(name,"Sset")==0) return SsetEnum;
    17351737              else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
    1736               else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    1737               else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    17381738         else stage=15;
    17391739   }
    17401740   if(stage==15){
    1741               if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
     1741              if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
     1742              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
     1743              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    17421744              else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
    17431745              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;

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

    r27484 r27780  
    1111                default_id                              = NaN;
    1212                covariance                              = NaN;
     13                timecovariance                  = NaN;
    1314                stochastictimestep   = 0;
    1415                randomflag                              = 1;
    … …  
    4041                        end
    4142
    42                         %Check that covariance matrix is positive definite
    43                         try
    44                                 chol(self.covariance);
    45                         catch
    46                                 error('md.stochasticforcing.covariance is not positive definite');
     43
     44                        if(numel(size(self.covariance)==3))
     45                                numtcovmat = numel(self.covariance(1,1,:)); %number of covariance matrices in time
     46                                lsCovmats = {};
     47                                for ii=[1:numtcovmat] %loop over 3rd dimension
     48                                        lsCovmats{ii} = self.covariance(:,:,ii);
     49                        %Check that covariance matrix is positive definite
     50                        try
     51                                chol(self.covariance(:,:,ii));
     52                        catch
     53                                error('an entry in md.stochasticforcing.covariance is not positive definite');
     54                                        end
     55                                end
     56                        elseif(numel(size(self.covariance)==2))
     57                                numtcovmat = 1; %number of covariance matrices in time
     58                                lsCovmats = {self.covariance};
     59                        %Check that covariance matrix is positive definite
     60                        try
     61                                chol(self.covariance);
     62                        catch
     63                                error('md.stochasticforcing.covariance is not positive definite');
     64                        end
    4765                        end
    4866
    … …  
    176194                                if(indPwarma~=-1)
    177195                                        if(md.basalforcings.arma_timestep~=md.hydrology.arma_timestep)
    178                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
    179                                                 if any(crossentries~=0)
    180                                                         error('BasalforcingsDeepwaterMeltingRatearma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     196                                                for(ii=[1:numel(lsCovmats)])
     197                                                        covm = lsCovmats{ii};
     198                                                        crossentries = reshape(covm(1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
     199                                                        if any(crossentries~=0)
     200                                                                error('BasalforcingsDeepwaterMeltingRatearma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     201                                                        end
    181202                                                end
    182203                                        end
    183204                                elseif(indSdarma~=-1)
    184205                                        if(md.frontalforcings.sd_arma_timestep~=md.basalforcings.arma_timestep)
    185                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
    186                                                 if any(crossentries~=0)
    187                                                         error('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     206                                                for(ii=[1:numel(lsCovmats)])
     207                                                        covm = lsCovmats{ii};
     208                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
     209                                                        if any(crossentries~=0)
     210                                                                error('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     211                                                        end
    188212                                                end
    189213                                        end
    190214                                elseif(indSMBarma~=-1)
    191215                                        if(md.smb.arma_timestep~=md.basalforcings.arma_timestep)
    192                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
    193                                                 if any(crossentries~=0)
    194                                                         error('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     216                                                for(ii=[1:numel(lsCovmats)])
     217                                                        covm = lsCovmats{ii};
     218                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
     219                                                        if any(crossentries~=0)
     220                                                                error('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     221                                                        end
    195222                                                end
    196223                                        end
    197224                                elseif(indTFarma~=-1)
    198225                                        if(md.frontalforcings.arma_timestep~=md.basalforcings.arma_timestep)
    199                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
    200                                                 if any(crossentries~=0)
    201                                                         error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     226                                                for(ii=[1:numel(lsCovmats)])
     227                                                        covm = lsCovmats{ii};
     228                                                        crossentries = reshape(covm(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
     229                                                        if any(crossentries~=0)
     230                                                                error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     231                                                        end
    202232                                                end
    203233                                        end
    … …  
    206236                                if(indSdarma~=-1)
    207237                                        if(md.frontalforcings.sd_arma_timestep~=md.hydrology.arma_timestep)
    208                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
    209                                                 if any(crossentries~=0)
    210                                                         error('FrontalForcingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     238                                                for(ii=[1:numel(lsCovmats)])
     239                                                        covm = lsCovmats{ii};
     240                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
     241                                                        if any(crossentries~=0)
     242                                                                error('FrontalForcingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     243                                                        end
    211244                                                end
    212245                                        end
    213246                                elseif(indSMBarma~=-1)
    214247                                        if(md.smb.arma_timestep~=md.hydrology.arma_timestep)
    215                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
    216                                                 if any(crossentries~=0)
    217                                                         error('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     248                                                for(ii=[1:numel(lsCovmats)])
     249                                                        covm = lsCovmats{ii};
     250                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
     251                                                        if any(crossentries~=0)
     252                                                                error('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     253                                                        end
    218254                                                end
    219255                                        end
    220256                                elseif(indTFarma~=-1)
    221257                                        if(md.frontalforcings.arma_timestep~=md.hydrology.arma_timestep)
    222                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
    223                                                 if any(crossentries~=0)
    224                                                         error('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     258                                                for(ii=[1:numel(lsCovmats)])
     259                                                        covm = lsCovmats{ii};
     260                                                        crossentries = reshape(covm(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indPwarma-1)):sum(dimensions(1:indPwarma))),1,[]);
     261                                                        if any(crossentries~=0)
     262                                                                error('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance');
     263                                                        end
    225264                                                end
    226265                                        end
    … …  
    229268                                if(indSMBarma~=-1)
    230269                                        if(md.smb.arma_timestep~=md.frontalforcings.sd_arma_timestep)
    231                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma))),1,[]);
    232                                                 if any(crossentries~=0)
    233                                                         error('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
     270                                                for(ii=[1:numel(lsCovmats)])
     271                                                        covm = lsCovmats{ii};
     272                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma))),1,[]);
     273                                                        if any(crossentries~=0)
     274                                                                error('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
     275                                                        end
    234276                                                end
    235277                                        end
    236278                                elseif(indTFarma~=-1)
    237279                                        if(md.frontalforcings.sd_arma_timestep~=md.frontalforcings.arma_timestep)
    238                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
    239                                                 if any(crossentries~=0)
    240                                                         error('FrontalForcingsRignotarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
     280                                                for(ii=[1:numel(lsCovmats)])
     281                                                        covm = lsCovmats{ii};
     282                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
     283                                                        if any(crossentries~=0)
     284                                                                error('FrontalForcingsRignotarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
     285                                                        end
    241286                                                end
    242287                                        end
    … …  
    245290                                if(indTFarma~=-1)
    246291                                        if(md.smb.arma_timestep~=md.frontalforcings.arma_timestep)
    247                                                 crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
    248                                                 if any(crossentries~=0)
    249                                                         error('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance');
     292                                                for(ii=[1:numel(lsCovmats)])
     293                                                        covm = lsCovmats{ii};
     294                                                        crossentries = reshape(covm(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
     295                                                        if any(crossentries~=0)
     296                                                                error('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance');
     297                                                        end
    250298                                                end
    251299                                        end
    … …  
    256304                        md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
    257305                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
    258                         md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
     306                        md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot,numtcovmat]); %global covariance matrix
    259307                        md = checkfield(md,'fieldname','stochasticforcing.stochastictimestep','NaN',1,'Inf',1,'>=',md.timestepping.time_step);
    260308                        md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
     309                        if(numtcovmat>1) %check the time steps at which each covariance matrix starts to be applied
     310                                md = checkfield(md,'fieldname','stochasticforcing.timecovariance','NaN',1,'Inf',1,'>=',md.timestepping.start_time,'<=',md.timestepping.final_time,'size',[1,numtcovmat]);
     311                        end
    261312                        if(checkdefaults) %need to check the defaults
    262313                                md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0);
    … …  
    270321                        fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
    271322                        fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
    272                         fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
     323                        fielddisplay(self,'covariance',{'covariance matrix for within- and between-fields covariance (units must be squared field units)','multiple matrices can be concatenated along 3rd dimension to apply different covariances in time'});
     324                        fielddisplay(self,'timecovariance','starting dates at which covariances apply (only applicabe if multiple covariance matrices are prescribed)');
    273325                        fielddisplay(self,'stochastictimestep','timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)');
    274326                        fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    … …  
    325377                                end
    326378
     379                        if(numel(size(self.covariance)==3))
     380                                [nrow,ncol,numtcovmat] = size(self.covariance);
     381                                lsCovmats = {};
     382                                for ii=[1:numtcovmat] %loop over 3rd dimension
     383                                        lsCovmats{ii} = self.covariance(:,:,ii);
     384                                end
     385                                        if(md.timestepping.interp_forcing==1)
     386                                           disp('WARNING: md.timestepping.interp_forcing is 1, but be aware that there is no interpolation between covariance matrices');
     387                                           disp('         the changes between covariance matrices occur at the time steps specified in md.stochasticforcing.timecovariance');
     388                                        end
     389                                elseif(numel(size(self.covariance)==2))
     390                                [nrow,ncol] = size(self.covariance);
     391                                numtcovmat = 1; %number of covariance matrices in time
     392                                lsCovmats = {self.covariance};
     393                        end
     394   
    327395                                %Scaling covariance matrix (scale column-by-column and row-by-row)
    328396                                scaledfields = {'BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBarma','SMBforcing'}; %list of fields that need scaling *1/yts
    329                                 tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
    330                                 for i=1:num_fields
    331                                         if any(strcmp(scaledfields,self.fields(i)))
    332                                                 inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
    333                                                 for row=inds %scale rows corresponding to scaled field
    334                                                         tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
    335                                                 end
    336                                                 for col=inds %scale columns corresponding to scaled field
    337                                                         tempcovariance(:,col) = 1./yts*tempcovariance(:,col);
    338                                                 end
     397                                tempcovariance2d = zeros(numtcovmat,sum(nrow*ncol)); %covariance matrices in 2d array
     398                                % Loop over covariance matrices %
     399                                for kk=[1:numtcovmat]
     400                                        kkcov = self.covariance(:,:,kk); %extract covariance at index kk
     401                                        % Loop over the fields %
     402                                        for i=1:num_fields
     403                                                if any(strcmp(scaledfields,self.fields(i)))
     404                                                        inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
     405                                                        for row=inds %scale rows corresponding to scaled field
     406                                                                kkcov(row,:) = 1./yts*kkcov(row,:);
     407                                                        end
     408                                                        for col=inds %scale columns corresponding to scaled field
     409                                                                kkcov(:,col) = 1./yts*kkcov(:,col);
     410                                                        end
     411                                                end
     412                                        end
     413                                        % Save scaled covariance %
     414                                        for rr=[1:nrow]
     415                                                ind0 = 1+(rr-1)*ncol;
     416                                                tempcovariance2d(kk,ind0:ind0+ncol-1) = kkcov(rr,:);
    339417                                        end
    340418                                end
    … …  
    342420                                if isnan(self.default_id)
    343421                                        self.default_id = zeros(md.mesh.numberofelements,1);
     422                                end
     423                                %Set dummy timecovariance vector if a single covariance matrix is used
     424                                if(numtcovmat==1)
     425                                        self.timecovariance = [md.timestepping.start_time];
    344426                                end
    345427
    … …  
    349431                                WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed
    350432                                WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
    351                                 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
     433                                WriteData(fid,prefix,'data',numtcovmat,'name','md.stochasticforcing.num_timescovariance','format','Integer');
     434                                WriteData(fid,prefix,'data',tempcovariance2d,'name','md.stochasticforcing.covariance','format','DoubleMat');
     435                                WriteData(fid,prefix,'object',self,'fieldname','timecovariance','format','DoubleMat','scale',yts);
    352436                                WriteData(fid,prefix,'object',self,'fieldname','stochastictimestep','format','Double','scale',yts);
    353437                                WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');

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

    r27531 r27780  
    1919        self.default_id = np.nan
    2020        self.covariance = np.nan
     21        self.timecovariance = np.nan
    2122        self.stochastictimestep = 0
    2223        self.randomflag = 1
    … …  
    3334        s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)'))
    3435        s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'))
    35         s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)'))
     36        s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units),multiple matrices can be concatenated along 3rd dimension to apply different covariances in time'))
     37        s += '{}\n'.format(fielddisplay(self, 'timecovariance', 'starting dates at which covariances apply (only applicabe if multiple covariance matrices are prescribed)'))
    3638        s += '{}\n'.format(fielddisplay(self, 'stochastictimestep', 'timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)'))
    3739        s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
    … …  
    6971            print('      stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step')
    7072
    71         # Check that covariance matrix is positive definite (this is done internally by linalg)
    72         try:
    73             np.linalg.cholesky(self.covariance)
    74         except:
    75             raise TypeError('md.stochasticforcing.covariance is not positive definite')
     73        if(len(np.shape(self.covariance))==3):
     74           numtcovmat = np.shape(self.covariance)[2] #number of covariance matrices in time
     75           lsCovmats = []
     76           for ii in range(numtcovmat):
     77               lsCovmats.append(self.covariance[:,:,ii])
     78               try:
     79                   np.linalg.cholesky(self.covariance[:,:,ii])
     80               except:
     81                   raise TypeError('an entry in md.stochasticforcing.covariance is not positive definite')
     82        elif(len(np.shape(self.covariance))==2):
     83            numtcovmat = 1
     84            lsCovmats = [self.covariance]
     85            # Check that covariance matrix is positive definite (this is done internally by linalg)
     86            try:
     87                np.linalg.cholesky(self.covariance)
     88            except:
     89                raise TypeError('md.stochasticforcing.covariance is not positive definite')
    7690
    7791        # Check that all fields agree with the corresponding md class and if any field needs the default params
    … …  
    171185        if indBDWarma != -1:
    172186            if indPwarma != -1:
    173                 covsum = self.covariance[np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
    174                 if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
    175                     raise IOError('BasalforcingsDeepwaterMeltingRatarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
     187                for ii in range(len(lsCovmats)):
     188                    covm = lsCovmats[ii]
     189                    covsum = covm[np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
     190                    if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
     191                        raise IOError('BasalforcingsDeepwaterMeltingRatarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
    176192            elif indSdarma != -1:
    177                 covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
    178                 if md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
    179                     raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
     193                for ii in range(len(lsCovmats)):
     194                    covm = lsCovmats[ii]
     195                    covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
     196                    if md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
     197                        raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
    180198            elif indSMBarma != -1:
    181                 covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
    182                 if md.smb.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
    183                     raise IOError('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
     199                for ii in range(len(lsCovmats)):
     200                    covm = lsCovmats[ii]
     201                    covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
     202                    if md.smb.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
     203                        raise IOError('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
    184204            elif indTFarma != -1:
    185                 covsum = self.covariance[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
    186                 if md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
    187                     raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
     205                for ii in range(len(lsCovmats)):
     206                    covm = lsCovmats[ii]
     207                    covsum = covm[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
     208                    if md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep and np.any(covsum != 0):
     209                        raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
    188210        elif indPwarma != -1:
    189211            if indSdarma != -1:
    190                 covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
    191                 if md.frontalforcings.sd_arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
    192                     raise IOError('FrontalForingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance')
     212                for ii in range(len(lsCovmats)):
     213                    covm = lsCovmats[ii]
     214                    covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
     215                    if md.frontalforcings.sd_arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
     216                        raise IOError('FrontalForingsSubglacialDischargearma and hydrologyarmapw have different arma_timestep and non-zero covariance')
    193217            elif indSMBarma != -1:
    194                 covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
    195                 if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
    196                     raise IOError('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
     218                for ii in range(len(lsCovmats)):
     219                    covm = lsCovmats[ii]
     220                    covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
     221                    if md.smb.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
     222                        raise IOError('SMBarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
    197223            elif indTFarma != -1:
    198                 covsum = self.covariance[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
    199                 if md.frontalforcings.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
    200                     raise IOError('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
     224                for ii in range(len(lsCovmats)):
     225                    covm = lsCovmats[ii]
     226                    covsum = covm[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indPwarma]).astype(int):np.sum(dimensions[0:indPwarma + 1]).astype(int)]
     227                    if md.frontalforcings.arma_timestep != md.hydrology.arma_timestep and np.any(covsum != 0):
     228                        raise IOError('FrontalForcingsRignotarma and hydrologyarmapw have different arma_timestep and non-zero covariance')
    201229        elif indSdarma != -1:
    202230            if indSMBarma != -1:
    203                 covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)]
    204                 if md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep and np.any(covsum != 0):
    205                     raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance')
     231                for ii in range(len(lsCovmats)):
     232                    covm = lsCovmats[ii]
     233                    covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)]
     234                    if md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep and np.any(covsum != 0):
     235                        raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance')
    206236            elif indTFarma != -1:
    207                 covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
    208                 if md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0):
    209                     raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
     237                for ii in range(len(lsCovmats)):
     238                    covm = lsCovmats[ii]
     239                    covsum = covm[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
     240                    if md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0):
     241                        raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
    210242        elif indSMBarma != -1:
    211243            if indTFarma != -1:
    212                 covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
    213                 if md.smb.arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0):
    214                     raise IOError('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
     244                for ii in range(len(lsCovmats)):
     245                    covm = lsCovmats[ii]
     246                    covsum = covm[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
     247                    if md.smb.arma_timestep != md.frontalforcings.arma_timestep and np.any(covsum != 0):
     248                        raise IOError('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
    215249
    216250        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
    217251        md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', self.supportedstochforcings())
    218         md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot])  # global covariance matrix
     252        md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot, numtcovmat])  # global covariance matrix
    219253        md = checkfield(md, 'fieldname', 'stochasticforcing.stochastictimestep', 'NaN', 1,'Inf', 1, '>=', md.timestepping.time_step)
    220254        md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
     255        if(numtcovmat>1): #check the time steps at which each covariance matrix starts to be applie
     256            md = checkfield(md, 'fieldname', 'stochasticforcing.timecovariance', 'NaN', 1, 'Inf', 1, '>=',md.timestepping.start_time,'<=',md.timestepping.final_time,'size',[1,numtcovmat])  # global covariance matrix
    221257        if checkdefaults:
    222258            md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
    … …  
    260296                elif field == 'FrictionWaterPressure' and ispwHydroarma:
    261297                    dimensions[ind] = md.hydrology.num_basins
     298           
     299            if(len(np.shape(self.covariance))==3):
     300                nrow,ncol,numtcovmat = np.shape(self.covariance)
     301                lsCovmats = []
     302                for ii in range(numtcovmat):
     303                    lsCovmats.append(self.covariance[:,:,ii])
     304                if(md.timestepping.interp_forcing==1):
     305                    print('WARNING: md.timestepping.interp_forcing is 1, but be aware that there is no interpolation between covariance matrices')
     306                    print('         the changes between covariance matrices occur at the time steps specified in md.stochasticforcing.timecovariance')
     307            elif(len(np.shape(self.covariance))==2):
     308                nrow,ncol = np.shape(self.covariance)
     309                numtcovmat = 1
     310                lsCovmats = [self.covariance]
    262311
    263312            # Scaling covariance matrix (scale column-by-column and row-by-row)
    264313            scaledfields = ['BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBarma', 'SMBforcing']  # list of fields that need scaling * 1 / yts
    265             tempcovariance = np.copy(self.covariance)
    266             for i in range(num_fields):
    267                 if self.fields[i] in scaledfields:
    268                     inds = range(int(np.sum(dimensions[0:i])), int(np.sum(dimensions[0:i + 1])))
    269                     for row in inds:  # scale rows corresponding to scaled field
    270                         tempcovariance[row, :] = 1 / yts * tempcovariance[row, :]
    271                     for col in inds:  # scale columns corresponding to scaled field
    272                         tempcovariance[:, col] = 1 / yts * tempcovariance[:, col]
     314            tempcovariance2d = np.zeros((numtcovmat,nrow*ncol))
     315            # Loop over covariance matrices #
     316            for kk in range(numtcovmat):
     317                kkcov = lsCovmats[kk]
     318                # Loop over the fields #
     319                for i in range(num_fields):
     320                    if self.fields[i] in scaledfields:
     321                        inds = range(int(np.sum(dimensions[0:i])), int(np.sum(dimensions[0:i + 1])))
     322                        for row in inds:  # scale rows corresponding to scaled field
     323                            kkcov[row, :] = 1 / yts * kkcov[row, :]
     324                        for col in inds:  # scale columns corresponding to scaled field
     325                            kkcov[:, col] = 1 / yts * kkcov[:, col]
     326                # Save scaled covariance #
     327                for rr in range(nrow):
     328                    ind0 = rr*ncol
     329                    tempcovariance2d[kk,ind0:ind0+ncol] = np.copy(kkcov[rr,:])
     330                   
    273331            # Set dummy default_id vector if defaults not used
    274332            if np.any(np.isnan(self.default_id)):
    275333                self.default_id = np.zeros(md.mesh.numberofelements)
     334            # Set dummy timecovariance vector if a single covariance matrix is used
     335            if(numtcovmat==1):
     336                self.timecovariance = np.array([md.timestepping.start_time])
    276337            # Reshape dimensions as column array for marshalling
    277338            dimensions = dimensions.reshape(1, len(dimensions))
    … …  
    282343            WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'data', self.default_id - 1, 'format', 'IntMat', 'mattype', 2)  #12Nov2021 make sure this is zero-indexed!
    283344            WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer')
    284             WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
     345            WriteData(fid, prefix, 'data', numtcovmat, 'name', 'md.stochasticforcing.num_timescovariance', 'format', 'Integer')
     346            WriteData(fid, prefix, 'data', tempcovariance2d, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
     347            WriteData(fid, prefix, 'object', self, 'fieldname', 'timecovariance', 'format', 'DoubleMat', 'scale', yts)
    285348            WriteData(fid, prefix, 'object', self, 'fieldname', 'stochastictimestep', 'format', 'Double', 'scale', yts)
    286349            WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
Note: See TracChangeset for help on using the changeset viewer.