updated problem idea

mlmc/src/problem/StochasticEllipticProblem.cpp

 #include "problem/StochasticEllipticProblem.hpp"
-
-
-template<>
-Tensor StochasticLaplace1DT<Vector>::Permeability(const cell &c) const {
-    if (fineSample && !coarseSample)
-        return this->fineSample->operator()(c(), 0) * One;
-    else if (coarseSample && !fineSample)
-        return this->coarseSample->operator()(c(), 0) * One;
-    else
-        Exit("Entry in sample not found")
-}
-
-template<>
-Tensor StochasticLaplace2DT<Vector>::Permeability(const cell &c) const {
-    if (fineSample && !coarseSample)
-        return this->fineSample->operator()(c(), 0) * One;
-    else if (coarseSample && !fineSample)
-        return this->coarseSample->operator()(c(), 0) * One;
-    else
-    Exit("Entry in sample not found")
-}
\ No newline at end of file

mlmc/src/problem/StochasticEllipticProblem.hpp

@@ -4,65 +4,68 @@
 #include "problem/StochasticProblem.hpp"
 
 
-template<typename Sample=Vector>
-class IStochasticEllipticProblemT : public IStochasticProblemT<Sample> {
+class IStochasticEllipticProblem : public IStochasticProblem {
 public:
-    IStochasticEllipticProblemT(GeneratorNames generatorNames, Meshes &meshes)
-        : IStochasticProblemT<Sample>(generatorNames, meshes) {}
+    IStochasticEllipticProblem(GeneratorNames genNames, Meshes &meshes)
+        : IStochasticProblem(genNames, meshes) {}
 
     virtual Scalar Solution(const Point &x) const = 0;
 
-    virtual VectorField Flux(int, const Point &x) const = 0;
+    virtual VectorField Flux(const Point &x) const = 0;
 
-    virtual Scalar Load(int, const Point &x) const = 0;
+    virtual Scalar Load(const Point &x) const = 0;
 
     virtual Tensor Permeability(const cell &c) const = 0;
 };
 
-typedef IStochasticEllipticProblemT<> IStochasticEllipticProblem;
-
-template<typename Sample=Vector>
-class StochasticLaplace1DT : public IStochasticEllipticProblemT<Sample> {
+class StochasticLaplace1D : public IStochasticEllipticProblem {
 public:
-    StochasticLaplace1DT(Meshes &meshes) :
-        IStochasticEllipticProblemT<Sample>(GeneratorNames{"CirculantEmbedding1D"},
-                                            meshes) {};
+    StochasticLaplace1D(Meshes &meshes) :
+        IStochasticEllipticProblem(GeneratorNames{"CirculantEmbedding"},
+                                   meshes) {};
 
-    Scalar Solution(const Point &x) const override { return 1 - x[0]; }
+    Scalar Solution(const Point &x) const override {
+        return 1 - x[0];
+    }
 
-    VectorField Flux(int, const Point &x) const override {
+    VectorField Flux(const Point &x) const override {
         return VectorField(1.0, 0.0);
     }
 
-    Scalar Load(int, const Point &x) const override { return 0.0; }
+    Scalar Load(const Point &x) const override {
+        return 0.0;
+    }
 
-    Tensor Permeability(const cell &c) const override;
+    Tensor Permeability(const cell &c) const override {
+        return this->generators.find("CirculantEmbedding")->
+            second->EvalTensorSample(c);
+    }
 
     string Name() const override { return "StochasticLaplace1D"; }
 };
 
-typedef StochasticLaplace1DT<> StochasticLaplace1D;
-
-template<typename Sample=Vector>
-class StochasticLaplace2DT : public IStochasticEllipticProblemT<Sample> {
+class StochasticLaplace2D : public IStochasticEllipticProblem {
 public:
-    StochasticLaplace2DT(Meshes &meshes) :
-        IStochasticEllipticProblemT<Sample>(GeneratorNames{"CirculantEmbedding2D"},
-                                            meshes) {}
+    StochasticLaplace2D(Meshes &meshes) :
+        IStochasticEllipticProblem(GeneratorNames{"CirculantEmbedding"},
+                                   meshes) {}
 
-    Scalar Solution(const Point &x) const override { return -x[1]; }
+    Scalar Solution(const Point &x) const override {
+        return -x[1];
+    }
 
-    VectorField Flux(int, const Point &x) const override {
+    VectorField Flux(const Point &x) const override {
         return VectorField(0.0, -1.0);
     }
 
-    Scalar Load(int, const Point &x) const override { return 0.0; }
+    Scalar Load(const Point &x) const override { return 0.0; }
 
-    Tensor Permeability(const cell &c) const override;
+    Tensor Permeability(const cell &c) const override {
+        return this->generators.find("CirculantEmbedding")->
+            second->EvalTensorSample(c);
+    }
 
     string Name() const override { return "StochasticLaplace2D"; }
 };
 
-typedef StochasticLaplace2DT<> StochasticLaplace2D;
-
 #endif //STOCHASTICELLIPTICPROBLEM_HPP

mlmc/src/problem/StochasticProblem.cpp

@@ -3,30 +3,11 @@
 #include "stochastics/HybridFluxGenerator.hpp"
 
 
-template<typename Sample>
-SampleGenerator<Sample> *IStochasticProblemT<Sample>
-::createGenerator(int l, GeneratorName genName) {
-    if (genName == "CirculantEmbedding1D")
-        return new CirculantEmbedding1D(l, meshes);
-    if (genName == "CirculantEmbedding2D")
-        return new CirculantEmbedding2D(l, meshes);
-//    if (generatorName == "HybridFluxGenerator1D")
-//        return new HybridFluxGenerator1D(l, meshes);
-//    if (generatorName == "HybridFluxGenerator2D")
-//        return new HybridFluxGenerator2D(l, meshes);
+SampleGenerator *IStochasticProblem::createGenerator(GeneratorName genName,
+                                                     Meshes &meshes) {
+    if (genName == "CirculantEmbedding")
+        return new MultilevelCirculantEmbedding(meshes);
+    if (genName == "HybridFluxGenerator")
+        return new MultilevelHybridFluxGenerator(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

@@ -2,7 +2,7 @@
 #define STOCHASTICPROBLEM_HPP
 
 #include "Algebra.h"
-#include "montecarlo/SampleID.hpp"
+#include "utility/SampleID.hpp"
 #include "stochastics/SampleGenerator.hpp"
 #include "main/Utils.hpp"
 
@@ -11,47 +11,37 @@ typedef std::string GeneratorName;
 
 typedef std::vector<GeneratorName> GeneratorNames;
 
-template<typename Sample=Vector>
-class IStochasticProblemT {
+typedef std::map<GeneratorName, SampleGenerator *> SampleGenerators;
+
+class IStochasticProblem {
 private:
     int verbose = 1;
 
     GeneratorNames genNames;
 
-    SampleGenerator<Sample> *createGenerator(int l, GeneratorName genName);
-
-    Meshes &meshes;
+    SampleGenerator *createGenerator(GeneratorName genName, Meshes &meshes);
 
 protected:
-    typedef std::map<int, std::vector<SampleGenerator<Sample> *>> SampleGenerators;
-
     SampleGenerators generators;
 
-    Sample *fineSample = nullptr;
-
-    Sample *coarseSample = nullptr;
-
 public:
-    IStochasticProblemT(GeneratorNames genNames, Meshes &meshes) :
-        genNames(genNames), meshes(meshes) {
-        config.get("problemVerbose", verbose);
-        for (int l = meshes.pLevel(); l <= meshes.Level(); l++)
-            for (auto &name : genNames)
-                generators[l].push_back(createGenerator(l, name));
+    IStochasticProblem(GeneratorNames genNames, Meshes &meshes) : genNames(genNames) {
+        config.get("ProblemVerbose", verbose);
+
+        for (auto &name : genNames)
+            generators[name] = createGenerator(name, meshes);
     }
 
-    virtual ~IStochasticProblemT() {
-        if (fineSample) delete fineSample;
-        if (coarseSample) delete coarseSample;
-        for (auto &[level, gens]: generators)
-            for (auto &gen : gens)
-                delete gen;
+    virtual ~IStochasticProblem() {
+        for (auto &gen: generators) {
+            delete gen.second;
+        }
         generators.clear();
     };
 
     void DrawSample(SampleID id) {
-        for (auto &gens : generators[id.level])
-            gens->DrawSample(id, fineSample, coarseSample);
+        for (auto &gen : generators)
+            gen.second->DrawSample(id);
     }
 
     virtual void PrintInfo() const {
@@ -60,17 +50,7 @@ public:
                        PrintInfoEntry("Generators", vec2str(genNames)));
     };
 
-    Sample *GetFineSample() {
-        return fineSample;
-    }
-
-    Sample *GetCoarseSample() {
-        return coarseSample;
-    }
-
     virtual string Name() const = 0;
 };
 
-typedef IStochasticProblemT<> IStochasticProblem;
-
 #endif

mlmc/src/problem/StochasticReactionProblem.hpp

@@ -5,15 +5,15 @@
 #include "StochasticHybridFlux.hpp"
 
 
-template<typename Sample>
-class StochasticReactionProblem : public StochasticProblem<Sample> {
+class IStochasticReactionProblem : public IStochasticProblem {
 protected:
     double convection = 1.0;
     double diffusion = 0.001;
     double reaction = 1.0;
 
 public:
-    StochasticReactionProblem() {
+    IStochasticReactionProblem(GeneratorNames genNames, Meshes &meshes) :
+        IStochasticProblem(genNames, meshes) {
         config.get("Diffusion", diffusion);
         config.get("Convection", convection);
         config.get("Reaction", reaction);
@@ -31,13 +31,13 @@ public:
         return reaction;
     }
 
-    virtual ~StochasticReactionProblem() = default;
+    virtual ~IStochasticReactionProblem() = default;
 
     virtual Scalar Concentration(double t, const Point &x) const {
         return 0.0;
     }
 
-    virtual Scalar Load(int sd, const Point &x) const {
+    virtual Scalar Load(const Point &x) const {
         return 0.0;
     }
 
@@ -45,28 +45,25 @@ public:
         return zero;
     }
 
-    virtual VectorField Convection(int sd, const Point &x) const {
+    virtual VectorField Convection(const Point &x) const {
         return convection * zero;
     }
 
-    virtual VectorField CellConvection(int sd, const cell &c,
-                                       const Point &x) const {
-        return Convection(sd, x);
+    virtual VectorField CellConvection(const cell &c, const Point &x) const {
+        return Convection(x);
     }
 
-    virtual double FaceConvection(int sd, const cell &c, int f,
-                                  const VectorField &N, const Point &x) const {
-        return Convection(sd, x) * N;
+    virtual double FaceConvection(const cell &c,
+                                  int f,
+                                  const VectorField &N,
+                                  const Point &x) const {
+        return Convection(x) * N;
     }
 
-    virtual Tensor Diffusion(int sd, const Point &x) const {
+    virtual Tensor Diffusion(const Point &x) const {
         return diffusion * Zero;
     }
 
-    virtual double Diffusion(const Point &x) const {
-        return diffusion * 0.0;
-    }
-
     virtual Scalar Reaction(const Point &, double c) const {
         return reaction * 0.0;
     }
@@ -85,21 +82,21 @@ inline Scalar u0(const Point &x) {
     else return 0.0;
 }
 
-template<typename Sample>
-class ExponentialReaction2D : public StochasticReactionProblem<Sample> {
+class ExponentialReaction2D : public IStochasticReactionProblem {
 public:
 
-    ExponentialReaction2D() : StochasticReactionProblem<Sample>() {};
+    ExponentialReaction2D(Meshes &meshes) :
+        IStochasticReactionProblem(GeneratorNames{}, meshes) {};
 
     Scalar Concentration(double t, const Point &x) const override {
         return u0(x);
     }
 
-    VectorField Convection(int sd, const Point &x) const override {
+    VectorField Convection(const Point &x) const override {
         return VectorField(0.0, -1.0);
     }
 
-    Tensor Diffusion(int sd, const Point &x) const override {
+    Tensor Diffusion(const Point &x) const override {
         return this->diffusion * One;
     }
 
@@ -114,13 +111,13 @@ public:
     string Name() const override { return "ExponentialReaction2D"; }
 };
 
-template<typename Sample>
-class LogisticReaction2D : public StochasticReactionProblem<Sample> {
+class LogisticReaction2D : public IStochasticReactionProblem {
 private:
     double r0 = 1.0;
     double r1 = 0.0;
 public:
-    LogisticReaction2D() {
+    LogisticReaction2D(Meshes &meshes) :
+        IStochasticReactionProblem(GeneratorNames{}, meshes) {
         config.get("Reaction0", r0);
         config.get("Reaction1", r1);
     }
@@ -129,11 +126,11 @@ public:
         return u0(x);
     }
 
-    VectorField Convection(int sd, const Point &x) const override {
+    VectorField Convection(const Point &x) const override {
         return VectorField(0.0, -1.0);
     }
 
-    Tensor Diffusion(int sd, const Point &x) const override {
+    Tensor Diffusion(const Point &x) const override {
         return this->diffusion * One;
     }
 
@@ -148,33 +145,34 @@ public:
     string Name() const override { return "LogisticReaction2D"; }
 };
 
-template<typename Sample>
-class HybridExponentialReaction2D : public ExponentialReaction2D<Sample>,
+class HybridExponentialReaction2D : public ExponentialReaction2D,
                                     public StochasticHybridFlux {
 
-    VectorField CellConvection(int sd, const cell &c, const Point &x) const override {
+    VectorField CellConvection(const cell &c, const Point &x) const override {
         return this->convection * hybrid->EvaluateCellFlux(*flux, c);
     }
 
-    double FaceConvection(int sd, const cell &c, int f,
-                          const VectorField &N, const Point &x) const override {
+    double FaceConvection(const cell &c,
+                          int f,
+                          const VectorField &N,
+                          const Point &x) const override {
         return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
     }
 
     string Name() const override { return "HybridExponentialReaction2D"; }
 };
 
-template<typename Sample>
-class HybridLogisticReaction2D : public LogisticReaction2D<Sample>,
+class HybridLogisticReaction2D : public LogisticReaction2D,
                                  public StochasticHybridFlux {
 
-    VectorField CellConvection(int sd, const cell &c,
-                               const Point &x) const override {
+    VectorField CellConvection(const cell &c, const Point &x) const override {
         return this->convection * hybrid->EvaluateCellFlux(*flux, c);
     }
 
-    double FaceConvection(int sd, const cell &c, int f,
-                          const VectorField &N, const Point &x) const override {
+    double FaceConvection(const cell &c,
+                          int f,
+                          const VectorField &N,
+                          const Point &x) const override {
         return this->convection * hybrid->EvaluateNormalFlux(*flux, c, f);
     }
 

mlmc/src/problem/StochasticTransportProblem.hpp

 #ifndef MLMC_STOCHASTICTRANSPORTPROBLEM_HPP
 #define MLMC_STOCHASTICTRANSPORTPROBLEM_HPP
 
-#include "StochasticProblem.hpp"
-#include "StochasticHybridFlux.hpp"
+#include "problem/StochasticProblem.hpp"
 
 
-template<typename Sample=Vector>
-class IStochasticTransportProblemT : public IStochasticProblemT<Sample> {
-protected:
-    VectorField fluxField = Point(0.0, -1.0, 0.0);
-    bool exactAverage = false;  // Todo maybe
-
+class IStochasticTransportProblem : public IStochasticProblem {
 public:
-    bool rhs = false;
-
-    explicit IStochasticTransportProblemT() {
-        // Deal with b
-    }
-
-    virtual ~IStochasticTransportProblemT() = default;
-
-    bool RHS() const { return rhs; }
-
-    virtual double ut(double, const Point &x) const = 0;
-
-    virtual double u0(const Point &x) const { return ut(0, x); }
+    IStochasticTransportProblem(GeneratorNames genNames, Meshes &meshes)
+        : IStochasticProblem(genNames, meshes) {}
 
+    // Todo rmv
     virtual bool Dirichlet(const Point &x) const { return false; }
 
-    virtual VectorField B(const Point &x) const {
-        return fluxField;
-    }
+    virtual Scalar Solution(double t, const Point &x) const = 0;
 
-    virtual VectorField cellB(const cell &c, const Point &x) const {
-        return B(x);
-    }
+    // Todo cell and point necessary?
+    virtual VectorField CellFlux(const cell &c, const Point &x) const = 0;
 
-    virtual double faceB(const cell &c, int f, const VectorField &N,
-                         const Point &x) const {
-        return B(x) * N;
-    }
+    // Todo cell and point necessary?
+    virtual Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                                  const Point &x) const = 0;
 };
 
-typedef IStochasticTransportProblemT<> IStochasticTransportProblem;
-
-template<typename Sample=Vector>
-class StochasticPollution1D : public IStochasticTransportProblemT<Sample>,
-                              public StochasticHybridFlux {
+class StochasticPollution1D : public IStochasticTransportProblem {
 public:
-    StochasticPollution1D() {}
+    StochasticPollution1D(Meshes &meshes) :
+        IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
+                                    meshes) {}
 
-    double ut(double t, const Point &x) const override {
-        if (abs(1.0 * t - x[0] + 0.5) > 0.06251)
-            return 0.0;
-        else
-            return 1.0;
+    Scalar Solution(double t, const Point &x) const override {
+        if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
+        return 1.0;
     }
 
-    VectorField cellB(const cell &c, const Point &x) const override {
-        return hybrid->EvaluateCellFlux(*flux, c);
+    VectorField CellFlux(const cell &c, const Point &x) const override {
+        return this->generators.find("HybridFluxGenerator")->
+            second->EvalVectorFieldSample(c);
     }
 
-    double faceB(const cell &c, int f,
-                 const VectorField &N, const Point &x) const override {
-        return hybrid->EvaluateNormalFlux(*flux, c, f);
+    Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                          const Point &x) const override {
+        return this->generators.find("HybridFluxGenerator")->
+            second->EvalScalarSample(face, c);
     }
 
     string Name() const override { return "StochasticPollution1D"; }
 };
 
-template<typename Sample=Vector>
-class StochasticPollutionCosHat1D : public IStochasticTransportProblemT<Sample>,
-                                    public StochasticHybridFlux {
+class StochasticPollutionCosHat1D : public IStochasticTransportProblem {
 private:
     double amplitude = 1.00;
     double cc = 6.0;
 
 public:
+    StochasticPollutionCosHat1D(Meshes &meshes) :
+        IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
+                                    meshes) {}
 
-    StochasticPollutionCosHat1D() {}
-
-    double ut(double t, const Point &x) const override {
+    Scalar Solution(double t, const Point &x) const override {
         Point midPoint = Point(0.2, 0.0);
         double r = dist(midPoint, x);
         if (r < 1 / cc)
@@ -88,52 +65,60 @@ public:
         return 0.0;
     }
 
-    VectorField cellB(const cell &c, const Point &x) const override {
-        return hybrid->EvaluateCellFlux(*flux, c);
+    VectorField CellFlux(const cell &c, const Point &x) const override {
+        return this->generators.find("HybridFluxGenerator")->
+            second->EvalVectorFieldSample(c);
     }
 
-    double faceB(const cell &c, int f,
-                 const VectorField &N, const Point &x) const override {
-        return hybrid->EvaluateNormalFlux(*flux, c, f);
+    Scalar FaceNormalFlux(const cell &c, int face,
+                          const VectorField &N, const Point &x) const override {
+        return this->generators.find("HybridFluxGenerator")->
+            second->EvalScalarSample(face, c);
     }
 
     string Name() const override { return "StochasticPollutionCosHat1D"; }
 };
 
-template<typename Sample=std::pair<Vector, Vector>>
-class StochasticPollution2D : public IStochasticTransportProblemT<Sample>,
-                              public StochasticHybridFlux {
+class StochasticPollution2D : public IStochasticTransportProblem {
+private:
+    double d = 1.0 / 16.0;
+
 public:
-    double ut(double t, const Point &x) const override {
-        double d = 1.0 / 16.0;
+    StochasticPollution2D(Meshes &meshes) :
+        IStochasticTransportProblem(GeneratorNames{"HybridFluxGenerator"},
+                                    meshes) {}
+
+    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))
             return 1.0;
-        else return 0.0;
+        return 0.0;
     }
 
-    VectorField cellB(const cell &c, const Point &x) const override {
-        return hybrid->EvaluateCellFlux(*flux, c);
+    VectorField CellFlux(const cell &c, const Point &x) const override {
+        return this->generators.find("HybridFluxGenerator")->
+            second->EvalVectorFieldSample(c);
     }