Index: /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp (revision 26747) @@ -76,5 +76,8 @@ iomodel->FetchDataToInput(inputs,elements,"md.sampling.kappa",SamplingKappaEnum); iomodel->FetchDataToInput(inputs,elements,"md.sampling.beta",SamplingBetaEnum,0.); + iomodel->FetchDataToInput(inputs,elements,"md.sampling.tau",SamplingTauEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleNoiseEnum,0.); + if(iomodel->solution_enum==TransientSolutionEnum) iomodel->FetchDataToInput(inputs,elements,"md.sampling.phi",SamplingPhiEnum,0.); }/*}}}*/ @@ -84,8 +87,6 @@ char** requestedoutputs = NULL; - parameters->AddObject(iomodel->CopyConstantObject("md.sampling.tau",SamplingTauEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.sampling.alpha",SamplingAlphaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.sampling.robin",SamplingRobinEnum)); - if(solution_enum==TransientSolutionEnum) parameters->AddObject(iomodel->CopyConstantObject("md.sampling.phi",SamplingPhiEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.sampling.seed",SamplingSeedEnum)); @@ -273,8 +274,11 @@ element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* newsample = xNew(numnodes); + IssmDouble* tau = xNew(numnodes); + element->GetInputListOnNodes(&tau[0],SamplingTauEnum); /*Use the dof list to index into the solution vector: */ for(int i=0;i(newsample[i])) _error_("NaN found in solution vector"); @@ -287,4 +291,5 @@ /*Free ressources:*/ xDelete(newsample); + xDelete(tau); xDelete(doflist); }/*}}}*/ @@ -418,2 +423,48 @@ *pMff=Mff; }/*}}}*/ +void SamplingAnalysis::UpdateTransientSample(FemModel * femmodel){ + + for(int j=0;jelements->Size();j++){ + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + UpdateTransientSample(element); + } + +} +void SamplingAnalysis::UpdateTransientSample(Element * element){ + + /*Intermediaries */ + IssmDouble phi, sample, noise; + + /*Fetch number vertices for this element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize new sample*/ + IssmDouble* sample_new = xNew(numvertices); + + /*Retrieve all inputs and parameters*/ + Input* sample_input=element->GetInput(SampleOldEnum); _assert_(sample_input); + Input* noise_input=element->GetInput(SampleNoiseEnum); _assert_(noise_input); + Input* phi_input=element->GetInput(SamplingPhiEnum); _assert_(phi_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + + /*Get input values at gauss points*/ + sample_input->GetInputValue(&sample,gauss); + noise_input->GetInputValue(&noise,gauss); + phi_input->GetInputValue(&phi,gauss); + + /*Get new sample*/ + sample_new[iv] = phi*sample + noise; + + } + + element->AddInput(SampleEnum,sample_new,element->GetElementType()); + + /*Clean up and return*/ + xDelete(sample_new); + delete gauss; + +} Index: /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h (revision 26746) +++ /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h (revision 26747) @@ -40,4 +40,7 @@ void MassMatrix(Matrix** pMff,FemModel* femmodel); + void UpdateTransientSample(FemModel* femmodel); + void UpdateTransientSample(Element * element); + }; #endif Index: /issm/trunk-jpl/src/c/cores/sampling_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sampling_core.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/cores/sampling_core.cpp (revision 26747) @@ -5,4 +5,6 @@ #include "../classes/classes.h" #include "../solutionsequences/solutionsequences.h" +#include "../analyses/analyses.h" // new +#include "../modules/modules.h" void sampling_core(FemModel* femmodel){ @@ -30,6 +32,32 @@ /*Generate random sample*/ if(VerboseSolution()) _printf0_(" call computational core\n"); + SamplingAnalysis* analysis = new SamplingAnalysis(); femmodel->SetCurrentConfiguration(SamplingAnalysisEnum); + + if(solution_type==TransientSolutionEnum){ + InputDuplicatex(femmodel,SampleEnum,SampleOldEnum); + + int seed; + femmodel->parameters->FindParam(&seed,SamplingSeedEnum); + if(seed>=0){ + int step; + femmodel->parameters->FindParam(&step,StepEnum); + seed = seed + 13923272*step; // change default seed for transient simulations (by considering an arbitrary shift based on the step number) + femmodel->parameters->SetParam(seed,SamplingSeedEnum); + } + + } + solutionsequence_sampling(femmodel); + + if(solution_type==TransientSolutionEnum){ + + InputDuplicatex(femmodel,SampleEnum,SampleNoiseEnum); + + analysis->UpdateTransientSample(femmodel); + + } + + delete analysis; if(save_results){ @@ -46,3 +74,5 @@ /*profiler*/ femmodel->profiler->Stop(SAMPLINGCORE); + + } Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 26747) @@ -150,5 +150,5 @@ femmodel->parameters->FindParam(&issampling,TransientIssamplingEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); - femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); + femmodel->parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum); #if defined(_HAVE_OCEAN_ ) Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26746) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26747) @@ -352,10 +352,8 @@ syn keyword cConstant ModelnameEnum syn keyword cConstant SamplingAlphaEnum -syn keyword cConstant SamplingPhiEnum syn keyword cConstant SamplingNumRequestedOutputsEnum syn keyword cConstant SamplingRequestedOutputsEnum syn keyword cConstant SamplingRobinEnum syn keyword cConstant SamplingSeedEnum -syn keyword cConstant SamplingTauEnum syn keyword cConstant SaveResultsEnum syn keyword cConstant SolidearthPartitionIceEnum @@ -803,6 +801,10 @@ syn keyword cConstant RheologyBbarAbsGradientEnum syn keyword cConstant SampleEnum +syn keyword cConstant SampleOldEnum +syn keyword cConstant SampleNoiseEnum syn keyword cConstant SamplingBetaEnum syn keyword cConstant SamplingKappaEnum +syn keyword cConstant SamplingPhiEnum +syn keyword cConstant SamplingTauEnum syn keyword cConstant SealevelEnum syn keyword cConstant SealevelGRDEnum @@ -1580,5 +1582,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1587,6 +1588,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1601,6 +1602,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1609,10 +1610,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1630,5 +1630,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1641,5 +1640,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1654,6 +1652,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1669,11 +1667,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1694,4 +1692,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26746) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26747) @@ -346,10 +346,8 @@ ModelnameEnum, SamplingAlphaEnum, - SamplingPhiEnum, SamplingNumRequestedOutputsEnum, SamplingRequestedOutputsEnum, SamplingRobinEnum, SamplingSeedEnum, - SamplingTauEnum, SaveResultsEnum, SolidearthPartitionIceEnum, @@ -534,5 +532,5 @@ TransientIsgroundinglineEnum, TransientIshydrologyEnum, - TransientIsmasstransportEnum, + TransientIsmasstransportEnum, TransientIsoceantransportEnum, TransientIsmovingfrontEnum, @@ -799,6 +797,10 @@ RheologyBbarAbsGradientEnum, SampleEnum, + SampleOldEnum, + SampleNoiseEnum, SamplingBetaEnum, SamplingKappaEnum, + SamplingPhiEnum, + SamplingTauEnum, SealevelEnum, SealevelGRDEnum, @@ -1036,5 +1038,5 @@ TemperatureSEMICEnum, ThermalforcingAutoregressionNoiseEnum, - ThermalforcingValuesAutoregressionEnum, + ThermalforcingValuesAutoregressionEnum, ThermalSpctemperatureEnum, ThicknessAbsGradientEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26747) @@ -354,10 +354,8 @@ case ModelnameEnum : return "Modelname"; case SamplingAlphaEnum : return "SamplingAlpha"; - case SamplingPhiEnum : return "SamplingPhi"; case SamplingNumRequestedOutputsEnum : return "SamplingNumRequestedOutputs"; case SamplingRequestedOutputsEnum : return "SamplingRequestedOutputs"; case SamplingRobinEnum : return "SamplingRobin"; case SamplingSeedEnum : return "SamplingSeed"; - case SamplingTauEnum : return "SamplingTau"; case SaveResultsEnum : return "SaveResults"; case SolidearthPartitionIceEnum : return "SolidearthPartitionIce"; @@ -805,6 +803,10 @@ case RheologyBbarAbsGradientEnum : return "RheologyBbarAbsGradient"; case SampleEnum : return "Sample"; + case SampleOldEnum : return "SampleOld"; + case SampleNoiseEnum : return "SampleNoise"; case SamplingBetaEnum : return "SamplingBeta"; case SamplingKappaEnum : return "SamplingKappa"; + case SamplingPhiEnum : return "SamplingPhi"; + case SamplingTauEnum : return "SamplingTau"; case SealevelEnum : return "Sealevel"; case SealevelGRDEnum : return "SealevelGRD"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26746) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26747) @@ -345,10 +345,8 @@ syn keyword juliaConstC ModelnameEnum syn keyword juliaConstC SamplingAlphaEnum -syn keyword juliaConstC SamplingPhiEnum syn keyword juliaConstC SamplingNumRequestedOutputsEnum syn keyword juliaConstC SamplingRequestedOutputsEnum syn keyword juliaConstC SamplingRobinEnum syn keyword juliaConstC SamplingSeedEnum -syn keyword juliaConstC SamplingTauEnum syn keyword juliaConstC SaveResultsEnum syn keyword juliaConstC SolidearthPartitionIceEnum @@ -796,6 +794,10 @@ syn keyword juliaConstC RheologyBbarAbsGradientEnum syn keyword juliaConstC SampleEnum +syn keyword juliaConstC SampleOldEnum +syn keyword juliaConstC SampleNoiseEnum syn keyword juliaConstC SamplingBetaEnum syn keyword juliaConstC SamplingKappaEnum +syn keyword juliaConstC SamplingPhiEnum +syn keyword juliaConstC SamplingTauEnum syn keyword juliaConstC SealevelEnum syn keyword juliaConstC SealevelGRDEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26747) @@ -360,10 +360,8 @@ else if (strcmp(name,"Modelname")==0) return ModelnameEnum; else if (strcmp(name,"SamplingAlpha")==0) return SamplingAlphaEnum; - else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum; else if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum; else if (strcmp(name,"SamplingRequestedOutputs")==0) return SamplingRequestedOutputsEnum; else if (strcmp(name,"SamplingRobin")==0) return SamplingRobinEnum; else if (strcmp(name,"SamplingSeed")==0) return SamplingSeedEnum; - else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum; @@ -383,10 +381,10 @@ else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum; else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum; + else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum; + else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum; - else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum; - else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; + if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; @@ -506,10 +504,10 @@ else if (strcmp(name,"Step")==0) return StepEnum; else if (strcmp(name,"Steps")==0) return StepsEnum; + else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; + else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; - else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum; - else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum; + if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum; else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum; else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum; @@ -629,10 +627,10 @@ else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum; + else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum; + else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum; - else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; - else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; + if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; @@ -752,10 +750,10 @@ else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum; else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum; + else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum; + else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum; - else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; - else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; + if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; @@ -823,6 +821,10 @@ else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sample")==0) return SampleEnum; + else if (strcmp(name,"SampleOld")==0) return SampleOldEnum; + else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum; else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum; else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum; + else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum; + else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum; else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; - else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum; - else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; + if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum; + else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum; + else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum; else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum; else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum; - else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; - else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; + if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; + else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; + else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; - else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; - else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else stage=10; } if(stage==10){ - if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; + if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; + else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; + else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; - else if (strcmp(name,"Boundary")==0) return BoundaryEnum; - else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; + if (strcmp(name,"Boundary")==0) return BoundaryEnum; + else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; + else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; else if (strcmp(name,"Incremental")==0) return IncrementalEnum; - else if (strcmp(name,"Indexed")==0) return IndexedEnum; - else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"ElementInput")==0) return ElementInputEnum; + if (strcmp(name,"Indexed")==0) return IndexedEnum; + else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; + else if (strcmp(name,"ElementInput")==0) return ElementInputEnum; else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; @@ -1488,10 +1490,10 @@ else if (strcmp(name,"Regular")==0) return RegularEnum; else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; - else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; - else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum; + if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum; + else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum; else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26746) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26747) @@ -51,6 +51,5 @@ Vector* ug = NULL; Vector* uf = NULL; - Vector* old_uf = NULL; // previous solution vector (for transient) - Vector* pf = NULL; + Vector* pf = NULL; Vector* df = NULL; Vector* ys=NULL; @@ -62,24 +61,9 @@ /*parameters:*/ - int solution_type, alpha, step, seed, nsize; - IssmDouble phi, tau; + int alpha, seed, nsize; /*Recover parameters: */ - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&seed,SamplingSeedEnum); femmodel->parameters->FindParam(&alpha,SamplingAlphaEnum); - femmodel->parameters->FindParam(&tau,SamplingTauEnum); - - /*Recover parameters for transient simulation: */ - if(solution_type==TransientSolutionEnum) - { - femmodel->parameters->FindParam(&step,StepEnum); - femmodel->parameters->FindParam(&phi,SamplingPhiEnum); - if(seed>=0) seed = seed + 13923272*step; // change default seed for transient simulations (by considering an arbitrary shif based on the step number) - - GetSolutionFromInputsx(&ug,femmodel); - Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters); - old_uf=uf->Duplicate(); - } /*CreateAnalysis*/ @@ -101,7 +85,4 @@ Reduceloadx(pf, Kfs, ys); delete Kfs; - - /*Copy old solution for transient run */ - if(solution_type==TransientSolutionEnum) uf->Copy(old_uf); /* Generate random RHS */ @@ -130,10 +111,4 @@ } - /* Divide results by tau */ - uf->Scale(1.0/tau); - - /* Update solution x_{t+1} = phi x_{t} + noise for transient */ - if(solution_type==TransientSolutionEnum) uf->AXPY(old_uf,phi); - /* Update input */ Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters); @@ -141,5 +116,5 @@ /*clean-up*/ - delete Kff; delete pf; delete df; delete uf; delete ys; delete old_uf; + delete Kff; delete pf; delete df; delete uf; delete ys; delete Ml; delete Mscale; delete analysis; Index: /issm/trunk-jpl/src/m/classes/sampling.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.m (revision 26746) +++ /issm/trunk-jpl/src/m/classes/sampling.m (revision 26747) @@ -9,5 +9,5 @@ tau = 0; beta = NaN; - phi = 0; + phi = NaN; alpha = 0; robin = 0; @@ -30,5 +30,5 @@ disp(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):')); fielddisplay(self,'kappa','coefficient of the identity operator'); - fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)'); + fielddisplay(self,'tau','scaling coefficient of the solution'); fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)'); @@ -47,12 +47,6 @@ function self = setdefaultparameters(self) % {{{ - %Scaling coefficient - self.tau=1; - %Apply Robin boundary conditions self.robin=0; - - %Temporal correlation factor - self.phi=0; %Exponent in fraction SPDE (default: 2, biLaplacian covariance @@ -77,10 +71,9 @@ md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]); if(md.sampling.robin) md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - end - md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); + end md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); @@ -91,7 +84,7 @@ WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','tau','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','phi','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean'); @@ -110,8 +103,8 @@ nu = self.alpha-1; - KAPPA = sqrt(8*nu)/lc; - TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2)); - md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1); - md.sampling.tau = TAU; + KAPPA = sqrt(8*nu)./lc; + TAU = sqrt(gamma(nu)./(gamma(self.alpha)*(4*pi)*KAPPA.^(2*nu).*sigma.^2)); + md.sampling.kappa = KAPPA.*ones(md.mesh.numberofvertices,1); + md.sampling.tau = TAU.*ones(md.mesh.numberofvertices,1); end % }}} @@ -121,5 +114,5 @@ writejsdouble(fid,[modelname '.sampling.tau'],self.tau); writejsdouble(fid,[modelname '.sampling.beta'],self.beta); - writejsdouble(fid,[modelname '.sampling.phi'],self.beta); + writejsdouble(fid,[modelname '.sampling.phi'],self.phi); writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha); writejsdouble(fid,[modelname '.sampling.robin'],self.robin); Index: /issm/trunk-jpl/test/NightlyRun/test134.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test134.m (revision 26746) +++ /issm/trunk-jpl/test/NightlyRun/test134.m (revision 26747) @@ -5,5 +5,7 @@ md = md.sampling.setparameters(md,2e5,1); md.sampling.seed = 100; +md.sampling.phi = zeros(md.mesh.numberofvertices,1); md.cluster=generic('name',oshostname(),'np',1); +md.cluster.np=1; md=solve(md,'smp');