Commit 224d8f40 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch '19-develop-new-test-problem' into 'feature'

Resolve "Develop new Test Problem"

Closes #19

See merge request mpp/mlmc!25
parents c02f0403 87e1a2b1
Pipeline #148285 failed with stages
in 17 minutes and 50 seconds
......@@ -53,7 +53,7 @@ mpitest-mlmc:
tags: [ docker ]
elliptic-exercises-mlmc:
elliptic-experiments-mlmc:
stage: experiments
variables:
GIT_STRATEGY: none
......@@ -74,6 +74,25 @@ elliptic-exercises-mlmc:
paths: [ Elliptic\ Experiments.html ]
elliptic-experiments-mlmc2:
stage: experiments
variables:
GIT_STRATEGY: none
RUN_EXPERIMENTS: ''
image: ${REGISTRY}/${IMAGE_NAME_MLMC}
timeout: '1h'
script:
- cd /mpp/notebooks
- ls -al
- jupyter nbconvert --ExecutePreprocessor.timeout=2400
--execute --to html Elliptic\ Experiments\ 2.ipynb
- cp Elliptic\ Experiments\ 2.html $CI_PROJECT_DIR/
dependencies: [ "build-mlmc" ]
tags: [ docker ]
artifacts:
paths: [ Elliptic\ Experiments\ 2.html ]
deploy-mlmc:
stage: deploy
variables:
......
......@@ -29,7 +29,6 @@ lambda = [0.15, 0.15]
smoothing = 1.0
# ----- Plotting -----
GeneratorPlotting = 0
PDESolverPlotting = 0
# ----- Verbose -----
......
# ----- Problem Settings -----
Experiment = MLMCExperiment
#Experiment = ConvergenceTest
#Problem = StochasticPollution1D
#Problem = StochasticPollutionCosHat1D
#Problem = StochasticPollution2D
#Problem = StochasticPollutionCosHat2D
Problem = StochasticPollutionMollifiedBar2D
Problem = StochasticGaussHat2D
Model = DGTransport
#Model = PGTransport
#degree = 0
#degree = 1
degree = 2
Functional = Mass
#Funcitonal = Energy
#Functional = Outflow
Overlap = dG1
flux_alpha = 1 #Upwind: 1, Central: 0
# Streamline diffusion parameter
degree = 2
plevel = 2
Functional = Mass
#Funcitonal = Energy
#Functional = Outflow
# ----- Time Series -----
useDGTimeIntegrator = 0;
......@@ -36,45 +34,30 @@ Keps = 1e-5
# ----- Multilevel Monte Carlo -----
maxLevel = 7
MCParallel = false
epsilon = 0.01
mcOnly = false
initLevels = [4, 5, 6]
initSampleAmount = [8, 4, 2]
uniformSampleAmount = 100
onlyFine = false
initLevels = [3, 4, 5]
initSampleAmount = [12, 6, 3]
# ----- Stochastic Field -----
StochasticField = LogNormal
mean = 1.0
Mean = 0.0
sigma = 1.0
norm_p = 2
lambda1 = 0.10
lambda2 = 0.10
smoothing = 1.9
evtol = 1e-10
# ----- Solver -----
LinearReduction = 1e-12
LinearEpsilon = 1e-10
LinearSteps = 3000
NewtonSteps = 1
NewtonLineSearchSteps = 0
lambda = [0.15, 0.15]
smoothing = 1.0
# ----- Plotting -----
GeneratorPlotting = 1
MCPlotting = 1
PDESolverPlotting = 1
# ----- Verbose -----
MCVerbose = 1
PDEVerbose = 1
MainVerbose = 1
MLMCVerbose = 1
ConfigVerbose = 0
LinearVerbose = -1
NewtonVerbose = 1
GeneratorVerbose = 1
TimeIntegratorVerbose = 1
# ----- Logging -----
TimeLevel = -1
MuteLevel = -1
DebugLevel = -1
\ No newline at end of file
MainVerbose = 1
MeshVerbose = 1
ConfigVerbose = 1
LinearVerbose = 0
NewtonVerbose = 0
AssembleVerbose = 0
PDESolverVerbose = 0
GeneratorVerbose = 0
......@@ -47,8 +47,7 @@ NewtonSteps = 1
NewtonLineSearchSteps = 0
# ----- Plotting -----
GeneratorPlotting = 1
MCPlotting = 1
PDESolverPlotting = 1
# ----- Verbose -----
MCVerbose = 1
......
#include "m++.hpp"
#include "Main.hpp"
#ifdef COMPLEX
#error undef COMPLEX in src/Compiler.h
#endif
int main(int argc, char **argv) {
Config::setSearchPath("../mlmc/");
Config::setConfigFileName("m++.conf");
Mpp::initialize(&argc, argv);
Config::setSearchPath("../mlmc/");
Config::setConfigFileName("m++.conf");
Mpp::initialize(&argc, argv);
if(!PowerOfTwo(PPM->Size(0))) {
Warning("Number of process is not a power of 2."
"This might lead to undefined behaviour.")
}
MainProgram mainProgram;
return mainProgram.Run();
MainProgram mainProgram;
return mainProgram.Run();
}
......@@ -4,6 +4,12 @@
#include "montecarlo/MultilevelMonteCarlo.hpp"
bool PowerOfTwo(int n) {
if(n==0) return false;
return (ceil(log2(n)) == floor(log2(n)));
}
class MainProgram {
private:
int verbose = 1;
......
add_library(GENERATORS STATIC
NormalDistribution.cpp
UniformDistribution.cpp
SampleGeneratorContainer.cpp
algorithms/CirculantEmbedding.cpp
algorithms/HybridFluxGenerator.cpp
algorithms/NormalDistribution.cpp
algorithms/UniformDistribution.cpp
algorithms/SymmetricCovariance.cpp
algorithms/HybridFluxGenerator.cpp
)
target_link_libraries(GENERATORS BASICS PDESOLVER sprng MPP_LIBRARIES)
......@@ -4,27 +4,27 @@
#define PI 3.14159265358979323846
// Todo Ziggurat algorithm
void NormalDistributionReal::drawSample(const SampleID &id) {
// Todo Alternative: Ziggurat algorithm
if (internalCounter) {
updateCounter();
return;
}
double randNum1 = randNumberGen->sprng();
double randNum2 = randNumberGen->sprng();
firstSample = sqrt(-2.0 * log(randNum1)) * cos(2.0 * PI * randNum2);
secondSample = sqrt(-2.0 * log(randNum1)) * sin(2.0 * PI * randNum2);
if (internalCounter) {
updateCounter();
return;
}
double randNum1 = uniformDist.DrawAndEvalSample(id);
double randNum2 = uniformDist.DrawAndEvalSample(id);
firstSample = sqrt(-2.0 * log(randNum1)) * cos(2.0 * PI * randNum2);
secondSample = sqrt(-2.0 * log(randNum1)) * sin(2.0 * PI * randNum2);
updateCounter();
}
void NormalDistributionReal::updateCounter() {
internalCounter = internalCounter + 1;
internalCounter = internalCounter % 2;
internalCounter = internalCounter + 1;
internalCounter = internalCounter % 2;
}
Scalar NormalDistributionReal::EvalSample() {
if (internalCounter)
return firstSample;
else
return secondSample;
if (internalCounter)
return firstSample;
else
return secondSample;
}
#ifndef NORMALDISTRIBUTION_HPP
#define NORMALDISTRIBUTION_HPP
#include "../SampleGenerator.hpp"
//#define USE_MPI
#include "sprng_cpp.h"
#define SEED 985456376
/*
* Todo
*/
//class CollocationUniformDist : public SampleGenerator<Scalar> {
// list {
//
//
//
// }
//
// drawSample(id)
// list(id)
//
//};
#include "SampleGenerator.hpp"
#include "UniformDistribution.hpp"
class NormalDistributionReal : public SampleGenerator<Scalar> {
private:
Sprng *randNumberGen;
UniformDistributionReal uniformDist;
Scalar firstSample;
......@@ -39,22 +17,11 @@ private:
void drawSample(const SampleID &id) override;
void updateCounter();;
void updateCounter();
public:
// seed + PPM->Proc(0)
NormalDistributionReal(Meshes &meshes,
int streamID = 0, // PPM->Proc(0)
int numStreams = 1, // PPM->Size(0)
int gType = 0,
int seed = SEED) :
SampleGenerator(meshes) {
this->randNumberGen = SelectType(gType);
this->randNumberGen->init_sprng(streamID, numStreams,
seed, SPRNG_DEFAULT);
}
~NormalDistributionReal() { randNumberGen->free_sprng(); };
NormalDistributionReal(Meshes &meshes) :
SampleGenerator(meshes), uniformDist(UniformDistributionReal(meshes, 0.0, 1.0)) {}
Scalar EvalSample() override;
......@@ -70,10 +37,8 @@ private:
Complex sample;
void drawSample(const SampleID &id) override {
normalDist.DrawSample(id);
Scalar real = normalDist.EvalSample();
normalDist.DrawSample(id);
Scalar imag = normalDist.EvalSample();
Scalar real = normalDist.DrawAndEvalSample(id);
Scalar imag = normalDist.DrawAndEvalSample(id);
sample = {real, imag};
}
......@@ -98,10 +63,8 @@ private:
RVector sample;
void drawSample(const SampleID &id) override {
for (int i = 0; i < sample.size(); i++) {
normalDist.DrawSample(id);
sample[i] = normalDist.EvalSample();
}
for (int i = 0; i < sample.size(); i++)
sample[i] = normalDist.DrawAndEvalSample(id);
}
public:
......@@ -134,10 +97,8 @@ private:
CVector sample;
void drawSample(const SampleID &id) override {
for (int i = 0; i < sample.size(); i++) {
normalDist.DrawSample(id);
sample[i] = normalDist.EvalSample();
}
for (int i = 0; i < sample.size(); i++)
sample[i] = normalDist.DrawAndEvalSample(id);
}
public:
......@@ -174,10 +135,8 @@ private:
RMatrix sample;
void drawSample(const SampleID &id) override {
for (int i = 0; i < sample.rows(); i++) {
normalDist.DrawSample(id);
sample.InsertRow(normalDist.EvalSample(), i);
}
for (int i = 0; i < sample.rows(); i++)
sample.InsertRow(normalDist.DrawAndEvalSample(id), i);
}
public:
......@@ -210,10 +169,8 @@ private:
CMatrix sample;
void drawSample(const SampleID &id) override {
for (int i = 0; i < sample.rows(); i++) {
normalDist.DrawSample(id);
sample.InsertRow(normalDist.EvalSample(), i);
}
for (int i = 0; i < sample.rows(); i++)
sample.InsertRow(normalDist.DrawAndEvalSample(id), i);
}
public:
......
......@@ -20,20 +20,18 @@ static double abs_compl(const Complex &complex) {
template<typename T>
class SampleGenerator {
protected:
int plotting = 0;
int verbose = 0;
virtual void drawSample(const SampleID &id) = 0;
int commSplit;
virtual void plotSample(const SampleID &id) {};
virtual void drawSample(const SampleID &id) = 0;
public:
Meshes &meshes;
Meshes &meshes; // todo make const ref or remove completely
explicit SampleGenerator(Meshes &meshes) : meshes(meshes) {
explicit SampleGenerator(Meshes &meshes)
: meshes(meshes), commSplit(meshes.CommSplit()) {
// Todo add prefix
config.get("GeneratorPlotting", plotting);
config.get("GeneratorVerbose", verbose);
}
......@@ -43,7 +41,6 @@ public:
mout.StartBlock(Name());
vout(1) << id.IdString() << endl;
drawSample(id);
if (plotting) plotSample(id);
mout.EndBlock(verbose == 0);
}
......@@ -55,29 +52,64 @@ public:
Exit("Not implemented")
}
virtual T EvalSample(const Point &) {
virtual T EvalSample(const Point &x) {
Exit("Not implemented")
}
virtual T EvalSample(int, const Point &) {
virtual T EvalSample(int n, const Point &x) {
Exit("Not implemented")
}
virtual T EvalSample(double, const Point &) {
virtual T EvalSample(double t, const Point &x) {
Exit("Not implemented")
}
virtual T EvalSample(const cell &) {
virtual T EvalSample(const cell &c) {
Exit("Not implemented")
}
virtual T EvalSample(int, const cell &) {
virtual T EvalSample(int n, const cell &c) {
Exit("Not implemented")
}
virtual T EvalSample(double, const cell &) {
virtual T EvalSample(double t, const cell &c) {
Exit("Not implemented")
}
T DrawAndEvalSample(const SampleID &id) {
DrawSample(id);
return EvalSample();
}
T DrawAndEvalSample(const SampleID &id, const Point &x) {
DrawSample(id);
return EvalSample(x);
}
T DrawAndEvalSample(const SampleID &id, int n, const Point &x) {
DrawSample(id);
return EvalSample(n, x);
}
T DrawAndEvalSample(const SampleID &id, double t, const Point &x) {
DrawSample(id);
return EvalSample(t, x);
}
T DrawAndEvalSample(const SampleID &id, const cell &c) {
DrawSample(id);
return EvalSample(c);
}
T DrawAndEvalSample(const SampleID &id, int n, const cell &c) {
DrawSample(id);
return EvalSample(n, c);
}
T DrawAndEvalSample(const SampleID &id, double t, const cell &c) {
DrawSample(id);
return EvalSample(t, c);
}
};
class ScalarDummy : public SampleGenerator<Scalar> {
......
#include "SampleGeneratorContainer.hpp"
void SampleGeneratorContainer::init(Meshes &meshes) {
scalarGenerator = new ScalarDummy(meshes);
complexGenerator = new ComplexDummy(meshes);
vectorFieldGenerator = new VectorFieldDummy(meshes);
tensorGenerator = new TensorDummy(meshes);
rVectorGenerator = new RVectorDummy(meshes);
cVectorGenerator = new CVectorDummy(meshes);
rMatrixGenerator = new RMatrixDummy(meshes);
cMatrixGenerator = new CMatrixDummy(meshes);
}
void SampleGeneratorContainer::initWithGenNames(Meshes &meshes, GeneratorNames names) {
void SampleGeneratorContainer::init(Meshes &meshes, GeneratorNames names) {
for (auto const &genName: names) {
if (genName == "NormalDistributionReal")
scalarGenerator = new NormalDistributionReal(meshes);
......@@ -27,6 +16,12 @@ void SampleGeneratorContainer::initWithGenNames(Meshes &meshes, GeneratorNames n
if (genName == "NormalDistributionCMatrix")
cMatrixGenerator = new NormalDistributionCMatrix(meshes);
if (genName == "UniformDistributionReal")
scalarGenerator = new UniformDistributionReal(meshes, 0.0, 1.0);
if (genName == "UniformDistributionRVector")
rVectorGenerator = new UniformDistributionRVector(meshes, 2, -1.0, 1.0);
if (genName == "CirculantEmbedding") {
if (meshes.dim() == 1)
tensorGenerator = new CirculantEmbedding1D(meshes);
......@@ -39,12 +34,20 @@ void SampleGeneratorContainer::initWithGenNames(Meshes &meshes, GeneratorNames n
}
void SampleGeneratorContainer::DrawSample(const SampleID &id) {
scalarGenerator->DrawSample(id);
complexGenerator->DrawSample(id);
vectorFieldGenerator->DrawSample(id);
tensorGenerator->DrawSample(id);
rVectorGenerator->DrawSample(id);
cVectorGenerator->DrawSample(id);
rMatrixGenerator->DrawSample(id);
cMatrixGenerator->DrawSample(id);
if (scalarGenerator != nullptr)
scalarGenerator->DrawSample(id);
if (complexGenerator != nullptr)
complexGenerator->DrawSample(id);
if (vectorFieldGenerator != nullptr)
vectorFieldGenerator->DrawSample(id);
if (tensorGenerator != nullptr)
tensorGenerator->DrawSample(id);
if (rVectorGenerator != nullptr)
rVectorGenerator->DrawSample(id);
if (cVectorGenerator != nullptr)
cVectorGenerator->DrawSample(id);
if (rMatrixGenerator != nullptr)
rMatrixGenerator->DrawSample(id);
if (cMatrixGenerator != nullptr)
cMatrixGenerator->DrawSample(id);
}
......@@ -4,42 +4,57 @@
#include "SampleGenerator.hpp"
#include "algorithms/CirculantEmbedding.hpp"
#include "algorithms/HybridFluxGenerator.hpp"
#include "algorithms/NormalDistribution.hpp"
#include "algorithms/UniformDistribution.hpp"
#include "NormalDistribution.hpp"
#include "UniformDistribution.hpp"
typedef std::vector<std::string> GeneratorNames;
class SampleGeneratorContainer {
void init(Meshes &meshes);
void initWithGenNames(Meshes &meshes, GeneratorNames names);
void init(Meshes &meshes, GeneratorNames names);
public:
SampleGenerator<Scalar> *scalarGenerator;
SampleGenerator<Scalar> *scalarGenerator = nullptr;
SampleGenerator<Complex> *complexGenerator;
SampleGenerator<Complex> *complexGenerator = nullptr;
SampleGenerator<VectorField> *vectorFieldGenerator;
SampleGenerator<VectorField> *vectorFieldGenerator = nullptr;
SampleGenerator<Tensor> *tensorGenerator;
SampleGenerator<Tensor> *tensorGenerator = nullptr;
SampleGenerator<RVector> *rVectorGenerator;
SampleGenerator<RVector> *rVectorGenerator = nullptr;
SampleGenerator<CVector> *cVectorGenerator;
SampleGenerator<CVector> *cVectorGenerator = nullptr;
SampleGenerator<RMatrix> *rMatrixGenerator;
SampleGenerator<RMatrix> *rMatrixGenerator = nullptr;
SampleGenerator<CMatrix> *cMatrixGenerator;
SampleGenerator<CMatrix> *cMatrixGenerator = nullptr;
SampleGeneratorContainer(GeneratorNames names, Meshes &meshes) {
init(meshes);
initWithGenNames(meshes, names);
}
SampleGeneratorContainer(GeneratorNames names, Meshes &meshes) {
init(meshes, names);
}
// todo clean this up
~SampleGeneratorContainer() {
if (scalarGenerator != nullptr)
delete scalarGenerator;
if (complexGenerator != nullptr)
delete complexGenerator;
if (vectorFieldGenerator != nullptr)
delete vectorFieldGenerator;
if (tensorGenerator != nullptr)