template based stochastic problems

mlmc/src/problem/CMakeLists.txt

 set(problem
+        problem/StochasticProblem.cpp
+        problem/StochasticEllipticProblem.cpp
         )
\ No newline at end of file

mlmc/src/problem/StochasticEllipticProblem.cpp

+#include "problem/StochasticEllipticProblem.hpp"
+
+
+template<>
+Tensor StochasticLaplace1DT<Vector>::Permeability(const cell &c) const {
+    mout << *this->fineSample << endl;
+    mout << this->fineSample->size() << endl;
+    return this->fineSample->operator()(c(), 0) * One;
+}
+
+template<>
+Tensor StochasticLaplace2DT<Vector>::Permeability(const cell &c) const {
+    return this->fineSample->operator()(c(), 0) * One;
+}
\ No newline at end of file

mlmc/src/problem/StochasticEllipticProblem.hpp

@@ -4,11 +4,11 @@
 #include "problem/StochasticProblem.hpp"
 
 
-class StochasticEllipticProblem : public StochasticProblem {
+template<typename Sample=Vector>
+class IStochasticEllipticProblemT : public IStochasticProblemT<Sample> {
 public:
-    StochasticEllipticProblem() = default;
-
-    ~StochasticEllipticProblem() override {};
+    IStochasticEllipticProblemT(GeneratorNames generatorNames, Meshes &meshes)
+        : IStochasticProblemT<Sample>(generatorNames, meshes) {}
 
     virtual Scalar Solution(const Point &x) const = 0;
 
@@ -16,21 +16,23 @@ public:
 
     virtual Scalar Load(int, const Point &x) const = 0;
 
-    virtual Tensor Permeability(const cell &c) const {
-        if(fieldSample == nullptr)
-            Exit("\nNo permeability sample has been set!\n")
-        return (*fieldSample)(c(), 0) * One;
-    }
+    virtual Tensor Permeability(const cell &c) const = 0;
 
+    // Todo remove
     virtual Tensor InversePermeability(const cell &c) const {
         Tensor permeabilityTensor = Permeability(c);
         return Invert(permeabilityTensor);
     }
 };
 
-class StochasticLaplace1D : public StochasticEllipticProblem {
+typedef IStochasticEllipticProblemT<> IStochasticEllipticProblem;
+
+template<typename Sample=Vector>
+class StochasticLaplace1DT : public IStochasticEllipticProblemT<Sample> {
 public:
-    StochasticLaplace1D() = default;
+    StochasticLaplace1DT(Meshes &meshes) :
+        IStochasticEllipticProblemT<Sample>(GeneratorNames{"CirculantEmbedding1D"},
+                                            meshes) {};
 
     Scalar Solution(const Point &x) const override { return 1 - x[0]; }
 
@@ -40,12 +42,19 @@ public:
 
     Scalar Load(int, const Point &x) const override { return 0.0; }
 
-    string Name() const override { return "Stochastic Laplace 1D"; }
+    Tensor Permeability(const cell &c) const override;
+
+    string Name() const override { return "StochasticLaplace1D"; }
 };
 
-class StochasticLaplace2D : public StochasticEllipticProblem {
+typedef StochasticLaplace1DT<> StochasticLaplace1D;
+
+template<typename Sample=Vector>
+class StochasticLaplace2DT : public IStochasticEllipticProblemT<Sample> {
 public:
-    StochasticLaplace2D() = default;
+    StochasticLaplace2DT(Meshes &meshes) :
+        IStochasticEllipticProblemT<Sample>(GeneratorNames{"CirculantEmbedding2D"},
+                                            meshes) {}
 
     Scalar Solution(const Point &x) const override { return -x[1]; }
 
@@ -55,7 +64,11 @@ public:
 
     Scalar Load(int, const Point &x) const override { return 0.0; }
 
-    string Name() const override { return "Stochastic Laplace 2D"; }
+    Tensor Permeability(const cell &c) const override;
+
+    string Name() const override { return "StochasticLaplace2D"; }
 };
 
+typedef StochasticLaplace2DT<> StochasticLaplace2D;
+
 #endif //STOCHASTICELLIPTICPROBLEM_HPP

mlmc/src/problem/StochasticProblem.cpp

+#include "StochasticProblem.hpp"
+#include "stochastics/CirculantEmbedding.hpp"
+#include "stochastics/HybridFluxGenerator.hpp"
+
+
+template<typename Sample>
+SampleGenerator<Sample> *IStochasticProblemT<Sample>
+::createGenerator(int l, GeneratorName generatorName) {
+    if (generatorName == "CirculantEmbedding1D")
+        return new CirculantEmbedding1D(l, meshes);
+    if (generatorName == "CirculantEmbedding2D")
+        return new CirculantEmbedding2D(l, meshes);
+//    if (generatorName == "HybridFluxGenerator1D")
+//        return new HybridFluxGenerator1D(l, meshes);
+//    if (generatorName == "HybridFluxGenerator2D")
+//        return new HybridFluxGenerator2D(l, meshes);
+    Exit("Generator not implemented")
+}
+
+template<>
+SampleGenerator<Vector> *IStochasticProblemT<Vector>
+::createGenerator(int l, GeneratorName generatorName) {
+    if (generatorName == "CirculantEmbedding1D")
+        return new CirculantEmbedding1D(l, meshes);
+    if (generatorName == "CirculantEmbedding2D")
+        return new CirculantEmbedding2D(l, meshes);
+//    if (generatorName == "HybridFluxGenerator1D")
+//        return new HybridFluxGenerator1D(l, meshes);
+//    if (generatorName == "HybridFluxGenerator2D")
+//        return new HybridFluxGenerator2D(l, meshes);
+    Exit("Generator not implemented")
+}
\ No newline at end of file

mlmc/src/problem/StochasticProblem.hpp

@@ -3,28 +3,54 @@
 
 #include "Algebra.h"
 #include "montecarlo/SampleID.hpp"
+#include "stochastics/SampleGenerator.hpp"
 
 
-class StochasticProblem {
+typedef std::string GeneratorName;
+
+typedef std::vector<GeneratorName> GeneratorNames;
+
+template<typename Sample=Vector>
+class IStochasticProblemT {
+private:
+    SampleGenerator<Sample> *createGenerator(int l, GeneratorName generatorName);
+
+    Meshes &meshes; // Todo: should be removed at some point
 protected:
-    Vector *fieldSample = nullptr;
-public:
+    typedef std::map<int, std::vector<SampleGenerator<Sample> *>> SampleGenerators;
 
-    StochasticProblem() = default;
+    SampleGenerators generators;
 
-    virtual ~StochasticProblem() = default;;
+    Sample *fineSample = nullptr;
+    Sample *coarseSample = nullptr;
 
-    void SetSample(Vector *sample) {
-        fieldSample = sample;
+public:
+    IStochasticProblemT(GeneratorNames generatorNames, Meshes &meshes) :
+        meshes(meshes) {
+        for (int l = meshes.pLevel(); l <= meshes.Level(); l++)
+            for (auto &name : generatorNames)
+                generators[l].push_back(createGenerator(l, name));
     }
 
+    virtual ~IStochasticProblemT() {
+        delete fineSample;
+        delete coarseSample;
+
+        for (auto &[level, gens]: generators)
+            for (auto &gen : gens)
+                delete gen;
+    };
+
     void DrawSample(SampleID id) {
-        // Sagt allen Generatoren sie sollen Sample erzeugen
-        // Falls in SampleID.coarse == true
-        //   dann vorgehaltenes sample
+        for (auto &gens : generators[id.level])
+            gens->DrawSample(id, fineSample, coarseSample);
     }
 
+    virtual void PrintInfo() const {};  // Todo make pure virtual at some point
+
     virtual string Name() const = 0;
 };
 
+typedef IStochasticProblemT<> IStochasticProblem;
+
 #endif

mlmc/src/problem/StochasticReactionProblem.hpp

@@ -5,11 +5,13 @@
 #include "StochasticHybridFlux.hpp"
 
 
-class StochasticReactionProblem : public StochasticProblem {
+template<typename Sample>
+class StochasticReactionProblem : public StochasticProblem<Sample> {
 protected:
     double convection = 1.0;
     double diffusion = 0.001;
     double reaction = 1.0;
+
 public:
     StochasticReactionProblem() {
         config.get("Diffusion", diffusion);
@@ -83,10 +85,11 @@ inline Scalar u0(const Point &x) {
     else return 0.0;
 }
 
-class ExponentialReaction2D : public StochasticReactionProblem {
+template<typename Sample>
+class ExponentialReaction2D : public StochasticReactionProblem<Sample> {
 public:
 
-    ExponentialReaction2D() : StochasticReactionProblem() {};
+    ExponentialReaction2D() : StochasticReactionProblem<Sample>() {};
 
     Scalar Concentration(double t, const Point &x) const override {
         return u0(x);
@@ -97,21 +100,22 @@ public:
     }
 
     Tensor Diffusion(int sd, const Point &x) const override {
-        return diffusion * One;
+        return this->diffusion * One;
     }
 
     Scalar Reaction(const Point &, double c) const override {
-        return reaction * c;
+        return this->reaction * c;
     }
 
     Scalar DerivativeReaction(const Point &, double c) const override {
-        return reaction;
+        return this->reaction;
     }
 
     string Name() const override { return "ExponentialReaction2D"; }
 };
 
-class LogisticReaction2D : public StochasticReactionProblem {
+template<typename Sample>
+class LogisticReaction2D : public StochasticReactionProblem<Sample> {
 private:
     double r0 = 1.0;
     double r1 = 0.0;
@@ -130,7 +134,7 @@ public:
     }
 
     Tensor Diffusion(int sd, const Point &x) const override {
-        return diffusion * One;
+        return this->diffusion * One;
     }
 
     Scalar Reaction(const Point &x, double c) const override {
@@ -144,32 +148,34 @@ public:
     string Name() const override { return "LogisticReaction2D"; }
 };
 
-class HybridExponentialReaction2D : public ExponentialReaction2D,
+template<typename Sample>
+class HybridExponentialReaction2D : public ExponentialReaction2D<Sample>,
                                     public StochasticHybridFlux {
 
     VectorField CellConvection(int sd, const cell &c, const Point &x) const override {
-        return convection * hybrid->EvaluateCellFlux(*flux, c);
+        return this->convection * hybrid->EvaluateCellFlux(*flux, c);
     }
 
     double FaceConvection(int sd, const cell &c, int f,
                           const VectorField &N, const Point &x) const override {
-        return convection * hybrid->EvaluateNormalFlux(*flux, c, f);
+        return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
     }
 
     string Name() const override { return "HybridExponentialReaction2D"; }
 };
 
-class HybridLogisticReaction2D : public LogisticReaction2D,
+template<typename Sample>
+class HybridLogisticReaction2D : public LogisticReaction2D<Sample>,
                                  public StochasticHybridFlux {
 
     VectorField CellConvection(int sd, const cell &c,
                                const Point &x) const override {
-        return convection * hybrid->EvaluateCellFlux(*flux, c);
+        return this->convection * hybrid->EvaluateCellFlux(*flux, c);
     }
 
     double FaceConvection(int sd, const cell &c, int f,
                           const VectorField &N, const Point &x) const override {
-        return convection * hybrid->EvaluateNormalFlux(*flux, c, f);
+        return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
     }
 
     string Name() const override { return "HybridLogisticReaction2D"; }

mlmc/src/problem/StochasticTransportProblem.hpp

@@ -5,18 +5,8 @@
 #include "StochasticHybridFlux.hpp"
 
 
-class StochasticTransportProblem : public StochasticProblem {
-    /*
-     * Stochastic Transport Problem can handle:
-     *  - Random right hand sides   ->
-     *  - Random initial conditions ->
-     *  - Random inflow conditions  ->
-     *  - Random coefficients       -> StochasticHybridFluxTransport
-     *              -> warum hat der hier nicht einfach einen generator
-     *
-     * Design your problem by deriving from
-     *
-     */
+template<typename Sample=Vector>
+class IStochasticTransportProblemT : public IStochasticProblemT<Sample> {
 protected:
     VectorField fluxField = Point(0.0, -1.0, 0.0);
     bool exactAverage = false;  // Todo maybe
@@ -24,11 +14,11 @@ protected:
 public:
     bool rhs = false;
 
-    explicit StochasticTransportProblem() {
+    explicit IStochasticTransportProblemT() {
         // Deal with b
     }
 
-    virtual ~StochasticTransportProblem() = default;
+    virtual ~IStochasticTransportProblemT() = default;
 
     bool RHS() const { return rhs; }
 
@@ -52,7 +42,10 @@ public:
     }
 };
 
-class StochasticPollution1D : public StochasticTransportProblem,
+typedef IStochasticTransportProblemT<> IStochasticTransportProblem;
+
+template<typename Sample=Vector>
+class StochasticPollution1D : public IStochasticTransportProblemT<Sample>,
                               public StochasticHybridFlux {
 public:
     StochasticPollution1D() {}
@@ -76,7 +69,8 @@ public:
     string Name() const override { return "StochasticPollution1D"; }
 };
 
-class StochasticPollutionCosHat1D : public StochasticTransportProblem,
+template<typename Sample=Vector>
+class StochasticPollutionCosHat1D : public IStochasticTransportProblemT<Sample>,
                                     public StochasticHybridFlux {
 private:
     double amplitude = 1.00;
@@ -106,7 +100,8 @@ public:
     string Name() const override { return "StochasticPollutionCosHat1D"; }
 };
 
-class StochasticPollution2D : public StochasticTransportProblem,
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticPollution2D : public IStochasticTransportProblemT<Sample>,
                               public StochasticHybridFlux {
 public:
     double ut(double t, const Point &x) const override {
@@ -128,7 +123,8 @@ public:
     string Name() const override { return "Stochastic Pollution 2D"; }
 };
 
-class StochasticPollutionCosHat2D : public StochasticTransportProblem,
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticPollutionCosHat2D : public IStochasticTransportProblemT<Sample>,
                                     public StochasticHybridFlux {
 private:
     double amplitude = 1.00;
@@ -157,7 +153,8 @@ public:
     string Name() const override { return "StochasticPollutionCosHat2D"; }
 };
 
-class StochasticPollutionMollifiedBar2D : public StochasticTransportProblem,
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticPollutionMollifiedBar2D : public IStochasticTransportProblemT<Sample>,
                                           public StochasticHybridFlux {
 private:
     double ee = 0.12;
@@ -208,11 +205,12 @@ public:
     string Name() const override { return "Stochastic Pollution 2D"; }
 };
 
-class StochasticInitialConditionsTransport1D : public StochasticTransportProblem {
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticInitialConditionsTransport1D : public IStochasticTransportProblemT<Sample> {
     double c = 1.0;
 public:
     explicit StochasticInitialConditionsTransport1D() {
-        rhs = false;
+        this->rhs = false;
     }
 
     double Amplitude(double r) const {
@@ -231,9 +229,10 @@ public:
     string Name() const override { return "RiemannTransport1D"; }
 };
 
-class StochasticInflowTransport1D : public StochasticTransportProblem {
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticInflowTransport1D : public IStochasticTransportProblemT<Sample> {
 public:
-    explicit StochasticInflowTransport1D() : StochasticTransportProblem() {}
+    explicit StochasticInflowTransport1D() {}
 
     double ut(double t, const Point &x) const override {
         if (x[0] == 0.0 && x[1] >= 1.0) return 1.0;
@@ -243,7 +242,8 @@ public:
     string Name() const override { return "InflowTest"; }
 };
 
-class StochasticCircleWave2D : public StochasticTransportProblem {
+template<typename Sample=std::pair<Vector, Vector>>
+class StochasticCircleWave2D : public IStochasticTransportProblemT<Sample> {
     /*
      * Todo Idea: create stochastic flux
      */