Changeset 26859


Ignore:
Timestamp:
02/08/22 08:59:38 (3 years ago)
Author:
vverjans
Message:

NEW: added user-defined timestep for stochasticity in stochastiforcings

Location:
issm/trunk-jpl/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r26836 r26859  
    545545      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
    546546      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
     547      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.stochastictimestep",StochasticForcingTimestepEnum));
    547548      iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
    548549      if(num_fields<1) _error_("no stochasticforcing fields found");

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

    r26837 r26859  
    2626   femmodel->parameters->FindParam(&covariance,&M,&N,StochasticForcingCovarianceEnum); _assert_(M==dimtot); _assert_(N==dimtot);
    2727
     28        /*Check if this is a timestep for new noiseterms computation*/
     29        bool isstepforstoch = false;
     30        IssmDouble time,dt,starttime,tstep_stoch;
     31   femmodel->parameters->FindParam(&time,TimeEnum);
     32   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     33   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     34   femmodel->parameters->FindParam(&tstep_stoch,StochasticForcingTimestepEnum);
     35               
     36        #ifndef _HAVE_AD_
     37   if((fmod(time,tstep_stoch)<fmod((time-dt),tstep_stoch)) || (time<=starttime+dt) || tstep_stoch==dt) isstepforstoch = true;
     38   #else
     39   _error_("not implemented yet");
     40   #endif
     41
    2842   /*Compute noise terms*/
    29    IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
    30    my_rank=IssmComm::GetRank();
    31    if(my_rank==0){
    32       int fixedseed;
    33       IssmDouble time,dt,starttime;
    34       femmodel->parameters->FindParam(&time,TimeEnum);
    35       femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    36       femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    37                 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
    38       if(randomflag) fixedseed=-1;
    39       else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
    40                 /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
    41       IssmDouble* temparray = NULL;
    42       multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
    43       for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
     43        IssmDouble* noiseterms = xNew<IssmDouble>(dimtot);
     44   if(isstepforstoch){
     45                my_rank=IssmComm::GetRank();
     46        if(my_rank==0){
     47           int fixedseed;
     48                        /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
     49           if(randomflag) fixedseed=-1;
     50           else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
     51                        /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
     52           IssmDouble* temparray = NULL;
     53           multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
     54           for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
     55                        xDelete<IssmDouble>(temparray);
     56        }
     57        ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     58                femmodel->parameters->SetParam(noiseterms,dimtot,StochasticForcingNoisetermsEnum);
     59        }
     60        else{
     61                IssmDouble* temparray = NULL;
     62                femmodel->parameters->FindParam(&temparray,&N,StochasticForcingNoisetermsEnum); _assert_(N==dimtot);
     63                for(int i=0;i<dimtot;i++) noiseterms[i] = temparray[i];
    4464                xDelete<IssmDouble>(temparray);
    45    }
    46    ISSM_MPI_Bcast(noiseterms,dimtot,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    47    
     65        }
     66
    4867        int i=0;
    4968   for(int j=0;j<numfields;j++){

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

    r26837 r26859  
    113113syn keyword cConstant CalvingMinthicknessEnum
    114114syn keyword cConstant CalvingTestSpeedfactorEnum
     115syn keyword cConstant CalvingUseParamEnum
     116syn keyword cConstant CalvingScaleThetaEnum
     117syn keyword cConstant CalvingAmpAlphaEnum
     118syn keyword cConstant CalvingMidpointEnum
     119syn keyword cConstant CalvingNonlinearLawEnum
    115120syn keyword cConstant ConfigurationTypeEnum
    116121syn keyword cConstant ConstantsGEnum
    … …  
    428433syn keyword cConstant StochasticForcingIsStochasticForcingEnum
    429434syn keyword cConstant StochasticForcingIsWaterPressureEnum
     435syn keyword cConstant StochasticForcingNoisetermsEnum
    430436syn keyword cConstant StochasticForcingNumFieldsEnum
    431437syn keyword cConstant StochasticForcingRandomflagEnum
     438syn keyword cConstant StochasticForcingTimestepEnum
    432439syn keyword cConstant SolidearthSettingsReltolEnum
    433440syn keyword cConstant SolidearthSettingsSelfAttractionEnum
    … …  
    12481255syn keyword cConstant CalvingLevermannEnum
    12491256syn keyword cConstant CalvingTestEnum
     1257syn keyword cConstant CalvingParameterizationEnum
    12501258syn keyword cConstant CalvingVonmisesEnum
    12511259syn keyword cConstant CfdragcoeffabsgradEnum

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

    r26850 r26859  
    427427        StochasticForcingIsStochasticForcingEnum,
    428428        StochasticForcingIsWaterPressureEnum,
     429        StochasticForcingNoisetermsEnum,
    429430        StochasticForcingNumFieldsEnum,
    430431        StochasticForcingRandomflagEnum,
     432        StochasticForcingTimestepEnum,
    431433        SolidearthSettingsReltolEnum,
    432434        SolidearthSettingsSelfAttractionEnum,

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

    r26837 r26859  
    115115                case CalvingMinthicknessEnum : return "CalvingMinthickness";
    116116                case CalvingTestSpeedfactorEnum : return "CalvingTestSpeedfactor";
     117                case CalvingUseParamEnum : return "CalvingUseParam";
     118                case CalvingScaleThetaEnum : return "CalvingScaleTheta";
     119                case CalvingAmpAlphaEnum : return "CalvingAmpAlpha";
     120                case CalvingMidpointEnum : return "CalvingMidpoint";
     121                case CalvingNonlinearLawEnum : return "CalvingNonlinearLaw";
    117122                case ConfigurationTypeEnum : return "ConfigurationType";
    118123                case ConstantsGEnum : return "ConstantsG";
    … …  
    430435                case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
    431436                case StochasticForcingIsWaterPressureEnum : return "StochasticForcingIsWaterPressure";
     437                case StochasticForcingNoisetermsEnum : return "StochasticForcingNoiseterms";
    432438                case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
    433439                case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
     440                case StochasticForcingTimestepEnum : return "StochasticForcingTimestep";
    434441                case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";
    435442                case SolidearthSettingsSelfAttractionEnum : return "SolidearthSettingsSelfAttraction";
    … …  
    12501257                case CalvingLevermannEnum : return "CalvingLevermann";
    12511258                case CalvingTestEnum : return "CalvingTest";
     1259                case CalvingParameterizationEnum : return "CalvingParameterization";
    12521260                case CalvingVonmisesEnum : return "CalvingVonmises";
    12531261                case CfdragcoeffabsgradEnum : return "Cfdragcoeffabsgrad";

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

    r26837 r26859  
    106106syn keyword juliaConstC CalvingMinthicknessEnum
    107107syn keyword juliaConstC CalvingTestSpeedfactorEnum
     108syn keyword juliaConstC CalvingUseParamEnum
     109syn keyword juliaConstC CalvingScaleThetaEnum
     110syn keyword juliaConstC CalvingAmpAlphaEnum
     111syn keyword juliaConstC CalvingMidpointEnum
     112syn keyword juliaConstC CalvingNonlinearLawEnum
    108113syn keyword juliaConstC ConfigurationTypeEnum
    109114syn keyword juliaConstC ConstantsGEnum
    … …  
    421426syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
    422427syn keyword juliaConstC StochasticForcingIsWaterPressureEnum
     428syn keyword juliaConstC StochasticForcingNoisetermsEnum
    423429syn keyword juliaConstC StochasticForcingNumFieldsEnum
    424430syn keyword juliaConstC StochasticForcingRandomflagEnum
     431syn keyword juliaConstC StochasticForcingTimestepEnum
    425432syn keyword juliaConstC SolidearthSettingsReltolEnum
    426433syn keyword juliaConstC SolidearthSettingsSelfAttractionEnum
    … …  
    12411248syn keyword juliaConstC CalvingLevermannEnum
    12421249syn keyword juliaConstC CalvingTestEnum
     1250syn keyword juliaConstC CalvingParameterizationEnum
    12431251syn keyword juliaConstC CalvingVonmisesEnum
    12441252syn keyword juliaConstC CfdragcoeffabsgradEnum

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

    r26838 r26859  
    115115              else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
    116116              else if (strcmp(name,"CalvingTestSpeedfactor")==0) return CalvingTestSpeedfactorEnum;
     117              else if (strcmp(name,"CalvingUseParam")==0) return CalvingUseParamEnum;
     118              else if (strcmp(name,"CalvingScaleTheta")==0) return CalvingScaleThetaEnum;
     119              else if (strcmp(name,"CalvingAmpAlpha")==0) return CalvingAmpAlphaEnum;
     120              else if (strcmp(name,"CalvingMidpoint")==0) return CalvingMidpointEnum;
     121              else if (strcmp(name,"CalvingNonlinearLaw")==0) return CalvingNonlinearLawEnum;
    117122              else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
    118123              else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
    … …  
    132137              else if (strcmp(name,"CumGmtslc")==0) return CumGmtslcEnum;
    133138              else if (strcmp(name,"CumGmslc")==0) return CumGmslcEnum;
    134               else if (strcmp(name,"DamageC1")==0) return DamageC1Enum;
     139         else stage=2;
     140   }
     141   if(stage==2){
     142              if (strcmp(name,"DamageC1")==0) return DamageC1Enum;
    135143              else if (strcmp(name,"DamageC2")==0) return DamageC2Enum;
    136144              else if (strcmp(name,"DamageC3")==0) return DamageC3Enum;
    137145              else if (strcmp(name,"DamageC4")==0) return DamageC4Enum;
    138146              else if (strcmp(name,"Damage")==0) return DamageEnum;
    139          else stage=2;
    140    }
    141    if(stage==2){
    142               if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
     147              else if (strcmp(name,"DamageEquivStress")==0) return DamageEquivStressEnum;
    143148              else if (strcmp(name,"DamageEvolutionNumRequestedOutputs")==0) return DamageEvolutionNumRequestedOutputsEnum;
    144149              else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
    … …  
    255260              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    256261              else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum;
    257               else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
    258266              else if (strcmp(name,"InversionDxmin")==0) return InversionDxminEnum;
    259267              else if (strcmp(name,"InversionGatol")==0) return InversionGatolEnum;
    260268              else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
    261269              else if (strcmp(name,"InversionGrtol")==0) return InversionGrtolEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"InversionGttol")==0) return InversionGttolEnum;
     270              else if (strcmp(name,"InversionGttol")==0) return InversionGttolEnum;
    266271              else if (strcmp(name,"InversionIncompleteAdjoint")==0) return InversionIncompleteAdjointEnum;
    267272              else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum;
    … …  
    378383              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    379384              else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
    380               else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
    381389              else if (strcmp(name,"SolidearthPartitionOcean")==0) return SolidearthPartitionOceanEnum;
    382390              else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum;
    383391              else if (strcmp(name,"SolidearthNpartOcean")==0) return SolidearthNpartOceanEnum;
    384392              else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
     393              else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
    389394              else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum;
    390395              else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum;
    … …  
    439444              else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
    440445              else if (strcmp(name,"StochasticForcingIsWaterPressure")==0) return StochasticForcingIsWaterPressureEnum;
     446              else if (strcmp(name,"StochasticForcingNoiseterms")==0) return StochasticForcingNoisetermsEnum;
    441447              else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
    442448              else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
     449              else if (strcmp(name,"StochasticForcingTimestep")==0) return StochasticForcingTimestepEnum;
    443450              else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    444451              else if (strcmp(name,"SolidearthSettingsSelfAttraction")==0) return SolidearthSettingsSelfAttractionEnum;
    … …  
    499506              else if (strcmp(name,"SmbLapseRatePos")==0) return SmbLapseRatePosEnum;
    500507              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    501               else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    502512              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    503513              else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
    … …  
    506516              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    507517              else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
     518              else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
    512519              else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum;
    513520              else if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum;
    … …  
    622629              else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
    623630              else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
    624               else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
    625635              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
    626636              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
    … …  
    629639              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    630640              else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
     641              else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
    635642              else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum;
    636643              else if (strcmp(name,"BasalforcingsMeltrateFactor")==0) return BasalforcingsMeltrateFactorEnum;
    … …  
    745752              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    746753              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    747               else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    748758              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    749759              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    … …  
    752762              else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
    753763              else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
     764              else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
    758765              else if (strcmp(name,"NGia")==0) return NGiaEnum;
    759766              else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum;
    … …  
    868875              else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
    869876              else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
    870               else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
    871881              else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
    872882              else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum;
    … …  
    875885              else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
    876886              else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
     887              else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
    881888              else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
    882889              else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
    … …  
    991998              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    992999              else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
    993               else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
    9941004              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    9951005              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    … …  
    9981008              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
    9991009              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
     1010              else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    10041011              else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
    10051012              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    … …  
    11141121              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    11151122              else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    1116               else if (strcmp(name,"Vy")==0) return VyEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Vy")==0) return VyEnum;
    11171127              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    11181128              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    … …  
    11211131              else if (strcmp(name,"Vz")==0) return VzEnum;
    11221132              else if (strcmp(name,"VzFS")==0) return VzFSEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"VzHO")==0) return VzHOEnum;
     1133              else if (strcmp(name,"VzHO")==0) return VzHOEnum;
    11271134              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    11281135              else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
    … …  
    12371244              else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    12381245              else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
    1239               else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
    12401250              else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
    12411251              else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum;
    … …  
    12441254              else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum;
    12451255              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     1256              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    12501257              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    12511258              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
    … …  
    12801287              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    12811288              else if (strcmp(name,"CalvingTest")==0) return CalvingTestEnum;
     1289              else if (strcmp(name,"CalvingParameterization")==0) return CalvingParameterizationEnum;
    12821290              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    12831291              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    … …  
    13591367              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    13601368              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    1361               else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    13621373              else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum;
    13631374              else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
    … …  
    13671378              else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
    13681379              else if (strcmp(name,"GroundedAreaScaled")==0) return GroundedAreaScaledEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
     1380              else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
    13731381              else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
    13741382              else if (strcmp(name,"Gset")==0) return GsetEnum;
    … …  
    14821490              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    14831491              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    1484               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    14851496              else if (strcmp(name,"None")==0) return NoneEnum;
    14861497              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
    … …  
    14901501              else if (strcmp(name,"OceantransportAnalysis")==0) return OceantransportAnalysisEnum;
    14911502              else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
     1503              else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
    14961504              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
    14971505              else if (strcmp(name,"Open")==0) return OpenEnum;
    … …  
    16051613              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    16061614              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    1607               else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    16081619              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    16091620              else if (strcmp(name,"Tria")==0) return TriaEnum;
    … …  
    16131624              else if (strcmp(name,"Vertex")==0) return VertexEnum;
    16141625              else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum;
    1615          else stage=14;
    1616    }
    1617    if(stage==14){
    1618               if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
     1626              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
    16191627              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    16201628              else if (strcmp(name,"Vertices")==0) return VerticesEnum;

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

    r26848 r26859  
    1111                default_id                              = NaN;
    1212                covariance                              = NaN;
     13                stochastictimestep   = 0;
    1314                randomflag                              = 1;
    1415        end
    … …  
    3536
    3637                        num_fields = numel(self.fields);
     38                        if(self.stochastictimestep==0)
     39                                md.stochasticforcing.stochastictimestep = md.timestepping.time_step; %by default: stochastictimestep set to ISSM time step
     40                        end
    3741
    3842                        %Check that covariance matrix is positive definite
    … …  
    116120                                indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with other ar timesteps
    117121                                dimensions(indSMBar) = md.smb.num_basins;
     122                                if(md.smb.ar_timestep<self.stochastictimestep)
     123                                        error('SMBautoregression cannot have a timestep shorter than stochastictimestep');
     124                                end
    118125                        end
    119126                        if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
    120127                                indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with other ar timesteps
    121128                                dimensions(indTFar) = md.frontalforcings.num_basins;
     129                                if(md.frontalforcings.ar_timestep<self.stochastictimestep)
     130                                        error('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep');
     131                                end
    122132                        end
    123133                        if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression'))
    124134                                indBDWar        = find(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression')); %index of BDWar, now check for consistency with other ar timesteps
    125135                                dimensions(indBDWar) = md.basalforcings.num_basins;
     136                                if(md.basalforcings.ar_timestep<self.stochastictimestep)
     137                                        error('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep');
     138                                end
    126139                        end
    127140                        size_tot = sum(dimensions);
    … …  
    154167                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
    155168                        md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
     169                        md = checkfield(md,'fieldname','stochasticforcing.stochastictimestep','NaN',1,'Inf',1,'>=',md.timestepping.time_step);
    156170                        md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
    157171                        if(checkdefaults) %need to check the defaults
    … …  
    167181                        fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
    168182                        fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
     183                        fielddisplay(self,'stochastictimestep','timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)');
    169184                        fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    170185                        disp('Available fields:');
    … …  
    182197                                return
    183198                        else
     199
     200                                if(self.stochastictimestep==0)
     201                                        disp('      stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step');
     202                                        self.stochastictimestep = md.timestepping.time_step; %by default: stochastictimestep set to ISSM time step
     203                                end
    184204
    185205                                %Retrieve dimensionality of each field
    … …  
    225245                                WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
    226246                                WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
     247                                WriteData(fid,prefix,'object',self,'fieldname','stochastictimestep','format','Double','scale',yts);
    227248                                WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
    228249                        end

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

    r26855 r26859  
    1818        self.default_id = np.nan
    1919        self.covariance = np.nan
     20        self.stochastictimestep = 0
    2021        self.randomflag = 1
    2122
    … …  
    3233        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)'))
    3334        s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)'))
     35        s += '{}\n'.format(fielddisplay(self, 'stochastictimestep', 'timestep at which new stochastic noise terms are generated (default: md.timestepping.time_step)'))
    3436        s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
    3537        s += 'Available fields:\n'
    … …  
    5860
    5961        num_fields = len(self.fields)
     62        if(self.stochastictimestep==0):
     63            md.stochasticforcing.stochastictimestep = md.timestepping.time_step #by default: stochastictimestep set to ISSM time step
     64            print('      stochasticforcing.stocahstictimestep not specified: set to md.timestepping.time_step')
    6065
    6166        # Check that covariance matrix is positive definite (this is done internally by linalg)
    … …  
    119124            indSMBar = self.fields.index('SMBautoregression')  # Index of SMBar, now check for consistency with other timesteps
    120125            dimensions[indSMBar] = md.smb.num_basins
     126            if(md.smb.ar_timestep<self.stochastictimestep):
     127                raise TypeError('SMBautoregression cannot have a timestep shorter than stochastictimestep')
    121128        if ('FrontalForcingsRignotAutoregression' in self.fields):
    122129            indTFar = self.fields.index('FrontalForcingsRignotAutoregression')  # Index of TFar, now check for consistency with other timesteps
    123130            dimensions[indTFar] = md.frontalforcings.num_basins
     131            if(md.frontalforcings.ar_timestep<self.stochastictimestep):
     132                raise TypeError('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep')
    124133        if ('BasalforcingsDeepwaterMeltingRateAutoregression' in self.fields):
    125134            indBDWar = self.fields.index('BasalforcingsDeepwaterMeltingRateAutoregression')  # Index of BDWar, now check for consistency with other timesteps
    126135            dimensions[indTFar] = md.basalforcings.num_basins
     136            if(md.basalforcings.ar_timestep<self.stochastictimestep):
     137                raise TypeError('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep')
    127138        size_tot = np.sum(dimensions)
    128139
    … …  
    143154        md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', self.supportedstochforcings())
    144155        md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot])  # global covariance matrix
     156        md = checkfield(md, 'fieldname', 'stochasticforcing.stochastictimestep', 'NaN', 1,'Inf', 1, '>=', md.timestepping.time_step)
    145157        md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
    146158        if (checkdefaults):
    … …  
    164176        else:
    165177            num_fields = len(self.fields)
     178            if(self.stochastictimestep==0):
     179                md.stochasticforcing.stochastictimestep = md.timestepping.time_step #by default: stochastictimestep set to ISSM time step
    166180            # Retrieve dimensionality of each field
    167181            dimensions = self.defaultdimension * np.ones((num_fields))
    … …  
    197211            WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer')
    198212            WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
     213            WriteData(fid, prefix, 'object', self, 'fieldname', 'stochastictimestep', 'format', 'Double', 'scale', yts)
    199214            WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
    200215    # }}}

  • issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m

    r26535 r26859  
    4343                ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],...
    4444                ['./' basename '-' ncdate '.nc'],...
     45                '/media/vincent/TOSH4TB/GeorgiaTech/DataSearch/BedMachine/BedMachineGreenland-2021-04-20.nc',...
    4546                };
    4647
    … …  
    8788if strcmp(string,'mask') | strcmp(string,'source'),
    8889        %Need nearest neighbor to avoid interpolation between 0 and 2
    89         output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
     90        %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
     91        output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV
    9092else
    91         output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
     93        %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
     94        output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV
    9295end
    9396
Note: See TracChangeset for help on using the changeset viewer.