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

major refactoring

parent 25ced1ba
Experiment = MLMCExperiment
#Experiment = ConvergenceTest
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = DGFEM
#sign = 1
#penalty = 20
#Problem = StochasticLaplace1D
#Problem = DeterministicLaplace1D
Problem = StochasticLaplace2D
degree = 1
plevel = 2
functional = L2
#functional = Energy
......@@ -16,34 +18,16 @@ functional = L2
#functional = Outflow
#functional = Goal
StochasticField = LogNormal
#StochasticField = Gaussian
#Experiment = ConvergenceTest
Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
# Multilevel Monte Carlo
maxLevel = 9
epsilon = 0.01
mcOnly = false
initLevels = 3, 4, 5
initSampleAmount = 12, 6, 3
uniformSampleAmount = 500
epsilon = 0.005
mcOnly = false
maxLevel = 9
plevel = 2
degree = 1
MainVerbose = 1
MCVerbose = 2
MLMCVerbose = 2
PDEVerbose = 0
GeneratorVerbose = 0
GeneratorPlotting = 1
MCPlotting = 0
# Stochastic Field
StochasticField = LogNormal
mean = 1.0
sigma = 1.0
norm_p = 2
......@@ -52,21 +36,29 @@ lambda2 = 0.15
smoothing = 1.8
evtol = 1e-10
# Other
Distribution = RCB
LinearSolver = GMRES
Preconditioner = LIB_PS
# Solver
LinearReduction = 1e-12
LinearEpsilon = 1e-10
LinearSteps = 3000
NewtonSteps = 1
NewtonLineSearchSteps = 0
ConfigVerbose = 0
TimeLevel = -1
# Plotting
GeneratorPlotting = 1
MCPlotting = 1
A
# Verbose
MCVerbose = 1
PDEVerbose = 1
MainVerbose = 1
MLMCVerbose = 1
ConfigVerbose = -1
LinearVerbose = -1
MuteLevel = -1
NewtonVerbose = -1
GeneratorVerbose = -1
# Logging
TimeLevel = -1
MuteLevel = -1
DebugLevel = -1
\ No newline at end of file
......@@ -11,10 +11,7 @@ int main(int argc, char **argv) {
Config::setConfigFileName("m++.conf");
DPO dpo(&argc, argv);
std::unique_ptr<MainProgram> mainProgram =
std::make_unique<MainProgram>();
std::unique_ptr<MainProgram> mainProgram = std::make_unique<MainProgram>();
mainProgram->Initialize();
mainProgram->Run();
return 0;
return mainProgram->Run();
}
......@@ -10,13 +10,12 @@ void MainProgram::Initialize() {
problem = createStochasticProblem();
assemble = createAssemble();
// plots = make_unique<MultilevelPlotter>(*meshes);
//
// mlmc = make_unique<MultilevelMonteCarlo>(initLevels, initSampleAmount, meshes,
// stochField, assemble, plots);
mlmc = make_unique<MultilevelMonteCarlo>(initLevels, initSampleAmount,
meshes.get(), stochField.get(),
assemble.get());
}
void MainProgram::Run() {
int MainProgram::Run() {
this->PrintInfo();
meshes->PrintInfo();
// disc->PrintInfo();
......@@ -38,6 +37,7 @@ void MainProgram::Run() {
mlmc->ShowKurtosisWarning();
mlmc->ShowExponentResults();
}
return 0;
}
void MainProgram::PrintInfo() {
......@@ -80,23 +80,33 @@ void MainProgram::checkValues() {
shared_ptr<Meshes> MainProgram::createMeshes() {
if (problemName.find("1D") != string::npos)
return make_unique<Meshes>("Line", pLevel, maxLevel);
return make_shared<Meshes>("Line", pLevel, maxLevel);
if (problemName.find("2D") != string::npos)
return make_shared<Meshes>("UnitSquare", pLevel, maxLevel);
Exit("\nMesh not found in " + problemName + "\n")
}
shared_ptr<Discretization> MainProgram::createDiscretization() {
if (modelName == "LagrangeFEM" && degree == 1)
return make_shared<Discretization>("linear", 1, meshes->dim());
if (modelName == "LagrangeFEM" && degree == 2)
return make_shared<Discretization>("serendipity", 1, meshes->dim());
if (modelName == "MixedFEM" || modelName == "HybridFEM")
return make_shared<Discretization>("RT0_P0", 1, meshes->dim());
if (modelName == "DGFEM")
return make_shared<DGDiscretization>(new DGDoF(degree));
if (modelName == "DGTransport")
return make_shared<DGDiscretization>(new DGDoF(degree));
shared_ptr<IDiscretizationT<>> MainProgram::createDiscretization() {
if (modelName == "LagrangeFEM") {
if (degree == 1)
return make_shared<DiscretizationT<>>(*meshes, "P1");
if (degree == 2)
return make_shared<DiscretizationT<>>(*meshes, "P2");
}
if (modelName == "MixedFEM" || modelName == "HybridFEM") {
return make_shared<DiscretizationT<>>(*meshes, "RT0_P0");
}
if (modelName == "DGFEM" || modelName == "DGTransport") {
if (degree == 0)
return make_shared<DGDiscretization>(*meshes, "DG_P0");
if (degree == 1)
return make_shared<DGDiscretization>(*meshes, "DG_P1");
if (degree == 2)
return make_shared<DGDiscretization>(*meshes, "DG_P2");
}
Exit("\nDiscretization not found in " + modelName + "\n")
}
......@@ -170,3 +180,17 @@ shared_ptr<Assemble> MainProgram::createAssemble() {
Exit("\nAssembling for this problem with this discretization not implemented\n")
}
void MainProgram::headMLMCExperiment() {
Vout(1) << endl << "**********************************************************"
<< endl << "*** Run Multilevel Monte Carlo Method ***"
<< endl << "**********************************************************"
<< endl;
}
void MainProgram::headConvergenceTest() {
Vout(1) << endl << "**********************************************************"
<< endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
<< endl << "**********************************************************"
<< endl;
}
......@@ -4,7 +4,7 @@
#include "assemble/AssembleHeader.h"
#include "problem/StochasticEllipticProblem.hpp"
#include "MultilevelMonteCarlo.hpp"
#include "discretization/DGDiscretization.h"
#include "discretization/DGDiscretization.hpp"
#include "stochastics/StochasticField.h"
#include "utils/MultilevelPlotter.hpp"
#include <memory>
......@@ -29,14 +29,13 @@ public:
int uniformSampleAmount = -1;
shared_ptr<Meshes> meshes;
shared_ptr<Discretization> disc;
shared_ptr<IDiscretizationT<>> disc;
shared_ptr<StochasticField> stochField;
shared_ptr<StochasticProblem> problem;
shared_ptr<Assemble> assemble;
shared_ptr<MultilevelPlotter> plots;
shared_ptr<MultilevelMonteCarlo> mlmc;
unique_ptr<MultilevelMonteCarlo> mlmc;
MainProgram() {
config.get("MainVerbose", verbose);
......@@ -55,7 +54,7 @@ public:
checkValues();
}
void Run();
int Run();
void Initialize();
......@@ -66,7 +65,7 @@ private:
shared_ptr<Meshes> createMeshes();
shared_ptr<Discretization> createDiscretization();
shared_ptr<IDiscretizationT<>> createDiscretization();
shared_ptr<StochasticField> createStochasticField();
......@@ -74,19 +73,9 @@ private:
shared_ptr<Assemble> createAssemble();
static void headConvergenceTest() {
mout << endl << "**********************************************************"
<< endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
<< endl << "**********************************************************"
<< endl;
}
void headConvergenceTest();
static void headMLMCExperiment() {
mout << endl << "**********************************************************"
<< endl << "*** Run Multilevel Monte Carlo Method ***"
<< endl << "**********************************************************"
<< endl;
}
void headMLMCExperiment();
};
#endif //MLMC_MLMCMAIN_H
......@@ -8,9 +8,9 @@ MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool onlyFineLeve
auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
if (ellipticAssemble != nullptr)
return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes,
stochasticField, ellipticAssemble, plots);
stochasticField, ellipticAssemble);
if (transportAssemble != nullptr)
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes, plots,
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes,
stochasticField, transportAssemble);
Exit("Wrong Assembly Class")
}
......@@ -21,6 +21,7 @@ void MultilevelMonteCarlo::initializeMapMonteCarlo() {
bool onlyFineLevel = (i == 0) || (mcOnly);
int l = initLevels[i], M = initSampleAmount[i];
mapMonteCarlo[l] = getMonteCarlo(l, M, onlyFineLevel);
mapMonteCarlo[l]->Initialize();
}
}
......@@ -194,21 +195,21 @@ void MultilevelMonteCarlo::PrintInfo() {
}
void MultilevelMonteCarlo::ShowMCTableHead() {
mout << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc]"
" V[Qf] kurtosis cost"
<< endl << line() << endl;
Vout(1) << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc]"
" V[Qf] kurtosis cost"
<< endl << line() << endl;
}
void MultilevelMonteCarlo::ShowResultsMLMCRunHead() {
mout << endl
<< " eps value MLMC cost l"
" M"
<< endl << line() << endl;
Vout(1) << endl
<< " eps value MLMC cost l"
" M"
<< endl << line() << endl;
}
void MultilevelMonteCarlo::ShowTableBottom() {
mout << line() << endl;
Vout(1) << line() << endl;
}
string MultilevelMonteCarlo::line() {
......@@ -220,11 +221,11 @@ void MultilevelMonteCarlo::ShowMCTable() {
ShowMCTableHead();
for (auto &mc : mapMonteCarlo) {
int l = mc.first;
mout << setw(2) << l << setw(6) << mc.second->M
<< setw(12) << mc.second->avgY
<< setw(12) << mc.second->avgQ << setw(12) << mc.second->varY
<< setw(12) << mc.second->varQ << setw(12) << mc.second->kurtosis
<< setw(12) << mc.second->avgCost << endl;
Vout(1) << setw(2) << l << setw(6) << mc.second->M
<< setw(12) << mc.second->avgY
<< setw(12) << mc.second->avgQ << setw(12) << mc.second->varY
<< setw(12) << mc.second->varQ << setw(12) << mc.second->kurtosis
<< setw(12) << mc.second->avgCost << endl;
}
ShowTableBottom();
}
......@@ -232,36 +233,36 @@ void MultilevelMonteCarlo::ShowMCTable() {
void MultilevelMonteCarlo::ShowKurtosisWarning() {
for (auto &mc : mapMonteCarlo) {
if (mc.second->kurtosis > 100.0) {
mout << endl << "WARNING: kurtosis on level " << setw(2) << mc.first
<< " is " << setw(6) << mc.second->kurtosis << endl
<< "indicates MLMC correction is dominated by a few rare paths!"
<< endl;
Vout(1) << endl << "WARNING: kurtosis on level " << setw(2) << mc.first
<< " is " << setw(6) << mc.second->kurtosis << endl
<< "indicates MLMC correction is dominated by a few rare paths!"
<< endl;
}
}
}
void MultilevelMonteCarlo::ShowExponentResults() {
estimateExponents();
mout << endl << "alpha = " << setw(10) << alpha
<< " (exponent for the weak convergence of the FEM)"
<< endl << " beta = " << setw(10) << beta
<< " (exponent for the decay of the variance)"
<< endl << "gamma = " << setw(10) << gamma
<< " (exponent for the growth of the cost)"
<< endl << "numErr = " << setw(10) << numErr
<< " (estimated numerical error)"
<< endl << "statErr = " << setw(10) << statErr
<< " (estimated statistical error)"
<< endl << "error = " << setw(10) << totalErr
<< " (estimated total error)"
<< endl;
Vout(1) << endl << "alpha = " << setw(10) << alpha
<< " (exponent for the weak convergence of the FEM)"
<< endl << " beta = " << setw(10) << beta
<< " (exponent for the decay of the variance)"
<< endl << "gamma = " << setw(10) << gamma
<< " (exponent for the growth of the cost)"
<< endl << "numErr = " << setw(10) << numErr
<< " (estimated numerical error)"
<< endl << "statErr = " << setw(10) << statErr
<< " (estimated statistical error)"
<< endl << "error = " << setw(10) << totalErr
<< " (estimated total error)"
<< endl;
}
void MultilevelMonteCarlo::ShowResultsMLMCRun(double eps) {
ShowResultsMLMCRunHead();
string level_str = vec2str(levels);
string sample_str = vec2str(numsamples);
mout << setw(6) << eps << setw(12) << value << setw(12) << cost
Vout(1) << setw(6) << eps << setw(12) << value << setw(12) << cost
<< setw(14) << level_str << setw(36) << sample_str << endl;
ShowTableBottom();
}
......@@ -271,8 +272,8 @@ void MultilevelMonteCarlo::ShowResultsOverEpsilon(const vector<double> &epsLst)
int ctr = 0;
for (auto &eps: epsLst) {
string level_str, sample_str;
mout << setw(6) << eps << setw(12) << valueOverEpsilon[ctr]
<< setw(12) << costOverEpsilon[ctr];
Vout(1) << setw(6) << eps << setw(12) << valueOverEpsilon[ctr]
<< setw(12) << costOverEpsilon[ctr];
for (unsigned long i = 0; i < levelsOverEpsilon[ctr].size(); i++) {
if (i != levelsOverEpsilon[ctr].size() - 1)
level_str += to_string(levelsOverEpsilon[ctr][i]) + " ";
......@@ -280,14 +281,14 @@ void MultilevelMonteCarlo::ShowResultsOverEpsilon(const vector<double> &epsLst)
level_str += to_string(levelsOverEpsilon[ctr][i]);
}
mout << setw(14) << level_str;
Vout(1) << setw(14) << level_str;
for (unsigned long i = 0; i < sampleAmountOverEpsilon[ctr].size(); i++) {
if (i != sampleAmountOverEpsilon[ctr].size() - 1)
sample_str += to_string(sampleAmountOverEpsilon[ctr][i]) + " ";
else
sample_str += to_string(sampleAmountOverEpsilon[ctr][i]);
}
mout << setw(36) << sample_str << endl;
Vout(1) << setw(36) << sample_str << endl;
ctr++;
}
ShowTableBottom();
......
......@@ -7,37 +7,10 @@
class MultilevelMonteCarlo {
private:
void initializeMapMonteCarlo();
void initializeValues();
MonteCarlo *getMonteCarlo(int l, int dM, bool onlyFineLevel);
void clearMapMonteCarlo();
void wrapUpResults(double &val,
double &mlmcCost,
vector<int> &listLevels,
vector<int> &samples);
void estimateExponents(bool excludeBaseLevel = true);
bool checkForNewLevel();
void updateSampleAmount(const double &eps);
void estimateNumericalError();
void estimateStatisticalError();
void appendLevel(double eps);
static string line();
public:
int plotting = 0;
int verbose = 0;
int verbose = -1;
int plotting = -1;
IndentedLogger *logger;
int maxLevel, pLevel;
......@@ -45,7 +18,6 @@ public:
map<int, MonteCarlo *> mapMonteCarlo;
Meshes *meshes;
MultilevelPlotter &plots;
StochasticField *stochasticField;
Assemble *assemble;
......@@ -62,13 +34,11 @@ public:
MultilevelMonteCarlo(vector<int> &initLevels,
vector<int> &initSampleAmount,
Meshes *meshes,
MultilevelPlotter &plots,
StochasticField *stochasticField,
Assemble *assemble) :
initLevels(initLevels),
initSampleAmount(initSampleAmount),
meshes(meshes),
plots(plots),
stochasticField(stochasticField),
assemble(assemble) {
......@@ -85,7 +55,6 @@ public:
}
~MultilevelMonteCarlo() {
delete logger;
clearMapMonteCarlo();
}
......@@ -97,11 +66,11 @@ public:
void PrintInfo();
static void ShowMCTableHead();
void ShowMCTableHead();
static void ShowResultsMLMCRunHead();
void ShowResultsMLMCRunHead();
static void ShowTableBottom();
void ShowTableBottom();
void ShowMCTable();
......@@ -112,6 +81,34 @@ public:
void ShowResultsMLMCRun(double eps);
void ShowResultsOverEpsilon(const vector<double> &epsLst);
private:
void initializeMapMonteCarlo();
void initializeValues();
MonteCarlo *getMonteCarlo(int l, int dM, bool onlyFineLevel);
void clearMapMonteCarlo();
void wrapUpResults(double &val,
double &mlmcCost,
vector<int> &listLevels,
vector<int> &samples);
void estimateExponents(bool excludeBaseLevel = true);
bool checkForNewLevel();
void updateSampleAmount(const double &eps);
void estimateNumericalError();
void estimateStatisticalError();
void appendLevel(double eps);
static string line();
};
#endif //MULTILEVELMONTECARLO_HPP
......@@ -11,13 +11,8 @@
class MonteCarlo {
private:
public:
int plotting = 0;
Plot *coarsePlot = nullptr;
Plot *finePlot = nullptr;
int verbose = 0;
IndentedLogger *logger = nullptr;
......@@ -43,38 +38,28 @@ public:
double fineQ = 0.0, coarseQ = 0.0;
double fineCost = 0.0, coarseCost = 0.0;
MonteCarlo(int l,
int dM,
bool onlyFineLevel,
Meshes *meshes,
MultilevelPlotter &plots,
StochasticField *stochFields) :
MonteCarlo(int l, int dM, bool onlyFineLevel,
Meshes *meshes, StochasticField *stochFields) :
l(l), dM(dM), onlyFineLevel(onlyFineLevel),
meshes(meshes), pLevel(meshes->pLevel()),
stochasticField(stochFields) {
config.get("MCPlotting", plotting);
config.get("functional", functional);
// coarsePlot = plots.GetPlot(l - 1);
// finePlot = plots.GetPlot(l);
config.get("MCVerbose", verbose);
config.get("Functional", functional);
logger = IndentedLogger::GetInstance();
}
void Method();
protected:
virtual void initialize() = 0;
virtual void Initialize() = 0;
protected:
virtual void method() = 0;
virtual void solvePDE(int l,
int m,
bool fineLevel,
double &valueQ,
double &cost,
virtual void solvePDE(int l, int m, bool fineLevel,
double &valueQ, double &cost,
Vector &solution) = 0;
void updateSums(double fineCost, double diffQ, double fineQ);
......
......@@ -3,7 +3,7 @@
using namespace std;
void MonteCarloElliptic::initialize() {
void MonteCarloElliptic::Initialize() {
finePermeability = 0.0, coarsePermeability = 0.0;
fineSolution = 0.0, coarseSolution = 0.0;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
......@@ -11,36 +11,37 @@ void MonteCarloElliptic::initialize() {
void MonteCarloElliptic::method() {
for (int m = M; m < M + dM; m++) {
// Refactor dir stuff
string sampleDirName = mkSampleDirName(m, true);
mkSampleDir(m, sampleDirName);
stochasticField->SetSampleDir(l, sampleDirName);
stochasticField->SetPlot(l, finePlot);