Commit cb1622cd authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Template Based Sample Generator

parent 266a8c90
......@@ -20,7 +20,6 @@ protected:
public:
explicit IStochasticTransportAssemble(IStochasticTransportProblem *problem) :
problem(problem) {
InfoEntries problemEntries = problem->GetInfoEntries();
};
void DrawSample(const SampleID &id) {
......
#include "IStochasticProblem.hpp"
#include "generators/CirculantEmbedding.hpp"
#include "generators/HybridFluxGenerator.hpp"
SampleGenerator *IStochasticProblem::createGenerator(GeneratorName genName,
Meshes &meshes) {
if (genName == "DummyGenerator")
return new DummyGenerator(meshes);
// if (genName == "UniformRandomNumberGenerator")
// return new UniformRandomNumberGenerator(meshes);
// if (genName == "NormalRandomNumberGenerator")
// return new NormalRandomNumberGenerator(meshes);
if (genName == "CirculantEmbedding")
return new MultilevelCirculantEmbedding(meshes);
if (genName == "HybridFluxGenerator")
return new MultilevelHybridFluxGenerator(meshes);
Exit("Generator not implemented")
}
#include "generators/HybridFluxGenerator.hpp"
\ No newline at end of file
......@@ -3,52 +3,76 @@
#include "Algebra.hpp"
#include "basics/Sample.hpp"
#include "generators/SampleGenerator.hpp"
#include "basics/Utilities.hpp"
#include "generators/CirculantEmbedding.hpp"
//#include "generators/HybridFluxGenerator.hpp"
#include "generators/SampleGenerator.hpp"
typedef std::string GeneratorName;
typedef std::vector<GeneratorName> GeneratorNames;
typedef std::vector<std::string> GeneratorNames;
typedef std::map<GeneratorName, SampleGenerator *> SampleGenerators;
class SampleGeneratorContainer {
void init(Meshes &meshes) {
scalarGenerator = new ScalarDummy(meshes);
complexGenerator = new ComplexDummy(meshes);
vectorFieldGenerator = new VectorFieldDummy(meshes);
tensorGenerator = new TensorDummy(meshes);
scalarSequence1DGenerator = new ScalarSequence1DDummy(meshes);
complexSequence1DGenerator = new ComplexSequence1DDummy(meshes);
scalarSequence2DGenerator = new ScalarSequence2DGenerator(meshes);
complexSequence2DGenerator = new ComplexSequence2DGenerator(meshes);
}
public:
SampleGenerator<Scalar> *scalarGenerator;
SampleGenerator<Complex> *complexGenerator;
SampleGenerator<VectorField> *vectorFieldGenerator;
SampleGenerator<Tensor> *tensorGenerator;
SampleGenerator<ScalarSequence1D> *scalarSequence1DGenerator;
SampleGenerator<ComplexSequence1D> *complexSequence1DGenerator;
SampleGenerator<ScalarSequence2D> *scalarSequence2DGenerator;
SampleGenerator<ComplexSequence2D> *complexSequence2DGenerator;
SampleGeneratorContainer(GeneratorNames names, Meshes &meshes) {
init(meshes);
for (auto const &genName: names) {
if (genName == "CirculantEmbedding")
tensorGenerator = new MultilevelCirculantEmbedding(meshes);
// if (genName == "HybridFluxGenerator")
// new MultilevelHybridFluxGenerator(meshes);
}
}
void DrawSample(const SampleID &id) {
scalarGenerator->DrawSample(id);
complexGenerator->DrawSample(id);
vectorFieldGenerator->DrawSample(id);
tensorGenerator->DrawSample(id);
scalarSequence1DGenerator->DrawSample(id);
complexSequence1DGenerator->DrawSample(id);
scalarSequence2DGenerator->DrawSample(id);
complexSequence2DGenerator->DrawSample(id);
}
};
class IStochasticProblem {
private:
int verbose = 1;
InfoEntries entries;
GeneratorNames genNames;
SampleGenerator *createGenerator(GeneratorName genName, Meshes &meshes);
protected:
SampleGenerators generators;
SampleGeneratorContainer genContainer;
public:
IStochasticProblem(GeneratorNames genNames, Meshes &meshes) : genNames(genNames) {
IStochasticProblem(Meshes &meshes, GeneratorNames genNames) :
genNames(genNames), genContainer(genNames, meshes) {
config.get("ProblemVerbose", verbose);
for (auto &name : genNames) {
entries.push_back(InfoEntry("Generator", name));
generators[name] = createGenerator(name, meshes);
}
}
virtual ~IStochasticProblem() {
for (auto &gen: generators) {
delete gen.second;
}
generators.clear();
};
void DrawSample(const SampleID &id) {
for (auto &gen : generators)
gen.second->DrawSample(id);
}
InfoEntries GetInfoEntries() {
return entries;
genContainer.DrawSample(id);
}
virtual void PrintInfo() const {
......
......@@ -6,8 +6,8 @@
class IStochasticEllipticProblem : public IStochasticProblem {
public:
IStochasticEllipticProblem(GeneratorNames genNames, Meshes &meshes)
: IStochasticProblem(genNames, meshes) {}
IStochasticEllipticProblem(Meshes &meshes, GeneratorNames genNames = {})
: IStochasticProblem(meshes, genNames) {}
virtual Scalar Solution(const Point &x) const = 0;
......@@ -21,8 +21,7 @@ public:
class StochasticLaplace1D : public IStochasticEllipticProblem {
public:
StochasticLaplace1D(Meshes &meshes) :
IStochasticEllipticProblem(GeneratorNames{"CirculantEmbedding"},
meshes) {};
IStochasticEllipticProblem(meshes, GeneratorNames{"CirculantEmbedding"}) {};
Scalar Solution(const Point &x) const override {
return 1 - x[0];
......@@ -37,8 +36,7 @@ public:
}
Tensor Permeability(const cell &c) const override {
return this->generators.find("CirculantEmbedding")->
second->EvalTensorSample(c);
return this->genContainer.tensorGenerator->EvalSample(c);
}
string Name() const override { return "StochasticLaplace1D"; }
......@@ -46,9 +44,7 @@ public:
class Laplace1D : public IStochasticEllipticProblem {
public:
Laplace1D(Meshes &meshes) :
IStochasticEllipticProblem(GeneratorNames{"DummyGenerator"},
meshes) {};
Laplace1D(Meshes &meshes) : IStochasticEllipticProblem(meshes) {};
Scalar Solution(const Point &x) const override {
return 1 - x[0];
......@@ -72,8 +68,7 @@ public:
class StochasticLaplace2D : public IStochasticEllipticProblem {
public:
StochasticLaplace2D(Meshes &meshes) :
IStochasticEllipticProblem(GeneratorNames{"CirculantEmbedding"},
meshes) {}
IStochasticEllipticProblem(meshes, GeneratorNames{"CirculantEmbedding"}) {}
Scalar Solution(const Point &x) const override {
return -x[1];
......@@ -86,8 +81,7 @@ public:
Scalar Load(const Point &x) const override { return 0.0; }
Tensor Permeability(const cell &c) const override {
return this->generators.find("CirculantEmbedding")->
second->EvalTensorSample(c);
return this->genContainer.tensorGenerator->EvalSample(c);
}
string Name() const override { return "StochasticLaplace2D"; }
......@@ -95,9 +89,7 @@ public:
class Laplace2D : public IStochasticEllipticProblem {
public:
Laplace2D(Meshes &meshes) :
IStochasticEllipticProblem(GeneratorNames{"DummyGenerator"},
meshes) {}
Laplace2D(Meshes &meshes) : IStochasticEllipticProblem(meshes) {}
Scalar Solution(const Point &x) const override {
return -x[1];
......
......@@ -37,8 +37,8 @@ protected:
double reaction = 1.0;
public:
IStochasticReactionProblem(GeneratorNames genNames, Meshes &meshes) :
IStochasticProblem(genNames, meshes) {
IStochasticReactionProblem(Meshes &meshes, GeneratorNames genNames={}) :
IStochasticProblem(meshes, genNames) {
config.get("Diffusion", diffusion);
config.get("Convection", convection);
config.get("Reaction", reaction);
......@@ -111,7 +111,7 @@ class ExponentialReaction2D : public IStochasticReactionProblem {
public:
ExponentialReaction2D(Meshes &meshes) :
IStochasticReactionProblem(GeneratorNames{}, meshes) {};
IStochasticReactionProblem(meshes) {};
Scalar Concentration(double t, const Point &x) const override {
return u0(x);
......@@ -142,7 +142,7 @@ private:
double r1 = 0.0;
public:
LogisticReaction2D(Meshes &meshes) :
IStochasticReactionProblem(GeneratorNames{}, meshes) {
IStochasticReactionProblem(meshes) {
config.get("Reaction0", r0);
config.get("Reaction1", r1);
}
......
......@@ -9,8 +9,8 @@ private:
bool rhs = true;
public:
IStochasticTransportProblem(GeneratorNames genNames, Meshes &meshes)
: IStochasticProblem(genNames, meshes) {}
IStochasticTransportProblem(Meshes &meshes, GeneratorNames genNames = {})
: IStochasticProblem(meshes, genNames) {}
bool RHS() const { return rhs; }
......@@ -30,8 +30,7 @@ public:
class StochasticPollution1D : public IStochasticTransportProblem {
public:
StochasticPollution1D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
meshes) {}
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
......@@ -39,14 +38,12 @@ public:
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalVectorFieldSample(c);
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalScalarSample(face, c);
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
string Name() const override { return "StochasticPollution1D"; }
......@@ -59,8 +56,7 @@ private:
public:
StochasticPollutionCosHat1D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
meshes) {}
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
Point midPoint = Point(0.2, 0.0);
......@@ -71,14 +67,12 @@ public:
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalVectorFieldSample(c);
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalScalarSample(face, c);
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
string Name() const override { return "StochasticPollutionCosHat1D"; }
......@@ -90,8 +84,7 @@ private:
public:
StochasticPollution2D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
meshes) {}
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
......@@ -100,14 +93,12 @@ public:
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalVectorFieldSample(c);
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalScalarSample(face, c);
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
string Name() const override { return "StochasticPollution2D"; }
......@@ -118,7 +109,7 @@ class CircleWave2D : public IStochasticTransportProblem {
double cc = 0.5;
public:
explicit CircleWave2D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"DummyGenerator"}, meshes) {
IStochasticTransportProblem(meshes) {
}
double Solution(double t, const Point &x) const override {
......@@ -154,8 +145,7 @@ private:
public:
StochasticPollutionCosHat2D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator1D"},
meshes) {}
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator1D"}) {}
double Solution(double t, const Point &x) const override {
Point midPoint = Point(0.5, 0.8);
......@@ -166,14 +156,12 @@ public:
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalVectorFieldSample(c);
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalScalarSample(face, c);
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
string Name() const override { return "StochasticPollutionCosHat2D"; }
......@@ -198,8 +186,7 @@ private:
public:
StochasticPollutionMollifiedBar2D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
meshes) {}
IStochasticTransportProblem(meshes, GeneratorNames{"HybridFluxGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
if (abs(midPoint[0] - x[0]) <= 0.5 * xl &&
......@@ -219,14 +206,12 @@ public:
}
VectorField CellFlux(const cell &c, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalVectorFieldSample(c);
return this->genContainer.vectorFieldGenerator->EvalSample(c);
}
Scalar FaceNormalFlux(const cell &c, int face,
const VectorField &N, const Point &x) const override {
return this->generators.find("HybridFluxGenerator")->
second->EvalScalarSample(face, c);
return this->genContainer.scalarGenerator->EvalSample(face, c);
}
string Name() const override { return "StochasticPollutionMollifiedBar2D"; }
......@@ -235,8 +220,8 @@ public:
class StochasticInitialConditionsTransport1D : public IStochasticTransportProblem {
public:
StochasticInitialConditionsTransport1D(Meshes &meshes) :
IStochasticTransportProblem(GeneratorNames{"NormalDistributionGenerator"},
meshes) {}
IStochasticTransportProblem(meshes,
GeneratorNames{"NormalDistributionGenerator"}) {}
Scalar Solution(double t, const Point &x) const override {
if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
......@@ -266,8 +251,7 @@ public:
mu(VectorField(0.5, 0.5, 0.0)),
sigma(Tensor()),
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(GeneratorNames{"DummyGenerator"},
meshes) {
IStochasticTransportProblem(meshes) {
sigma[0][0] = 0.01;
sigma[0][1] = 0.0;
sigma[1][0] = 0.0;
......@@ -304,8 +288,7 @@ private:
public:
StochasticGaussHat2D(Meshes &meshes) :
vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
IStochasticTransportProblem(GeneratorNames{"DummyGenerator"},
meshes) {}
IStochasticTransportProblem(meshes) {}
Scalar Solution(double t, const Point &x) const override {
// Todo 2D Gaussfunction for t = 0
......
......@@ -52,14 +52,14 @@ void CirculantEmbedding1D::generateCoarseSample(const SampleID &id,
fineSample = nullptr;
}
ComplexField1D CirculantEmbedding1D::generateField(const SampleID &id) {
ComplexSequence1D CirculantEmbedding1D::generateField(const SampleID &id) {
generator.DrawSample(id);
ComplexField1D normalDist = generator.EvalComplexField1DSample();
ComplexField1D Xf;
ComplexSequence1D normalDist = generator.EvalSample();
ComplexSequence1D Xf;
fineComplexField.clear();
Xf.resize(numberOfCellsEmbedded);
ComplexField1D product;
ComplexSequence1D product;
product.resize(numberOfCellsEmbedded);
for (int i = 0; i < numberOfCellsEmbedded; i++) {
product[i] = sqrtEigenvalues[i] * normalDist[i];
......@@ -206,9 +206,9 @@ void CirculantEmbedding2D::generateCoarseSample(const SampleID &id,
fineSample = nullptr;
}
ComplexField2D CirculantEmbedding2D::generateField(const SampleID &id) {
ComplexField2D Xf;
ComplexField2D Z;
ComplexSequence2D CirculantEmbedding2D::generateField(const SampleID &id) {
ComplexSequence2D Xf;
ComplexSequence2D Z;
Z.resize(numberOfCellsEmbedded[1]);
fineComplexField.clear();
Xf.resize(numberOfCellsEmbedded[1]);
......@@ -217,11 +217,11 @@ ComplexField2D CirculantEmbedding2D::generateField(const SampleID &id) {
Xf[i].resize(numberOfCellsEmbedded[0]);
for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
normalDist.DrawSample(id);
Z[i][j] = normalDist.EvalComplexSample();
Z[i][j] = normalDist.EvalSample();
}
}
ComplexField2D product;
ComplexSequence2D product;
product.resize(numberOfCellsEmbedded[1]);
for (int i = 0; i < numberOfCellsEmbedded[1]; i++) {
product[i].resize(numberOfCellsEmbedded[0]);
......
......@@ -69,9 +69,9 @@ public:
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public:
NormalDistributionComplex1DField generator;
NormalDistributionComplexSequence1D generator;
ComplexField1D fineComplexField;
ComplexSequence1D fineComplexField;
SqrtEigenValues1D sqrtEigenvalues;
......@@ -100,7 +100,7 @@ public:
Vector *&coarseSample) override;
//private:
ComplexField1D generateField(const SampleID &id);
ComplexSequence1D generateField(const SampleID &id);
void createCovarianceMatrix(ToeplitzRow &ar,
ToeplitzColumn &ac);
......@@ -120,7 +120,7 @@ public:
class CirculantEmbedding2D : public CirculantEmbedding, public CovarianceFunction2D {
public:
ComplexField2D fineComplexField;
ComplexSequence2D fineComplexField;
SqrtEigenValues2D sqrtEigenvalues;
......@@ -150,7 +150,7 @@ public:
Vector *&coarseSample) override;
//private:
ComplexField2D generateField(const SampleID &id);
ComplexSequence2D generateField(const SampleID &id);
void createCovarianceMatrix(BlockToeplitzRow &ar,
BlockToeplitzColumn &ac);
......@@ -168,7 +168,7 @@ public:
SqrtEigenValues2D getSqrtEV();
};
class MultilevelCirculantEmbedding : public SampleGenerator {
class MultilevelCirculantEmbedding : public SampleGenerator<Tensor> {
private:
PlotMap plotMap;
......@@ -212,7 +212,7 @@ public:
circulantEmbeddings.clear();
}
Tensor EvalTensorSample(const cell &c) override {
Tensor EvalSample(const cell &c) override {
if (fineSample && !coarseSample) return fineSample->operator()(c(), 0) * One;
if (coarseSample && !fineSample) return coarseSample->operator()(c(), 0) * One;
Exit("Error in Generator")
......
......@@ -8,121 +8,121 @@
#include "pdesolver/PDESolver.hpp"
class HybridFluxGenerator : public SampleGenerator {
private:
PDESolver *pdeSolver;
public:
HybridFluxGenerator(Meshes &meshes) : SampleGenerator(meshes) {
std::string problemName = "StochasticLaplace2D";
pdeSolver = PDESolverCreator(meshes).
WithModel("HybridElliptic").
WithProblem(problemName).Create();
}
};
class MultilevelHybridFluxGenerator : public SampleGenerator {
private:
MatrixGraphs cellMGraphs;
MatrixGraphs faceMGraphs;
Vector *faceFlux = nullptr;
Vector *faceValues = nullptr;
Vector *cellFlux = nullptr;
protected:
IStochasticEllipticProblem *problem;
HybridEllipticAssemble *assemble;
RTLagrangeDiscretization *disc;
Preconditioner *pc;
Solver *solver;
Newton *newton;
Meshes &meshes;
void drawSample(const SampleID &id) override {
if (!id.coarse) {
generateFineSample(id);
//class HybridFluxGenerator : public SampleGenerator<Scalar> {
//private:
// PDESolver *pdeSolver;
//
//public:
// HybridFluxGenerator(Meshes &meshes) : SampleGenerator(meshes) {
// std::string problemName = "StochasticLaplace2D";
// pdeSolver = PDESolverCreator(meshes).
// WithModel("HybridElliptic").
// WithProblem(problemName).Create();
// }
//};
//class MultilevelHybridFluxGenerator : public SampleGenerator {
//private:
// MatrixGraphs cellMGraphs;
// MatrixGraphs faceMGraphs;
//
// Vector *faceFlux = nullptr;
// Vector *faceValues = nullptr;
// Vector *cellFlux = nullptr;
//
//protected:
// IStochasticEllipticProblem *problem;
// HybridEllipticAssemble *assemble;
// RTLagrangeDiscretization *disc;
// Preconditioner *pc;
// Solver *solver;
// Newton *newton;
// Meshes &meshes;
//
// void drawSample(const SampleID &id) override {
// if (!id.coarse) {
// generateFineSample(id);
// if (plotting)
// plotter->PlotVector("flux", *cellFlux, 3,
// id.level.fine, "CellData");
} else {
generateCoarseSample(id);
// } else {
// generateCoarseSample(id);
// if (plotting)
// plotter->PlotVector("flux", *cellFlux, 3,
// id.level.coarse, "CellData");