Merge branch 'feature' into 'master'

Feature

See merge request mpp/mlmc!27

.gitlab-ci.yml

@@ -53,24 +53,44 @@ mpitest-mlmc:
   tags: [ docker ]
 
 
-elliptic-exercises-mlmc:
+elliptic-experiments-mlmc:
   stage: experiments
   variables:
     GIT_STRATEGY: none
     RUN_EXPERIMENTS: ''
   image: ${REGISTRY}/${IMAGE_NAME_MLMC}
   timeout: '1h'
-  only:
-    variables: [ $RUN_EXPERIMENTS == 'all', $RUN_EXPERIMENTS == 'elliptic' ]
+  #only:
+  #  variables: [ $RUN_EXPERIMENTS == 'all', $RUN_EXPERIMENTS == 'elliptic' ]
   script:
     - cd /mpp/notebooks
     - ls -al
     - jupyter nbconvert --ExecutePreprocessor.timeout=2400
       --execute --to html Elliptic\ Experiments.ipynb
+    - cp Elliptic\ Experiments.html $CI_PROJECT_DIR/
   dependencies: [ "build-mlmc" ]
   tags: [ docker ]
   artifacts:
-    paths: [ notebooks/Elliptic\ Experiments.html ]
+    paths: [ Elliptic\ Experiments.html ]
+
+
+elliptic-experiments-mlmc2:
+  stage: experiments
+  variables:
+    GIT_STRATEGY: none
+    RUN_EXPERIMENTS: ''
+  image: ${REGISTRY}/${IMAGE_NAME_MLMC}
+  timeout: '1h'
+  script:
+    - cd /mpp/notebooks
+    - ls -al
+    - jupyter nbconvert --ExecutePreprocessor.timeout=2400
+      --execute --to html Elliptic\ Experiments\ 2.ipynb
+    - cp Elliptic\ Experiments\ 2.html $CI_PROJECT_DIR/
+  dependencies: [ "build-mlmc" ]
+  tags: [ docker ]
+  artifacts:
+    paths: [ Elliptic\ Experiments\ 2.html ]
 
 
 deploy-mlmc:

CMakeLists.txt

@@ -2,9 +2,11 @@ cmake_minimum_required(VERSION 3.5.1)
 project(MLMC)
 
 set(SPACE_DIM 2 CACHE STRING "SPACE_DIM")
+set(USE_FFTW ON CACHE STRING "USE_FFTW")
+set(USE_SPLITTED_COMMS ON CACHE STRING "USE_SPLITTED_COMMS")
 set(USE_SPACETIME OFF CACHE STRING "USE_SPACETIME")
-set(COMPILER_OPTIMIZE -O0 CACHE STRING "COMPILER_OPTIMIZE")
-#set(COMPILER_OPTIMIZE -O3 CACHE STRING "COMPILER_OPTIMIZE")
+#set(COMPILER_OPTIMIZE -O0 CACHE STRING "COMPILER_OPTIMIZE")
+set(COMPILER_OPTIMIZE -O3 CACHE STRING "COMPILER_OPTIMIZE")
 set(NO_DEPRECATED OFF CACHE STRING "NO_DEPRECATED")
 set(AFFINE_LINEAR_TRAFO ON CACHE STRING "AFFINE_LINEAR_TRAFO")
 set(BUILD_TESTS OFF CACHE STRING "BUILD_TESTS")

mlmc/conf/mlmc_elliptic.conf

 # ----- Problem Settings -----
-Experiment = MLMCExperiment
-#Experiment = ConvergenceTest
-
 #Problem = StochasticLaplace1D
 Problem = StochasticLaplace2D
-
 Model = LagrangeElliptic
 #Model = MixedElliptic
 #Model = HybridElliptic
 #Model = DGElliptic
-
 degree = 1
-plevel = 2
-
+#degree = 2
 Quantity = L2
 #Functional = Energy
 #Funcitonal = H1
@@ -21,44 +15,30 @@ Quantity = L2
 
 # ----- Multilevel Monte Carlo -----
 maxLevel = 7
+MCParallel = false
 epsilon = 0.01
 onlyFine = false
 initLevels = [3, 4, 5]
 initSampleAmount = [12, 6, 3]
-uniformSampleAmount = 500
 
 # ----- Stochastic Field -----
 StochasticField = LogNormal
-mean = 0.0
+Mean = 0.0
 sigma = 1.0
-norm_p = 2
 lambda = [0.15, 0.15]
-smoothing = 1.8
-evtol = 1e-10
-
-# ----- Solver -----
-LinearReduction = 1e-12
-LinearEpsilon = 1e-10
-LinearSteps = 3000
-NewtonSteps = 1
-NewtonLineSearchSteps = 0
+smoothing = 1.0
 
 # ----- Plotting -----
-GeneratorPlotting = 0
-MCPlotting = 0
+PDESolverPlotting = 0
 
 # ----- Verbose -----
 MCVerbose = 1
-PDEVerbose = 1
+MLMCVerbose = 1
 MainVerbose = 1
 MeshVerbose = 1
-MLMCVerbose = 1
 ConfigVerbose = 1
-LinearVerbose = 1
-NewtonVerbose = 1
+LinearVerbose = 0
+NewtonVerbose = 0
+AssembleVerbose = 0
+PDESolverVerbose = 0
 GeneratorVerbose = 0
-
-# ----- Logging -----
-TimeLevel = -1
-MuteLevel = -1
-DebugLevel = -1
\ No newline at end of file

mlmc/conf/mlmc_transport.conf

 # ----- Problem Settings -----
-Experiment = MLMCExperiment
-#Experiment = ConvergenceTest
-
 #Problem = StochasticPollution1D
 #Problem = StochasticPollutionCosHat1D
 #Problem = StochasticPollution2D
 #Problem = StochasticPollutionCosHat2D
-Problem = StochasticPollutionMollifiedBar2D
-
+Problem = StochasticGaussHat2D
 Model = DGTransport
 #Model = PGTransport
+#degree = 0
+#degree = 1
+degree = 2
+Functional = Mass
+#Funcitonal = Energy
+#Functional = Outflow
 
 Overlap = dG1
 flux_alpha = 1  #Upwind: 1, Central: 0
 # Streamline diffusion parameter
 
-degree = 2
-plevel = 2
 
-Functional = Mass
-#Funcitonal = Energy
-#Functional = Outflow
+
 
 # ----- Time Series -----
 useDGTimeIntegrator = 0;
@@ -36,45 +34,30 @@ Keps = 1e-5
 
 # ----- Multilevel Monte Carlo -----
 maxLevel = 7
+MCParallel = false
 epsilon = 0.01
-mcOnly = false
-initLevels = [4, 5, 6]
-initSampleAmount = [8, 4, 2]
-uniformSampleAmount = 100
+onlyFine = false
+initLevels = [3, 4, 5]
+initSampleAmount = [12, 6, 3]
 
 # ----- Stochastic Field -----
 StochasticField = LogNormal
-mean = 1.0
+Mean = 0.0
 sigma = 1.0
-norm_p = 2
-lambda1 = 0.10
-lambda2 = 0.10
-smoothing = 1.9
-evtol = 1e-10
-
-# ----- Solver -----
-LinearReduction = 1e-12
-LinearEpsilon = 1e-10
-LinearSteps = 3000
-NewtonSteps = 1
-NewtonLineSearchSteps = 0
+lambda = [0.15, 0.15]
+smoothing = 1.0
 
 # ----- Plotting -----
-GeneratorPlotting = 1
-MCPlotting = 1
+PDESolverPlotting = 1
 
 # ----- Verbose -----
 MCVerbose = 1
-PDEVerbose = 1
-MainVerbose = 1
 MLMCVerbose = 1
-ConfigVerbose = 0
-LinearVerbose = -1
-NewtonVerbose = 1
-GeneratorVerbose = 1
-TimeIntegratorVerbose = 1
-
-# ----- Logging -----
-TimeLevel = -1
-MuteLevel = -1
-DebugLevel = -1
\ No newline at end of file
+MainVerbose = 1
+MeshVerbose = 1
+ConfigVerbose = 1
+LinearVerbose = 0
+NewtonVerbose = 0
+AssembleVerbose = 0
+PDESolverVerbose = 0
+GeneratorVerbose = 0

mlmc/conf/mlmc_vs_sg.conf

@@ -47,8 +47,7 @@ NewtonSteps = 1
 NewtonLineSearchSteps = 0
 
 # ----- Plotting -----
-GeneratorPlotting = 1
-MCPlotting = 1
+PDESolverPlotting = 1
 
 # ----- Verbose -----
 MCVerbose = 1

mlmc/src/CMakeLists.txt

@@ -2,6 +2,7 @@ include_directories(basics)
 add_subdirectory(basics)
 
 include_directories(generators)
+include_directories(generators/algorithms)
 add_subdirectory(generators)
 
 include_directories(problems)

mlmc/src/Main.cpp

 #include "m++.hpp"
 #include "Main.hpp"
 
+
 #ifdef COMPLEX
 #error undef COMPLEX in src/Compiler.h
 #endif
 
 int main(int argc, char **argv) {
-    Config::setSearchPath("../mlmc/");
-    Config::setConfigFileName("m++.conf");
-    Mpp::initialize(&argc, argv);
+  Config::setSearchPath("../mlmc/");
+  Config::setConfigFileName("m++.conf");
+  Mpp::initialize(&argc, argv);
+
+  if(!PowerOfTwo(PPM->Size(0))) {
+    Warning("Number of process is not a power of 2."
+            "This might lead to undefined behaviour.")
+  }
 
-    MainProgram mainProgram;
-    return mainProgram.Run();
+  MainProgram mainProgram;
+  return mainProgram.Run();
 }

mlmc/src/Main.hpp

@@ -4,30 +4,37 @@
 #include "montecarlo/MultilevelMonteCarlo.hpp"
 
 
+
+bool PowerOfTwo(int n) {
+  if(n==0) return false;
+  return (ceil(log2(n)) == floor(log2(n)));
+}
+
 class MainProgram {
 private:
-    int verbose = 1;
+  int verbose = 1;
 
 public:
-    MultilevelMonteCarlo mlmc;
-
-    MainProgram() : mlmc(MultilevelMonteCarlo()){
-        config.get("MainVerbose", verbose);
-    }
-
-    int Run() {
-        mout.StartBlock("MLMC Experiment");
-        mout << "Run" << endl;
-        mlmc.Method();
-        mout.EndBlock();
-        mout << endl;
-
-        mlmc.PrintMCResults();
-        mlmc.PrintExponentResults();
-        mlmc.PrintMLMCResults();
-
-        return 0;
-    }
+  MultilevelMonteCarlo mlmc;
+
+  MainProgram() : mlmc(MultilevelMonteCarlo()) {
+    config.get("MainVerbose", verbose);
+  }
+
+  int Run() {
+    config.PrintInfo();
+    mout.StartBlock("MLMC Experiment");
+    mout << "Run" << endl;
+    mlmc.Method();
+    mout.EndBlock();
+    mout << endl;
+
+    mlmc.PrintMCResults();
+    mlmc.PrintExponentResults();
+    mlmc.PrintMLMCResults();
+
+    return 0;
+  }
 };
 
 #endif //MAIN_HPP

mlmc/src/basics/CMakeLists.txt

 add_library(BASICS STATIC
         Level.cpp
-        Sample.cpp
-        PlotMap.cpp
         Utilities.cpp
         )
-target_link_libraries(BASICS fftw3 MPP_LIBRARIES)
\ No newline at end of file
+target_link_libraries(BASICS MPP_LIBRARIES)
\ No newline at end of file

mlmc/src/basics/Level.cpp

 #include "Level.hpp"
 
-
-void Level::Update(int level) {
-    fine = level;
-    coarse = level - 1;
-    if (coarse < 0) coarse = 0;
-    pLevel = coarse;
-}
-
-bool Level::operator==(const Level &level) const {
-    if (this->fine != level.fine)
-        return false;
-    if (this->coarse != level.coarse)
-        return false;
-    if (this->pLevel != level.pLevel)
-        return false;
-    return true;
-}
-
-bool Level::operator==(int level) const {
-    if (this->fine != level)
-        return false;
-    return true;
-}
-
-bool Level::operator!=(const Level &level) const {
-    return !(*this == level);
-}
-
-bool Level::operator!=(int level) const {
-    return !(*this == level);
-}
-
-bool Level::operator<=(const Level &level) const {
-    if (this->fine < level.fine) return true;
-    else return false;
-}
-
-bool Level::operator>=(const Level &level) const {
-    if (this->fine > level.fine) return true;
-    else return false;
-}
-
-bool Level::operator<(const Level &level) const {
-    if (this->fine < level.fine) return true;
-    else return false;
-}
-
-bool Level::operator>(const Level &level) const {
-    if (this->fine > level.fine) return true;
-    else return false;
-}
-
-bool Level::operator<(int level) const {
-    if (this->fine < level) return true;
-    else return false;
-}
-
-bool Level::operator>(int level) const {
-    if (this->fine > level) return true;
-    else return false;
-}
-
-bool Level::operator<=(int level) const {
-    if (this->fine <= level) return true;
-    else return false;
-}
-
-bool Level::operator>=(int level) const {
-    if (this->fine >= level) return true;
-    else return false;
-}

mlmc/src/basics/Level.hpp

 #ifndef LEVEL_HPP
 #define LEVEL_HPP
 
+#include "Logging.hpp"
 #include <vector>
 #include <map>
 
 
-/*
- * Todo
- * Only Fine true: pLevel = fine = coarse
- * Only Fine false: plevel = coarse, fine = coarse + 1
- *
- * Add commSplit to Level
- */
-
-struct Level {
-    int fine;
-    int coarse;
-    int pLevel;
-
-    Level() {}
-
-    explicit Level(int level) {
-        Update(level);
-    }
-
-    void Update(int level);
-
-    int ActualLevel(bool getCoarse) const {
-        if (!getCoarse) return fine;
-        else return coarse;
-    };
-
-    bool operator==(const Level &level) const;
-
-    bool operator==(int level) const;
-
-    bool operator!=(const Level &level) const;
-
-    bool operator!=(int level) const;
-
-    bool operator<=(const Level &level) const;
-
-    bool operator>=(const Level &level) const;
-
-    bool operator<(const Level &level) const;
-
-    bool operator>(const Level &level) const;
-
-    bool operator<(int level) const;
-
-    bool operator>(int level) const;
-
-    bool operator<=(int level) const;
-
-    bool operator>=(int level) const;
-};
-
 /*
  * Todo:
- * * make second key template for degree
+ * * idea: make second key template for degree
  *   to realize MultidegreeMonteCarlo
- * * also key (thus Level) should include proc information
+ * && ?
  */
 
 template<typename T>
 struct LevelMap {
 protected:
-    std::map<Level, T> _levelMap;
+    std::map<int, T> _levelMap;
 
 public:
     LevelMap() {};
 
-    LevelMap(std::initializer_list<std::pair<Level, T>> levelMap) {
-        for (std::pair<Level, T> pair : levelMap) {
+    LevelMap(std::initializer_list<std::pair<int, T>> levelMap) {
+        for (std::pair<int, T> pair : levelMap) {
             _levelMap.emplace(pair.first, T(pair.second));
         }
     }
@@ -80,7 +30,7 @@ public:
     std::vector<int> GetLevelVector() const {
         std::vector<int> levelVector;
         for (auto &&[level, entry]: _levelMap)
-            levelVector.push_back(level.fine);
+            levelVector.push_back(level);
         return levelVector;
     }
 
@@ -123,7 +73,7 @@ public:
 
     auto clear() { _levelMap.clear(); }
 
-    auto insert(std::pair<Level, T> pair) { _levelMap.insert(pair); }
+    auto insert(std::pair<int, T> pair) { _levelMap.insert(pair); }
 
     auto begin() { return _levelMap.begin(); }
 
@@ -141,11 +91,20 @@ public:
 
     auto rend() const { return _levelMap.end(); }
 
-    auto operator[](Level level) { return _levelMap[level]; }
+    auto operator[](int level) { return _levelMap[level]; }
 
-    auto at(Level level) { return _levelMap.at(level); }
+    auto at(int level) { return _levelMap.at(level); }
 
-    auto find(Level level) { return _levelMap.find(level); }
+    auto find(int level) { return _levelMap.find(level); }
+
+    friend Logging &operator<<(Logging &s, const LevelMap<T> &levelMap) {
+        s << "[";
+        for(auto &[level, entry] : levelMap) {
+            if(level != levelMap.rbegin()->first) s << entry << ", ";
+            else s << entry;
+        }
+        return s << "]";
+    }
 };
 
 #endif //LEVEL_HPP

mlmc/src/basics/PlotMap.cpp

-#include "PlotMap.hpp"
-
-// Todo bring implementations over here!
\ No newline at end of file

mlmc/src/basics/PlotMap.hpp

-#ifndef PLOTMAP_HPP
-#define PLOTMAP_HPP
-
-#include "Sample.hpp"
-
-#include "Plot.hpp"
-#include "dof/LagrangeDoF.hpp"
-
-#include <sys/stat.h>
-
-
-/*
- * Todo make singelton
- */
-
-struct PlotMap : public  std::map<int, Plot *> {
-    const std::string defaultDir = "data/vtk/";
-
-    std::string currentDir = defaultDir;
-
-    std::string currentDirShort = "";
-
-    std::string currentFileShort = "";
-
-    std::string currentFile = "";
-
-    PlotMap() {};
-
-    ~PlotMap() {
-        for (auto &[level, plot]:*this)
-            delete plot;
-    }
-
-    PlotMap(std::allocator<std::pair<int, Plot *>> plotMap) :
-        std::map<int, Plot *>(plotMap) {};
-
-    void UpdateMap(SampleSolution &dataToPlot) {
-        int level = GetLevelInteger(dataToPlot);
-        UpdateMap(level, dataToPlot.U.GetMesh());
-    }
-
-    void UpdateMap(int newLevel, const Mesh &mesh) {
-        if (this->find(newLevel) == this->end())
-            (*this)[newLevel] = new Plot(mesh);
-        else
-            return;
-    }
-
-    int GetLevelInteger(const SampleSolution &dataToPlot) {
-        int level;
-        if (dataToPlot.id.coarse)
-            level = dataToPlot.id.level.coarse;
-        else
-            level = dataToPlot.id.level.fine;
-        return level;
-    }
-
-    string checkData(const SampleSolution &dataToPlot) {
-        return dataToPlot.U.Name();
-    }
-
-    void VtkPlot(const SampleSolution &dataToPlot) {
-        string dataType = checkData(dataToPlot);
-        int level = GetLevelInteger(dataToPlot);
-        int dataDim = dataToPlot.U.sizeDoF();
-        UpdateDirectory(dataToPlot);
-        UpdateMap(level, dataToPlot.U.GetMesh());