introduced sample id, adpated to new samle generator and new stochastic field. minor other changes

mlmc/src/montecarlo/MonteCarlo.cpp

@@ -22,12 +22,10 @@ void MonteCarlo::Method() {
     mout.EndBlock(verbose == 0);
 }
 
-void MonteCarlo::SolvePDE(int l, int m, bool fine,
-                          double &valueQ, double &cost,
-                          Vector &solution) {
+void MonteCarlo::SolvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) {
     mout.StartBlock("Solve PDE");
-    vout(1) << "l=" << l << " m=" << m << " fine=" << fine << endl;
-    solvePDE(l, m, fine, valueQ, cost, solution);
+    vout(1) << id.Str() << endl;
+    solvePDE(id, valueQ, cost, solution);
     vout(1) << "Q=" << valueQ << " cost=" << cost << endl;
     mout.EndBlock(verbose == 0);
 }
@@ -60,31 +58,3 @@ void MonteCarlo::updateStatistic() {
     kurtosis = (avgY4 - 4.0 * avgY3 * avgY + 6.0 * avgY2 * avgY * avgY -
         3.0 * avgY * avgY * avgY * avgY) / varY / varY;
 }
-
-string MonteCarlo::mkSampleDirName(int m, bool fineLevel) {
-    if (fineLevel)
-        return buildName(m, "sample");
-    else
-        return buildName(m, "sample", "coarse");
-}
-
-void MonteCarlo::mkSampleDir(int m, const string &sampleDirName) {
-    string pathToData = "data/vtk/";
-    char *sampleDirNameCharArr =
-        const_cast<char *> ((pathToData.append(sampleDirName)).c_str());
-    mkdir(sampleDirNameCharArr, 0777);
-}
-
-string MonteCarlo::buildName(int m, const string &name, const string &suffix) {
-    string returnName = name;
-    if (!suffix.empty())
-        returnName = returnName.append("_").append(suffix);
-    if (l != -1)
-        returnName = returnName.append("_").append(to_string(l));
-    if (m != -1)
-        returnName = returnName.append("_").append(to_string(m));
-    return returnName;
-}
-
-
-

mlmc/src/montecarlo/MonteCarlo.hpp

@@ -2,13 +2,6 @@
 #define MLMC_MC_HPP
 
 #include "main/MultilevelPlotter.hpp"
-#include "stochastics/StochasticField.hpp"
-
-
-struct Statistics {
-
-
-};
 
 
 class MonteCarlo {
@@ -22,7 +15,6 @@ public:
     string functional = "L2";
 
     Meshes *meshes = nullptr;
-    StochasticField *stochasticField = nullptr;
 
     int M = 0, dM = 0;
 
@@ -38,43 +30,33 @@ public:
     double fineQ = 0.0, coarseQ = 0.0;
     double fineCost = 0.0, coarseCost = 0.0;
 
-    MonteCarlo(int l, int dM, bool onlyFineLevel,
-               Meshes *meshes, StochasticField *stochFields) :
+    MonteCarlo(int l, int dM, bool onlyFineLevel, Meshes *meshes) :
         l(l), dM(dM), onlyFineLevel(onlyFineLevel),
-        meshes(meshes), pLevel(meshes->pLevel()),
-        stochasticField(stochFields) {
+        meshes(meshes), pLevel(meshes->pLevel()) {
 
         config.get("MCPlotting", plotting);
         config.get("MCVerbose", verbose);
         config.get("Functional", functional);
     }
 
+    virtual ~MonteCarlo() {};
+
     void Method();
 
-    void SolvePDE(int l, int m, bool fine,
-                  double &valueQ, double &cost,
-                  Vector &solution);
+    void SolvePDE(SampleID id, double &valueQ, double &cost, Vector &solution);
 
     virtual void Initialize() = 0;
 
 protected:
     virtual void method() = 0;
 
-    virtual void solvePDE(int l, int m, bool fine,
-                          double &valueQ, double &cost,
-                          Vector &solution) = 0;
+    virtual void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) = 0;
 
     void updateSums(double fineCost, double diffQ, double fineQ);
 
     void updateAvg();
 
     void updateStatistic();
-
-    string mkSampleDirName(int m, bool fineLevel);
-
-    static void mkSampleDir(int m, const string &sampleDirName);
-
-    string buildName(int m, const string &name, const string &suffix = "");
 };
 
 #endif //MLMC_MC_HPP
\ No newline at end of file

mlmc/src/montecarlo/MonteCarloElliptic.cpp

@@ -4,53 +4,43 @@
 using namespace std;
 
 void MonteCarloElliptic::Initialize() {
-    fineInput = 0.0, coarseInput = 0.0;
     fineSolution = 0.0, coarseSolution = 0.0;
     fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
 }
 
 void MonteCarloElliptic::method() {
     Initialize();
+    SampleID id;
     for (int m = M; m < M + dM; m++) {
-        string sampleDirName = mkSampleDirName(m, true);
-        plotter->SetDirectory(sampleDirName);
+        id = SampleID(l, m, false);
+        plotter->SetDirectory(id);
+        assemble->problem->DrawSample(id);
+        SolvePDE(id, fineQ, fineCost, fineSolution);
+
+        if (!onlyFineLevel) {
+            id = SampleID(l, m, false);
+            plotter->SetDirectory(id);
+            assemble->problem->DrawSample(id);
+            SolvePDE(SampleID(), coarseQ, coarseCost, coarseSolution);
+        } else
+            coarseQ = 0.0, coarseCost = 0.0;
 
-        stochasticField->GetFineSample(l, fineInput);
-        assemble->problem->SetSample(&fineInput);
-        SolvePDE(l, m, true, fineQ, fineCost, fineSolution);
+        updateSums(fineCost, fineQ - coarseQ, fineQ);
 
-        if (plotting)
+        if (plotting) {
             plotter->PlotVector("u", fineSolution,
                                 1, l, "VertexData");
-
-        if (onlyFineLevel)
-            coarseQ = 0.0, coarseCost = 0.0;
-
-        else {
-            sampleDirName = mkSampleDirName(m, false);
-            plotter->SetDirectory(sampleDirName);
-
-            stochasticField->GetCoarseSample(l, fineInput, coarseInput);
-            assemble->problem->SetSample(&coarseInput);
-            SolvePDE(l, 0, false,
-                     coarseQ, coarseCost, coarseSolution);
-
-            if (plotting)
-                plotter->PlotVector("u", coarseSolution,
-                                    1, l - 1, "VertexData");
+            plotter->PlotVector("u", coarseSolution,
+                                1, l - 1, "VertexData");
         }
-        updateSums(fineCost, fineQ - coarseQ, fineQ);
     }
 }
 
-void MonteCarloElliptic::solvePDE(int l, int m, bool fineLevel,
-                                  double &valueQ, double &cost, Vector &solution) {
-    // Todo refactor
-    Preconditioner *pc = GetPC("SuperLU");
-    Solver solver(pc, "GMRES");
-    Newton newton(solver);
-
-    newton.operator()(*assemble, solution);
+void MonteCarloElliptic::solvePDE(SampleID id,
+                                  double &valueQ,
+                                  double &cost,
+                                  Vector &solution) {
+    newton->operator()(*assemble, solution);
     cost = solution.pSize();
 
     if (functional == "L2")

mlmc/src/montecarlo/MonteCarloElliptic.hpp

@@ -3,6 +3,9 @@
 
 #include "MonteCarlo.hpp"
 #include "assemble/LagrangeEllipticAssemble.hpp"
+#include "solver/Solver.h"
+#include "solver/Newton.h"
+#include "dof/BasicDoFs.hpp"
 
 
 class MonteCarloElliptic : public MonteCarlo {
@@ -11,40 +14,37 @@ protected:
 
     void method() override;
 
-    void solvePDE(int l,
-                  int m,
-                  bool fineLevel,
-                  double &valueQ,
-                  double &cost,
-                  Vector &solution) override;
+    void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) override;
 
 public:
     LagrangeEllipticAssemble *assemble;
 
-    MatrixGraphs cellMatrixGraphs;
     MatrixGraphs solMatrixGraphs;
 
-    Vector fineInput;
-    Vector coarseInput;
     Vector fineSolution;
     Vector coarseSolution;
 
-    MonteCarloElliptic(int l,
-                       int dM,
-                       bool onlyFineLevel,
-                       Meshes *meshes,
-                       StochasticField *stochFields,
-                       LagrangeEllipticAssemble *assemble) :
-        MonteCarlo(l, dM, onlyFineLevel, meshes, stochFields),
+    Preconditioner *pc;
+    Solver *solver;
+    Newton *newton;
+
+    MonteCarloElliptic(int l, int dM, bool onlyFineLevel,
+                       Meshes *meshes, LagrangeEllipticAssemble *assemble) :
+        MonteCarlo(l, dM, onlyFineLevel, meshes),
         assemble(assemble),
-        cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
         solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)),
-        fineInput(Vector(cellMatrixGraphs[l - pLevel])),
-        coarseInput(Vector(cellMatrixGraphs[l - pLevel - 1])),
         fineSolution(Vector(solMatrixGraphs[l - pLevel])),
         coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) {
 
-        Initialize();
+        pc = GetPC("SuperLU");
+        solver = new Solver(pc, "GMRES");
+        newton = new Newton(*solver);
+    }
+
+    virtual ~MonteCarloElliptic() {
+        delete pc;
+        delete solver;
+        delete newton;
     }
 };
 

mlmc/src/montecarlo/MonteCarloTransport.cpp

@@ -11,44 +11,35 @@ void MonteCarloTransport::Initialize() {
 
 void MonteCarloTransport::method() {
     Initialize();
+    SampleID id;
     for (int m = M; m < M + dM; m++) {
-        string sampleDirName = mkSampleDirName(m, true);
-        plotter->SetDirectory(sampleDirName);
-
-        stochasticField->GetFineSample(l, fineInput);
-        assemble->problem->SetSample(&fineInput);
-        SolvePDE(l, m, true, fineQ, fineCost, fineSolution);
-
-//        if (plotting)
-//            plotter->PlotVector("u", fineSolution,
-//                                1, l, "VertexData");
+        id = SampleID{l, m, false};
+        plotter->SetDirectory(id);
+        assemble->problem->DrawSample(id);
+        SolvePDE(id, fineQ, fineCost, fineSolution);
 
         if (onlyFineLevel)
             coarseQ = 0.0, coarseCost = 0.0;
 
         else {
-            sampleDirName = mkSampleDirName(m, false);
-            plotter->SetDirectory(sampleDirName);
-
-            stochasticField->GetCoarseSample(l, fineInput, coarseInput);
-            assemble->problem->SetSample(&coarseInput);
-            SolvePDE(l - 1, m, false,
+            id = SampleID{l, m, true};
+            plotter->SetDirectory(id);
+            assemble->problem->DrawSample(id);
+            SolvePDE(SampleID(),
                      coarseQ, coarseCost, coarseSolution);
-
-//            if (plotting)
-//                plotter->PlotVector("u", coarseSolution,
-//                                    1, l - 1, "VertexData");
         }
         updateSums(fineCost, fineQ - coarseQ, fineQ);
     }
 }
 
-void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel,
-                                   double &valueQ, double &cost, Vector &solution) {
+void MonteCarloTransport::solvePDE(SampleID id,
+                                   double &valueQ,
+                                   double &cost,
+                                   Vector &solution) {
     Preconditioner *pc = GetPC("PointBlockGaussSeidel");
     Solver solver(pc, "GMRES");
 
-    TimeSeries timeSeries = getTimeSeries(fineLevel);
+    TimeSeries timeSeries = getTimeSeries(!id.coarse);
     assemble->resetTAssemble();
 
     timeInt->operator()(assemble, timeSeries, solution);

mlmc/src/montecarlo/MonteCarloTransport.hpp

@@ -6,6 +6,7 @@
 #include "timestepping/TimeIntegrator.hpp"
 #include "matrixgraph/CellMatrixGraph.hpp"
 #include "TimeSeries.h"
+#include "dof/BasicDoFs.hpp"
 
 
 class MonteCarloTransport : public MonteCarlo {
@@ -14,12 +15,7 @@ protected:
 
     void method() override;
 
-    void solvePDE(int l,
-                  int m,
-                  bool fineLevel,
-                  double &valueQ,
-                  double &cost,
-                  Vector &solution) override;
+    void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) override;
 
 private:
     double startTime = 0.0;
@@ -30,9 +26,9 @@ private:
 public:
     DGTransportAssemble *assemble;
 
-    Preconditioner* pc;
-    Solver* solver;
-    TimeIntegrator* timeInt;
+    Preconditioner *pc;
+    Solver *solver;
+    TimeIntegrator *timeInt;
 
     MatrixGraphs cellMatrixGraphs;
     MatrixGraphs faceMatrixGraphs;
@@ -43,25 +39,12 @@ public:
     Vector fineSolution;
     Vector coarseSolution;
 
-    /*
-     * Ultimativ braucht die MC nur noch level dM und baseLevel
-     * außerdem einmal Assemble! Die Assemble kennt das Problem und das
-     * Problem zieht ein neues Sample mit SampleID
-     *
-     *
-     * Meshes Objekt sollte sowieso hier erzeugt werden!
-     *
-     * Brauche Funktion CreateSampleID
-     */
-
-
     MonteCarloTransport(int l, int dM, bool baseLevel,
-                        Meshes *meshes, StochasticField *stochFields,
-                        DGTransportAssemble *assemble) :
-        MonteCarlo(l, dM, baseLevel, meshes, stochFields),
+                        Meshes *meshes, DGTransportAssemble *assemble) :
+        MonteCarlo(l, dM, baseLevel, meshes),
         assemble(assemble),
-        cellMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
-        faceMatrixGraphs(MatrixGraphs(*meshes, dof("face", 3))),
+        cellMatrixGraphs(MatrixGraphs(*meshes, dof(new CellDoF(3)))),
+        faceMatrixGraphs(MatrixGraphs(*meshes, dof(new FaceDoF(3)))),
         solMatrixGraphs(*meshes, *assemble->disc),
         fineInput(Vector(faceMatrixGraphs[l - pLevel])),
         coarseInput(Vector(faceMatrixGraphs[l - pLevel - 1])),
@@ -69,7 +52,7 @@ public:
         coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])),
         pc(GetPC("PointBlockGaussSeidel")),
         solver(new Solver(pc, "GMRES")),
-        timeInt(new TimeIntegrator(*solver)){
+        timeInt(new TimeIntegrator(*solver)) {
 
         config.get("startTime", startTime);
         config.get("endTime", endTime);

mlmc/src/montecarlo/MultilevelMonteCarlo.cpp

@@ -15,13 +15,11 @@ void MultilevelMonteCarlo::initializeMapMonteCarlo() {
 
 MonteCarlo *MultilevelMonteCarlo::createMonteCarlo(int l, int dM, bool onlyFineLevel) {
     auto ellipticAssemble = dynamic_cast<LagrangeEllipticAssemble *>(assemble);
-    auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
+//    auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
     if (ellipticAssemble != nullptr)
-        return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes,
-                                      stochasticField, ellipticAssemble);
-    if (transportAssemble != nullptr)
-        return new MonteCarloTransport(l, dM, onlyFineLevel, meshes,
-                                       stochasticField, transportAssemble);
+        return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes, ellipticAssemble);
+//    if (transportAssemble != nullptr)
+//        return new MonteCarloTransport(l, dM, onlyFineLevel, meshes, transportAssemble);
     Exit("Wrong Assembly Class")
 }
 

mlmc/src/montecarlo/MultilevelMonteCarlo.hpp

@@ -3,8 +3,11 @@
 
 #include "montecarlo/MonteCarloElliptic.hpp"
 #include "montecarlo/MonteCarloTransport.hpp"
+#include "main/Utils.hpp"
 
 
+typedef std::map<int, MonteCarlo *> MapMonteCarlo;
+
 class MultilevelMonteCarlo {
 public:
     int verbose = -1;
@@ -12,20 +15,19 @@ public:
 
     int maxLevel, pLevel;
     vector<int> &initLevels, &initSampleAmount;
-    map<int, MonteCarlo *> mapMonteCarlo;
+    MapMonteCarlo mapMonteCarlo;
 
     Meshes *meshes;
-    StochasticField *stochasticField;
     IAssemble *assemble;
 
     vector<int> levels = {};
     vector<int> samples = {};
-    vector<double> avgY {};
-    vector<double> avgQ {};
-    vector<double> varY {};
-    vector<double> varQ {};
-    vector<double> kurtosis {};
-    vector<int> avgCost {};
+    vector<double> avgY{};
+    vector<double> avgQ{};
+    vector<double> varY{};
+    vector<double> varQ{};
+    vector<double> kurtosis{};
+    vector<int> avgCost{};
     double value = 0.0, cost = 0.0;
     double epsilon = 0.0;
 
@@ -39,12 +41,10 @@ public:
     MultilevelMonteCarlo(vector<int> &initLevels,
                          vector<int> &initSampleAmount,
                          Meshes *meshes,
-                         StochasticField *stochasticField,
                          IAssemble *assemble) :
         initLevels(initLevels),
         initSampleAmount(initSampleAmount),
         meshes(meshes),
-        stochasticField(stochasticField),
         assemble(assemble) {
 
         config.get("MLMCPlotting", plotting);

mlmc/src/montecarlo/SampleID.hpp

@@ -5,14 +5,19 @@
 
 
 struct SampleID {
-    int level;
-    int m;
-    bool coarse;
+    int level = -1;
+    int m = -1;
+    bool coarse = false;
 
-    std::string ID() const {
+    SampleID() = default;
+
+    SampleID(int level, int m, int coarse) :
+        level(level), m(m), coarse(coarse) {}
+
+    std::string Str() const {
         return "l:" + std::to_string(level) +
-            "m:" + std::to_string(m) +
-            "c:" + std::to_string((int) coarse);
+            "-m:" + std::to_string(m) +
+            "-c:" + std::to_string((int) coarse);
     }
 };