Commit bb84c0fa authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch 'feature' into 'master'

Feature

See merge request mpp/mlmc!13
parents bd34db10 5247a17f
Pipeline #141950 passed with stages
in 9 minutes and 28 seconds
......@@ -4,11 +4,12 @@ stages:
- experiments
- deploy
variables:
OS: "ubuntu"
OS_VERSION: "18.04"
OS_VERSION_MLMC: "20.04"
REGISTRY: "ci.cluster.math.kit.edu"
IMAGE_NAME_MLMC: "mlmc-${CI_COMMIT_SHORT_SHA}-${OS}${OS_VERSION}"
IMAGE_NAME_MLMC: "mlmc-${CI_COMMIT_SHORT_SHA}-${OS}${OS_VERSION_MLMC}-${CI_PIPELINE_ID}"
build-mlmc:
......@@ -18,7 +19,7 @@ build-mlmc:
script:
- sed s/REGISTRY/${REGISTRY}/g docker/mlmc.baseimage >
docker/mlmc_registry.baseimage
- sed s/UBUNTUVERSION/${OS_VERSION}/g docker/mlmc_registry.baseimage >
- sed s/UBUNTUVERSION/${OS_VERSION_MLMC}/g docker/mlmc_registry.baseimage >
docker/${IMAGE_NAME_MLMC}.baseimage
- docker build
--build-arg UPSTREAM_COMMIT=${UPSTREAM_COMMIT}
......@@ -28,16 +29,26 @@ build-mlmc:
tags: [ shell ]
unit-test-mlmc:
ctest-mlmc:
stage: test
variables:
GIT_STRATEGY: none
image: ${REGISTRY}/${IMAGE_NAME_MLMC}
script:
- cd /mpp/build/
- ctest
dependencies: [ "build-mlmc" ]
tags: [ docker ]
mpitest-mlmc:
stage: test
variables:
GIT_STRATEGY: none
image: ${REGISTRY}/${IMAGE_NAME_MLMC}
script:
- cd /mpp/build/tests
- mpirun -n 4 TestMainProgramElliptic
# - mpirun -n 4 TestMainProgramTransport
# - mpirun -n 4 TestMultilevelPlotter
- cd /mpp/build/
- python3 mppyrun.py --mpi_tests=1 --mute=0
dependencies: [ "build-mlmc" ]
tags: [ docker ]
......@@ -68,10 +79,10 @@ deploy-mlmc:
GIT_STRATEGY: none
script:
- docker tag ${REGISTRY}/${IMAGE_NAME_MLMC}
${REGISTRY}/mlmc-${CI_COMMIT_TAG}-${OS}${OS_VERSION}
- docker push ${REGISTRY}/mlmc-${CI_COMMIT_TAG}-${OS}${OS_VERSION}
${REGISTRY}/mlmc-${CI_COMMIT_TAG}-${OS}${OS_VERSION_MLMC}
- docker push ${REGISTRY}/mlmc-${CI_COMMIT_TAG}-${OS}${OS_VERSION_MLMC}
only: [ tags ]
dependencies: [ "unit-test-mlmc" ]
dependencies: [ "ctest-mlmc", "mpitest-mlmc" ]
tags: [ shell ]
......
cmake_minimum_required(VERSION 3.5.1)
project(MLMC)
# Compiler options
#set(CXX_VERSION "17")
#set(COMPILER_VERSION "c++")
#set(COMPILER_VERSION "gnu++")
set(COMPILER_OPTIMIZE "-O0")
#set(COMPILER_OPTIMIZE "-O1")
#set(COMPILER_OPTIMIZE "-O2")
#set(COMPILER_OPTIMIZE "-O3")
#set(COMPILER_OPTIMIZE "-Ofast")
# Problem options
option(NO_DEPRECATED "Suppress deprecated warnings" ON)
option(PROBLEM_NO_TIME "Time independent problem" ON)
option(PROBLEM_1D "1 dimensional problem" OFF)
option(PROBLEM_2D "2 dimensional problem" ON)
option(AFFINE_LINEAR_TRAFO "Only affine linear transformations" ON)
# Constants
#set(CONST_GEOMETRIC_TOLERANCE "1e-10")
#set(CONST_NEAR_ZERO "1e-15")
#set(CONST_VERY_LARGE "1e30")
#set(CONST_INFTY "1e100")
# SuperLU options
option(USE_SUPERLU30 "Use SuperLU 3.0 (old version)" OFF)
# Test options
option(BUILD_TESTS "Build test executables" OFF)
option(BUILD_MLMC_TESTS "Build mlmc test executables" ON)
set(SPACE_DIM 2 CACHE STRING "SPACE_DIM")
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(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")
set(BUILD_MLMC_TESTS ON CACHE STRING "BUILD_MLMC_TESTS")
#---------------------------------------------------------------------------------------#
# Set path for mpp
set(PROJECT_MPP_DIR ${PROJECT_SOURCE_DIR}/mpp)
# include M++
include(mpp/CMakeLists.txt)
# Include settings for mpp
include(${PROJECT_MPP_DIR}/CMakeListsMpp.inc)
# Manage folder structure build folder
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/py)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/vtk)
file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/data/vtu)
# Link libraries
#find_package(FFTW3 REQUIRED) (Can't be found but still seems to link in)
link_directories(${PROJECT_SOURCE_DIR}/sprng5/lib)
# Include directories
include_directories(${PROJECT_SOURCE_DIR}/sprng5/include)
include_directories(${PROJECT_SOURCE_DIR}/mlmc/src)
# Subdirectory
# Subdirectories
add_subdirectory(mlmc/src)
# Libraries
set(MLMC_LIBRARIES MLMC sprng fftw3 m ${MPP_LIBRARIES})
# Executables
add_executable(MLMC-M++ mlmc/src/Main.cpp)
# Linking
target_link_libraries(MLMC-M++ ${MLMC_LIBRARIES})
target_link_libraries(MLMC-M++ MONTECARLO)
# Tests
if (BUILD_MLMC_TESTS)
add_subdirectory(tests)
include_directories(mlmc/tests/)
add_subdirectory(mlmc/tests/)
endif ()
......@@ -5,15 +5,15 @@ Experiment = MLMCExperiment
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = DGFEM
Model = LagrangeElliptic
#Model = MixedElliptic
#Model = HybridElliptic
#Model = DGElliptic
degree = 1
plevel = 2
Functional = L2
Quantity = L2
#Functional = Energy
#Funcitonal = H1
#Functional = Outflow
......@@ -56,7 +56,7 @@ MLMCVerbose = 1
ConfigVerbose = 1
LinearVerbose = 1
NewtonVerbose = 1
GeneratorVerbose = 1
GeneratorVerbose = 0
# ----- Logging -----
TimeLevel = -1
......
include(assemble/CMakeLists.txt)
include(main/CMakeLists.txt)
include(montecarlo/CMakeLists.txt)
include(problem/CMakeLists.txt)
include(stochastics/CMakeLists.txt)
add_library(MLMC STATIC ${assemble} ${main}
${montecarlo} ${problem} ${stochastics})
include_directories(basics)
add_subdirectory(basics)
include_directories(generators)
add_subdirectory(generators)
include_directories(problems)
add_subdirectory(problems)
include_directories(pdesolver)
add_subdirectory(pdesolver)
include_directories(montecarlo)
add_subdirectory(montecarlo)
#include "m++.hpp"
#include "main/MainProgram.hpp"
#include "Main.hpp"
#ifdef COMPLEX
#error undef COMPLEX in src/Compiler.h
......@@ -10,9 +10,6 @@ int main(int argc, char **argv) {
Config::setConfigFileName("m++.conf");
Mpp::initialize(&argc, argv);
std::unique_ptr<MainProgram> mainProgram = std::make_unique<MainProgram>();
mainProgram->Initialize();
int rc = mainProgram->Run();
return rc;
MainProgram mainProgram;
return mainProgram.Run();
}
#ifndef MAIN_HPP
#define MAIN_HPP
#include "montecarlo/MultilevelMonteCarlo.hpp"
class MainProgram {
private:
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;
}
};
#endif //MAIN_HPP
set(assemble
assemble/LagrangeEllipticAssemble.cpp
assemble/MixedEllipticAssemble.cpp
assemble/HybridEllipticAssemble.cpp
assemble/DGEllipticAssemble.cpp
assemble/DGTransportAssemble.cpp
assemble/PGReactionAssemble.cpp
assemble/DGReactionAssemble.cpp
)
\ No newline at end of file
#ifndef TUTORIAL_DGTRANSPORT_HPP
#define TUTORIAL_DGTRANSPORT_HPP
#include "timestepping/DGTAssemble.hpp"
#include "assemble/IStochasticTransportAssemble.hpp"
#include "problem/StochasticTransportProblem.hpp"
#include "TimeSeries.hpp"
#include "utility/Logging.hpp"
#include "elements/DGFaceElement.hpp"
class DGTransportAssemble : public IStochasticTransportAssemble {
protected:
DGDiscretization *disc;
double flux_alpha = 1.0;
double diffusion = 0.0;
public:
DGTransportAssemble(IDiscretization *disc, IStochasticTransportProblem *problem) :
IStochasticTransportAssemble(problem),
disc(dynamic_cast<DGDiscretization *>(disc)) {
config.get("flux_alpha", flux_alpha);
config.get("diffusion", diffusion);
}
IDiscretization *GetDisc() override {
return disc;
};
void PrintInfo() const {
mout.PrintInfo("Assemble", 1,
PrintInfoEntry("Name", Name()),
// PrintInfoEntry("Problem", problem->Name()),
PrintInfoEntry("Discretization", disc->DiscName()));
}
void PrintInfo(const Vector &u) const override {
// std::pair<double, double> rate = InFlowOutFlowRate(u);
// mout << "n=" << std::setw(3) << std::left << step
// << " t=" << std::setw(7) << std::left << t_
// << " Energy=" << std::setw(13) << std::left << Energy(u)
// << " Mass=" << std::setw(13) << std::left << Mass(u);
// if (abs(rate.first) >= 10e-10)
// mout << " IFR=" << std::setw(13) << std::left << rate.first;
// if (abs(rate.second) >= 10e-10)
// mout << " OFR=" << std::setw(13) << std::left << rate.second;
// mout << endl;
}
const char *Name() const override;
double Residual(const Vector &u, Vector &b) const override;
void Jacobi(const Vector &u, Matrix &A) const override;
void MassMatrix(Matrix &massMatrix) const override;
void SystemMatrix(Matrix &systemMatrix) const override;
void RHS(double t, Vector &rhs) const override;
double Energy(const Vector &u) const override;
void Energy(const cell &c, const Vector &u, double &energy) const override;
double Error(double t, const Vector &u) const;
std::pair<double, double> InFlowOutFlowRate(const Vector &u) const;
void SetInitialValue(Vector &u) const {
u = 0;
Point z[MaxNodalPoints];
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(*disc, u, c);
Elem.NodalPoints(z);
double lambda = 1e-9;
for (int j = 0; j < Elem.dimension(); ++j)
u(r, j) = problem->Solution(0, lambda * c() + (1 - lambda) * z[j]);
}
Accumulate(u);
}
void VtkPlotting(double t, const Vector &u) const {
// const Mesh &Mp = plot->GetMesh();
// const Mesh &Mu = u.GetMesh();
// if (Mp.CellCount() != Mu.CellCount()) return;
// vout (2) << " => VtkPlotting result of step " << step << ".\n";
// char filename[128];
// VtkPlotting_cell(t, u, filename);
}
void PrintMatrixInfo(Matrix &A, int diagonal) const {
for (cell c = A.cells(); c != A.cells_end(); ++c) {
DGElement elem(*disc, A, c);
DGRowEntries A_c(A, c, c);
for (int i = 0; i < elem.dimension(); ++i) {
for (int j = 0; j < elem.dimension(); ++j)
mout << A_c(i, j) << " ";
mout << endl;
}
mout << endl;
if (diagonal) continue;
for (int f = 0; f < c.Faces(); ++f) {
cell cf = c;
if (A.GetMesh().onBndDG(c, f))
continue;
else
cf = A.GetMesh().find_neighbour_cell(c, f);
DGRowEntries A_cf(A, c, cf);
for (int i = 0; i < elem.dimension(); ++i) {
for (int j = 0; j < elem.dimension(); ++j)
mout << A_cf(i, j) << " ";
mout << endl;
}
mout << endl;
}
mout << "------------------------" << endl;
}
}
void SetExactSolution(double t, Vector &u_ex) const;
void AssembleTransfer(TransferMatrix &TM) const override;
};
#endif //TUTORIAL_DGTRANSPORT_HPP
#ifndef ISTOCHASTICTRANSPORTASSEMBLE_HPP
#define ISTOCHASTICTRANSPORTASSEMBLE_HPP
#include "Assemble.hpp"
#include "problem/StochasticTransportProblem.hpp"
class IStochasticTransportAssemble : public TAssemble {
protected:
IStochasticTransportProblem *problem;
public:
explicit IStochasticTransportAssemble(IStochasticTransportProblem *problem) :
problem(problem) {
};
void DrawSample(const SampleID &id) {
problem->DrawSample(id);
}
IStochasticProblem *GetProblem() {
return problem;
};
virtual IDiscretization *GetDisc() = 0;
virtual ~IStochasticTransportAssemble() = default;
};
#endif //ISTOCHASTICTRANSPORTASSEMBLE_HPP
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
#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;
}
#ifndef LEVEL_HPP
#define LEVEL_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
* to realize MultidegreeMonteCarlo
* * also key (thus Level) should include proc information
*/
template<typename T>
struct LevelMap {
protected:
std::map<Level, T> _levelMap;
public:
LevelMap() {};
LevelMap(std::initializer_list<std::pair<Level, T>> levelMap) {
for (std::pair<Level, T> pair : levelMap) {
_levelMap.emplace(pair.first, T(pair.second));
}
}
std::vector<int> GetLevelVector() const {
std::vector<int> levelVector;
for (auto &&[level, entry]: _levelMap)
levelVector.push_back(level.fine);
return levelVector;
}
std::vector<double> GetQVector() const {
std::vector<double> qVector;
for (auto &&[level, entry]: _levelMap)
qVector.push_back(entry.Q);
return qVector;
}
std::vector<double> GetYVector() const {
std::vector<double> yVector;
for (auto &&[level, entry]: _levelMap)
yVector.push_back(entry.Y);
return yVector;