Refactored logger and plotting for sample generator + bug fix hybrid flux generator

mlmc/conf/mlmc.conf

@@ -26,7 +26,7 @@ StochasticField = LogNormal
 Experiment = MLMCExperiment
 #Experiment = MLMCOverEpsilon
 
-initLevels = 3, 4
+initLevels = 2, 3
 initSampleAmount = 1, 0
 
 epsilon = 0.01
@@ -34,7 +34,7 @@ mcOnly = true
 epsilon_lst = 0.05, 0.01, 0.005
 
 maxLevel = 4
-plevel = 2
+plevel = 1
 degree = 2
 
 #functional = L2
@@ -44,9 +44,12 @@ functional = Energy
 #functional = Goal
 #functional = Mass
 
+GeneratorVerbose = 1
 MCVerbose = 2
 MLMCVerbose = 2
 PDEVerbose = 1
+
+GeneratorPlotting = 4
 MCPlotting = 3
 MLMCPlotting = 2
 
@@ -79,13 +82,12 @@ NewtonLineSearchSteps = 0
 ConfigVerbose = 0
 TimeLevel = -1
 LinearVerbose = -1
-BaseSolverVerbose = -1
 MuteLevel = -1
 NewtonVerbose = -1
 DebugLevel = -1
 
 #Transport
-scaling = 8
+scaling = 16
 TimeSeries = uniform
 startTime = 0.0
 endTime = 1.0
@@ -99,7 +101,4 @@ Keps = 1e-5
 rkorder = -2  #Implicit MP-rule
 #rkorder = 1000  #Polynomial Krylov-Arnoldi method
 
-gamma = 0.01
-
-
-plotPermForTransport = true
\ No newline at end of file
+gamma = 0.01
\ No newline at end of file

mlmc/src/MonteCarloLogger.h

@@ -120,4 +120,44 @@ public:
     }
 };
 
+class GeneratorLogger : public MonteCarloLogger {
+private:
+    Date timestamp;
+    int verbose = 0;
+
+public:
+    GeneratorLogger() : MonteCarloLogger() {
+        ReadConfig(Settings, "GeneratorVerbose", verbose);
+        indent += defaultIndent;
+        innerIndent += defaultIndent;
+    }
+
+    void StartMethodMsg() override {
+        Date start;
+        timestamp = start;
+        vout (1) << indent << "<Generating input sample>" << endl;
+    }
+
+    void EndMethodMsg() override {
+        vout (1) << indent << "<End after " << Date() - timestamp << ">" << endl;
+    }
+
+    void LogMsgv1(const string &msg) override {
+        vout(1) << innerIndent << msg << endl;
+    }
+
+    void LogMsgv2(const string &msg) override {
+        vout(2) << innerIndent << msg << endl;
+    }
+
+    void LogMsgv1Flush(const string &msg) {
+        vout(1) << "\r" << innerIndent << msg << flush;
+    }
+
+    void LogMsgv2Flush(const string &msg) {
+        vout(2) << "\r" << innerIndent << msg << flush;
+    }
+};
+
+
 #endif //MLMC__MULTILEVELMONTECARLO__HPP

mlmc/src/mc/MonteCarloTransport.C

@@ -12,8 +12,8 @@ void MonteCarloTransport::Method() {
                              + " dM: " + to_string(dM)
                              + " onlyFineLevel: True");
 
-    Vector fineFlux(fluxMatrixGraphs[l - pLevel]);
-    Vector coarseFlux(fluxMatrixGraphs[l - pLevel - 1]);
+    Vector fineFlux(faceMatrixGraphs[l - pLevel]);
+    Vector coarseFlux(faceMatrixGraphs[l - pLevel - 1]);
     Vector coarseFluxUp(fineFlux);
     Vector fineSolution(solMatrixGraphs[l - pLevel]);
     Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
@@ -28,16 +28,20 @@ void MonteCarloTransport::Method() {
     assemble->plot = &plot;
 
     for (int m = M; m < M + dM; m++) {
+        string sampleDirName = mkSampleDirName(m);
         mkSampleDir(m);
+
+        stochasticField->generators[l]->SetSampleDir(sampleDirName);
         stochasticField->GetFineSample(l, fineFlux);
-        assemble->sampleDir = mkSampleDirName(m);
+
+        assemble->sampleDir = sampleDirName;
         assemble->problem->LoadNewSample(&fineFlux);
 
-        if (plotting > 1) {
-            string name = buildName(m, mkSampleDirName(m).append("/flux"));
-            plot.celldata(fineFlux, 3);
-            plot.vtk_cellvector(name);
-        }
+        //        if (plotting > 1) {
+//            string name = buildName(m, mkSampleDirName(m).append("/flux"));
+//            plot.celldata(fineFlux, 3);
+//            plot.vtk_cellvector(name);
+//        }
 
         solveTransport(l, true, fineQ, fineCost, fineSolution);
 
@@ -73,12 +77,15 @@ void MonteCarloTransport::solveTransport(int l,
                                          double &cost,
                                          Vector &solution) {
     assemble->logger->StartMethodMsg();
-    Solver solver;
-    Preconditioner *BJ = GetPC("PointBlockJacobi");
-    Solver IM(BJ);
-    DGTimeIntegrator timeIntegrator(solver, solver, solver, IM);
+
+    Solver s; // TODO RMV
+    Preconditioner *pc = GetPC("PointBlockGaussSeidel"); // TODO ADD MG AND ASSEMBLE
+    Solver solver(pc, "GMRES");
+    DGTimeIntegrator timeIntegrator(s, s, s, solver);
+
     TimeSeries timeSeries = getTimeSeries(fineLevel);
     assemble->resetTAssemble();
+
     timeIntegrator(*assemble, timeSeries, solution);
     cost = solution.pSize() * assemble->getStep();
 

mlmc/src/mc/MonteCarloTransport.h

@@ -25,6 +25,7 @@ public:
     double endTime = 0.0;
 
     MatrixGraphs fluxMatrixGraphs;
+    MatrixGraphs faceMatrixGraphs;
     CellMatrixGraphs solMatrixGraphs;
 
     Transfer *transfer;
@@ -34,6 +35,7 @@ public:
         MonteCarlo(l, dM, baseLevel, meshes, stochFields),
         assemble(assemble),
         fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
+        faceMatrixGraphs(MatrixGraphs(*meshes, dof("face", 3))),
         solMatrixGraphs(*meshes, *assemble->disc),
         plot(Plot((*meshes)[l], meshes->dim())) {
 

mlmc/src/stochastics/CirculantEmbedding.C

 #include "CirculantEmbedding.h"
 
 
-void CirculantEmbedding1D::GetFineSample(Vector &fineField) {
+void CirculantEmbedding1D::GenerateFineSample(Vector &fineField) {
     if (internalCounter == 0)
         fineComplexField = generateField();
 
@@ -150,9 +150,9 @@ vector<double> CirculantEmbedding1D::getSqrtEV() {
     embedCovarianceMatrix(ar, ac, br);
     computeEV(br, ev);
 
-    for (int i = 0; i < ev.size(); i++) {
-        if (ev[i] > 0) continue;
-        else if (ev[i] < 0 && ev[i] > evtol) ev[i] = 0.0;
+    for (double & i : ev) {
+        if (i > 0) continue;
+        else if (i < 0 && i > evtol) i = 0.0;
         else padding(ar, ac);
     }
 
@@ -162,7 +162,7 @@ vector<double> CirculantEmbedding1D::getSqrtEV() {
     return sqrt_ev_return;
 }
 
-void CirculantEmbedding2D::GetFineSample(Vector &fineField) {
+void CirculantEmbedding2D::GenerateFineSample(Vector &fineField) {
     if (internalCounter == 0)
         fineComplexField = generateField();
 
@@ -284,7 +284,7 @@ void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar,
 void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar,
                                                  vector<vector<double>> &ac,
                                                  vector<vector<double>> &br) {
-    int br_size = 2 * ar.size() - 1;
+    int br_size = (int) (2 * ar.size() - 1);
     br.resize(br_size);
     for (int i = 0; i < br.size(); i++) {
         br[i].reserve(br_size);

mlmc/src/stochastics/CirculantEmbedding.h

@@ -33,7 +33,7 @@ public:
     DeterministicPermeabilityGenerator(int l, Meshes &meshes)
         : CirculantEmbedding(l, meshes) {}
 
-    void GetFineSample(Vector &fineField) override {
+    void GenerateFineSample(Vector &fineField) override {
         fineField = 1.0;
     };
 
@@ -59,7 +59,7 @@ public:
         numberOfCellsEmbedded = sqrtEigenvalues.size();
     }
 
-    void GetFineSample(Vector &fineField) override;
+    void GenerateFineSample(Vector &fineField) override;
 
     void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
 
@@ -92,7 +92,7 @@ public:
         numberOfCellsEmbedded[1] = sqrtEigenvalues[0].size();
     }
 
-    void GetFineSample(Vector &fineField) override;
+    void GenerateFineSample(Vector &fineField) override;
 
     void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
 

mlmc/src/stochastics/HybridFluxGenerator.C

 #include "HybridFluxGenerator.h"
 
 
-void HybridFluxGenerator::GetFineSample(Vector &fineField) {
-    MatrixGraphs cellMatrixGraph(meshes, dof("cell", 1));
-    MatrixGraphs faceMatrixGraph(meshes, dof("face", 1));
-    Vector perm(cellMatrixGraph[l - meshes.pLevel()]);
-    Vector solution(faceMatrixGraph[l - meshes.pLevel()]);
-    permeabilityGenerator->GetFineSample(perm);
+void HybridFluxGenerator::GenerateFineSample(Vector &fineField) {
+    MatrixGraphs cellMg1(meshes, dof("cell", 1));
+    MatrixGraphs cellMg3(meshes, dof("cell", 3));
+    MatrixGraphs faceMg(meshes, dof("face", 1));
 
-    bool plotPermForTransport = false;
-    config.get("plotPermForTransport", plotPermForTransport);
-    if (plotPermForTransport) {
-        Plot plot(meshes[l], 2);
-        plot.celldata(perm, 1);
-        plot.vtk_celldata("permeability");
-    }
+    Vector &normalFlux = fineField;
+    Vector flux(cellMg3[l - meshes.pLevel()]);
+    Vector perm(cellMg1[l - meshes.pLevel()]);
+    Vector u(faceMg[l - meshes.pLevel()]);
 
+    permeabilityGenerator->GetFineSample(perm);
     assemble->problem->LoadNewSample(&perm);
 
-    Preconditioner *pc = GetPC(faceMatrixGraph, *assemble, "SuperLU");
-
+    Preconditioner *pc = GetPC(faceMg, *assemble, "SuperLU");
     Solver solver(pc, "GMRES");
-
     Newton newton(solver);
 
-    newton(*assemble, solution);
-    assemble->setFlux(solution, fineField);
+    newton(*assemble, u);
+    assemble->setFlux(u, flux);
+    assemble->SetNormalFlux(u, normalFlux);
+
+    PlotSample(perm, u, flux, normalFlux);
 }
 
 void HybridFluxGenerator::GetCoarseSample(const Vector &fineField,
@@ -41,4 +38,26 @@ void HybridFluxGenerator::GetCoarseSample(const Vector &fineField,
         for (int i = 0; i < coarseRow.n(); ++i)
             coarseField(coarseRow, i) /= c.Children();
     }
-}
\ No newline at end of file
+}
+
+void HybridFluxGenerator::PlotSample(const Vector &perm,
+                                     const Vector &u,
+                                     const Vector &flux,
+                                     const Vector &normalFlux) {
+    if (plotting >= 4) {
+//        string nameU = sampleDir + "/u_hybrid";
+//        assemble->plotU(*plot, u, nameU);
+    }
+    if (plotting >= 3) {
+        string namePerm = sampleDir + "/perm";
+        assemble->plotPerm(*plot, perm, namePerm);
+    }
+    if (plotting >= 2) {
+        string nameFlux = sampleDir +"/flux";
+        assemble->plotFlux(*plot, flux, nameFlux);
+    }
+    if (plotting >= 1) {
+        string nameNormalFlux = sampleDir +"/normalFlux";
+        // Todo Normal Flux
+    }
+}

mlmc/src/stochastics/HybridFluxGenerator.h

@@ -12,10 +12,14 @@ protected:
     StochasticEllipticProblem *problem = nullptr;
     HybridEllipticAssemble *assemble = nullptr;
     CirculantEmbedding *permeabilityGenerator = nullptr;
+    Plot *plot = nullptr;
 
 public:
     explicit HybridFluxGenerator(int l, Meshes &meshes)
-        : SampleGenerator(l, meshes) {}
+        : SampleGenerator(l, meshes) {
+        if(plotting > 0)
+            plot = new Plot(meshes[l], meshes.dim());
+    }
 
     ~HybridFluxGenerator() {
         delete disc;
@@ -24,7 +28,12 @@ public:
         delete permeabilityGenerator;
     };
 
-    void GetFineSample(Vector &fineField) override;
+    void PlotSample(const Vector &perm,
+                    const Vector &u,
+                    const Vector &flux,
+                    const Vector &normalFlux);
+
+    void GenerateFineSample(Vector &fineField) override;
 
     void GetCoarseSample(const Vector &fineField, Vector &coarseField) override;
 };

mlmc/src/stochastics/SampleGenerator.h

 #ifndef SAMPLEGENERATOR_H
 #define SAMPLEGENERATOR_H
 
+#include <utility>
+
 #include "m++.h"
+#include "MonteCarloLogger.h"
 
 
 class SampleGenerator {
+protected:
+    int plotting = 0;
+    GeneratorLogger *logger;
+    string sampleDir = "";
+
+    virtual void GenerateFineSample(Vector &fineField) = 0;
+
 public:
     int l;
     Meshes &meshes;
 
-    explicit SampleGenerator(int l, Meshes &meshes) : l(l), meshes(meshes) {}
+    explicit SampleGenerator(int l, Meshes &meshes) : l(l), meshes(meshes) {
+        config.get("GeneratorPlotting", plotting);
+        logger = new GeneratorLogger();
+    }
 
-    virtual void GetFineSample(Vector &fineField) = 0;
+    ~SampleGenerator() {
+        delete logger;
+    }
+
+    virtual void GetFineSample(Vector &fineField) {
+        logger->StartMethodMsg();
+        GenerateFineSample(fineField);
+        logger->EndMethodMsg();
+    }
 
     virtual void GetCoarseSample(const Vector &fineField, Vector &coarseField) = 0;
+
+    void SetSampleDir(const string& sampleDirName) {
+        sampleDir = sampleDirName;
+    }
 };
 
 #endif //SAMPLEGENERATOR_H