Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 26003) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 26004) @@ -430,4 +430,30 @@ AC_MSG_RESULT($HAVE_MELTING) dnl }}} +dnl with-Sampling{{{ + +AC_ARG_WITH([Sampling], + + AS_HELP_STRING([--with-Sampling = YES], [compile with Sampling capabilities (default is yes)]), + + [SAMPLING=$withval],[SAMPLING=yes]) + +AC_MSG_CHECKING(for Sampling capability compilation) + + +HAVE_SAMPLING=no + +if test "x$SAMPLING" = "xyes"; then + + HAVE_SAMPLING=yes + + AC_DEFINE([_HAVE_SAMPLING_],[1],[with Sampling capability]) + +fi + +AM_CONDITIONAL([SAMPLING], [test x$HAVE_SAMPLING = xyes]) + +AC_MSG_RESULT($HAVE_SAMPLING) + +dnl }}} dnl with-Sealevelrise{{{ AC_ARG_WITH([Sealevelrise], Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 26003) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 26004) @@ -256,4 +256,5 @@ ./solutionsequences/solutionsequence_fct.cpp \ ./solutionsequences/solutionsequence_schurcg.cpp \ + ./solutionsequences/solutionsequence_sampling.cpp \ ./solutionsequences/convergence.cpp \ ./classes/Options/Options.cpp \ @@ -548,4 +549,10 @@ endif #}}} +#Sampling sources {{{ +if SAMPLING +issm_sources += \ + ./cores/sampling_core.cpp \ + ./analyses/SamplingAnalysis.cpp +endif #Slr sources {{{ if SEALEVELRISE Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 26004) @@ -104,4 +104,7 @@ case MeltingAnalysisEnum : return new MeltingAnalysis(); #endif + #ifdef _HAVE_SAMPLING_ + case SamplingAnalysisEnum : return new SamplingAnalysis(); + #endif #ifdef _HAVE_SEALEVELRISE_ case SealevelchangeAnalysisEnum : return new SealevelchangeAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp (revision 26004) @@ -76,5 +76,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.sampling.kappa",SamplingKappaEnum); iomodel->FetchDataToInput(inputs,elements,"md.sampling.beta",SamplingBetaEnum,0.); - iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.); + //iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26004) @@ -832,4 +832,8 @@ case EsaSolutionEnum: analyses_temp[numanalyses++]=EsaAnalysisEnum; + break; + + case SamplingSolutionEnum: + analyses_temp[numanalyses++]=SamplingAnalysisEnum; break; @@ -4828,11 +4832,11 @@ IssmDouble* bslcice_partition_serial=NULL; IssmDouble* partitionice=NULL; - int npartice,nel; + int npartice,nel; Vector* bslchydro_partition=NULL; IssmDouble* bslchydro_partition_serial=NULL; IssmDouble* partitionhydro=NULL; - int nparthydro; - + int nparthydro; + /*Initialize temporary vector that will be used to sum barystatic components @@ -5105,5 +5109,5 @@ pUp->SetValues(gsize,indices,Up,ADD_VAL); pUp->Assemble(); - + /*Add RSL to Up to find the geoid: */ Index: /issm/trunk-jpl/src/c/classes/Profiler.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Profiler.h (revision 26003) +++ /issm/trunk-jpl/src/c/classes/Profiler.h (revision 26004) @@ -27,13 +27,14 @@ #define ESACORE 14 /*Profiling ESA */ #define SLRCORE 15 /*Profiling SLR */ -#define MPISERIAL 16 /*Profiling MPISerial */ -#define SEDLOOP 17 /*Profiling MPISerial */ -#define SEDMatrix 18 /*Profiling MPISerial */ -#define SEDUpdate 19 /*Profiling MPISerial */ -#define EPLLOOP 20 /*Profiling MPISerial */ -#define EPLMasking 21 /*Profiling MPISerial */ -#define EPLMatrices 22 /*Profiling MPISerial */ -#define EPLUpdate 23 /*Profiling MPISerial */ -#define MAXPROFSIZE 24 /*Used to initialize static arrays*/ +#define SAMPLINGCORE 16 /*Profiling SAMPLING */ +#define MPISERIAL 17 /*Profiling MPISerial */ +#define SEDLOOP 18 /*Profiling MPISerial */ +#define SEDMatrix 19 /*Profiling MPISerial */ +#define SEDUpdate 20 /*Profiling MPISerial */ +#define EPLLOOP 21 /*Profiling MPISerial */ +#define EPLMasking 22 /*Profiling MPISerial */ +#define EPLMatrices 23 /*Profiling MPISerial */ +#define EPLUpdate 24 /*Profiling MPISerial */ +#define MAXPROFSIZE 25 /*Used to initialize static arrays*/ Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 26004) @@ -1,5 +1,5 @@ /*!\file: CorePointerFromSolutionEnum.cpp * \brief: return type of analyses, number of analyses and core solution function. - */ + */ #ifdef HAVE_CONFIG_H @@ -82,4 +82,7 @@ #endif break; + case SamplingSolutionEnum: + solutioncore=&sampling_core; + break; default: Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 26003) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 26004) @@ -1,4 +1,4 @@ /* - * cores.h: + * cores.h: */ @@ -54,4 +54,5 @@ void bmb_core(FemModel* femmodel); void damage_core(FemModel* femmodel); +void sampling_core(FemModel* femmodel); /*sealevel change cores:*/ Index: /issm/trunk-jpl/src/c/cores/sampling_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sampling_core.cpp (revision 26004) +++ /issm/trunk-jpl/src/c/cores/sampling_core.cpp (revision 26004) @@ -0,0 +1,48 @@ +/*!\file: sampling_core.cpp + * \brief: core of the sampling solution + */ + +#include "../classes/classes.h" +#include "../solutionsequences/solutionsequences.h" + +void sampling_core(FemModel* femmodel){ + + /*Start profiler*/ + femmodel->profiler->Start(SAMPLINGCORE); + + /*parameters: */ + int numoutputs; + bool save_results; + int solution_type; + char** requested_outputs = NULL; + + /*activate configuration*/ + femmodel->SetCurrentConfiguration(SamplingAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&numoutputs,SamplingNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SamplingRequestedOutputsEnum); + + if(VerboseSolution()) _printf0_(" Generating random samples\n"); + + /*Generate random sample*/ + if(VerboseSolution()) _printf0_(" call computational core\n"); + femmodel->SetCurrentConfiguration(SamplingAnalysisEnum); + solutionsequence_sampling(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs = SampleEnum; + femmodel->RequestedOutputsx(&femmodel->results,&outputs,1); + } + + if(solution_type==SamplingSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} + + /*profiler*/ + femmodel->profiler->Stop(SAMPLINGCORE); +} Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26003) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26004) @@ -1307,4 +1307,6 @@ syn keyword cConstant RecoveryAnalysisEnum syn keyword cConstant RiftfrontEnum +syn keyword cConstant SamplingAnalysisEnum +syn keyword cConstant SamplingSolutionEnum syn keyword cConstant SIAApproximationEnum syn keyword cConstant SMBcomponentsEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26003) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26004) @@ -1306,4 +1306,6 @@ RecoveryAnalysisEnum, RiftfrontEnum, + SamplingAnalysisEnum, + SamplingSolutionEnum, SIAApproximationEnum, SMBcomponentsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26004) @@ -1309,4 +1309,6 @@ case RecoveryAnalysisEnum : return "RecoveryAnalysis"; case RiftfrontEnum : return "Riftfront"; + case SamplingAnalysisEnum : return "SamplingAnalysis"; + case SamplingSolutionEnum : return "SamplingSolution"; case SIAApproximationEnum : return "SIAApproximation"; case SMBcomponentsEnum : return "SMBcomponents"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26003) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26004) @@ -1339,4 +1339,6 @@ 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 if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum; else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; @@ -1365,10 +1367,10 @@ else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum; else if (strcmp(name,"SealevelchangeAnalysis")==0) return SealevelchangeAnalysisEnum; - else if (strcmp(name,"SealevelchangeSolution")==0) return SealevelchangeSolutionEnum; - else if (strcmp(name,"Seg")==0) return SegEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"SegInput")==0) return SegInputEnum; + if (strcmp(name,"SealevelchangeSolution")==0) return SealevelchangeSolutionEnum; + else if (strcmp(name,"Seg")==0) return SegEnum; + else if (strcmp(name,"SegInput")==0) return SegInputEnum; else if (strcmp(name,"Segment")==0) return SegmentEnum; else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26004) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp (revision 26004) @@ -0,0 +1,160 @@ +/*!\file: solutionsequence_sampling.cpp + * \brief: numerical core of linear solutions + */ + +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../analyses/analyses.h" +#include + +void GaussianVector(Vector* ppf,int seed){/*{{{*/ + + /*Intermediaries*/ + double rdnumber; + + /*Define random number generator*/ + std::normal_distribution distribution(0.0,1.0); // Define probability distribution as the standard normal distribution + if(seed<0){ + std::random_device rd; + seed = rd(); + } + else + { + int my_rank; + ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&my_rank); + seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number) + } + std::default_random_engine generator(seed); + + int *local_indices = NULL; + IssmDouble *local_vector = NULL; + ppf->GetLocalVector(&local_vector,&local_indices); + + int M; + ppf->GetLocalSize(&M); + for(int i=0;iSetValue(local_indices[i],rdnumber,INS_VAL); + } + ppf->Assemble(); + +}/*}}}*/ +void InsertSample(IssmDouble* u,Vector* usample,int nsize,int numsam){/*{{{*/ + + /*Intermediaries*/ + int rstart = numsam*nsize; + double value; + + /*Fill global vector*/ + for (int i=0;iGetValue(&value,i); + u[rstart+i] = value; + } + +}/*}}}*/ +void solutionsequence_sampling(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; // previous solution vector (for transient) + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys=NULL; + + SamplingAnalysis* analysis = NULL; + + Vector* Ml = NULL; // diagonal lumped mass matrix + Vector* Mscale = NULL; // square root of diagonal lumped factor matrix + + /*parameters:*/ + int solution_type, alpha, step, seed, nsize; + IssmDouble phi, tau; + + /*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); + 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*/ + + analysis = new SamplingAnalysis(); + + analysis->LumpedMassMatrix(&Ml,femmodel); //Create lumped mass matrix + + if(alpha%2!=0) analysis->LumpedKMatrix(&Mscale,femmodel); /* Compute square root of lump mass matrix (for alpha even) or stiffness matrix (for alpha odd) */ + else{ + Mscale=Ml->Duplicate(); + Ml->Copy(Mscale); + } + Mscale->Pow(0.5); + + femmodel->UpdateConstraintsx(); + SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes); + Reduceloadx(pf, Kfs, ys); + delete Kfs; + + /*Copy old solution for transient run */ + if(solution_type==TransientSolutionEnum) uf->Copy(old_uf); + + /* Generate random RHS */ + GaussianVector(pf,seed); + + /* Multiply random RHS by square root of mass matrix (for alpha even) or stiffness matrix (for alpha odd) */ + pf->PointwiseMult(pf,Mscale); + + /*Go solve SPDE */ + femmodel->profiler->Start(SOLVER); + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + femmodel->profiler->Stop(SOLVER); + + /* Iterate to compute a realization for alpha>2 */ + for(int i=3;i<=alpha;i+=2){ + + /*Create RHS */ + uf->Copy(pf); + pf->PointwiseMult(pf,Ml); + + /*Go solve SPDE */ + femmodel->profiler->Start(SOLVER); + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + femmodel->profiler->Stop(SOLVER); + + } + + /* 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); + InputUpdateFromSolutionx(femmodel,ug); + + /*clean-up*/ + delete Kff; delete pf; delete df; delete uf; delete ys; delete old_uf; + delete Ml; delete Mscale; + delete analysis; + delete ug; + +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 26003) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 26004) @@ -1,4 +1,4 @@ /* - * solutionsequences.h: + * solutionsequences.h: */ @@ -25,4 +25,5 @@ void solutionsequence_adjoint_linear(FemModel* femmodel); void solutionsequence_schurcg(FemModel* femmodel); +void solutionsequence_sampling(FemModel* femmodel); /*convergence*/ Index: /issm/trunk-jpl/src/m/solve/solve.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.m (revision 26003) +++ /issm/trunk-jpl/src/m/solve/solve.m (revision 26004) @@ -23,4 +23,5 @@ % - 'Esa' or 'esa' % - 'Sealevelchange' or 'slc' +% - 'Sampling' or 'smp' % % Extra options: @@ -75,4 +76,6 @@ elseif strcmpi(solutionstring,'slc') || strcmpi(solutionstring,'Sealevelchange') solutionstring = 'SealevelchangeSolution'; +elseif strcmpi(solutionstring,'smp') || strcmpi(solutionstring,'Sampling') + solutionstring = 'SamplingSolution'; else error(['solutionstring ' solutionstring ' not supported!']);