Commit 4cf1e03a authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

new feature fit in

parent d108ed62
......@@ -20,7 +20,7 @@ Functional = L2
#Functional = Goal
# ----- Multilevel Monte Carlo -----
maxLevel = 9
maxLevel = 7
epsilon = 0.01
mcOnly = false
initLevels = [3, 4, 5]
......@@ -49,14 +49,15 @@ GeneratorPlotting = 1
MCPlotting = 1
# ----- Verbose -----
MCVerbose = 1
PDEVerbose = 1
MainVerbose = 1
MLMCVerbose = 1
ConfigVerbose = -1
LinearVerbose = -1
NewtonVerbose = -1
GeneratorVerbose = -1
MCVerbose =1
PDEVerbose =1
MainVerbose =1
MeshVerbose =1
MLMCVerbose =1
ConfigVerbose =1
LinearVerbose = 1
NewtonVerbose = 1
GeneratorVerbose =1
# ----- Logging -----
TimeLevel = -1
......
......@@ -2,11 +2,11 @@
Experiment = MLMCExperiment
#Experiment = ConvergenceTest
Problem = StochasticPollution1D
#Problem = StochasticPollution1D
#Problem = StochasticPollutionCosHat1D
#Problem = StochasticPollution2D
#Problem = StochasticPollutionCosHat2D
#Problem = StochasticPollutionMollifiedBar2D
Problem = StochasticPollutionMollifiedBar2D
Model = DGTransport
#Model = PGTransport
......@@ -33,7 +33,7 @@ gamma = 0.01
Keps = 1e-5
# ----- Multilevel Monte Carlo -----
maxLevel = 9
maxLevel = 7
epsilon = 0.01
mcOnly = false
initLevels = [4, 5, 6]
......@@ -66,10 +66,11 @@ MCVerbose = 1
PDEVerbose = 1
MainVerbose = 1
MLMCVerbose = 1
ConfigVerbose = -1
ConfigVerbose = 0
LinearVerbose = -1
NewtonVerbose = -1
GeneratorVerbose = -1
NewtonVerbose = 1
GeneratorVerbose = 1
TimeIntegratorVerbose = 1
# ----- Logging -----
TimeLevel = -1
......
#include "m++.hpp"
#include "MainProgram.hpp"
#include "Parallel.h"
#ifdef COMPLEX
#error undef COMPLEX in src/Compiler.h
#endif
int main(int argc, char **argv) {
PPM = new ParallelProgrammingModel(&argc, argv);
logging = new Logging;
Config::setSearchPath("../mlmc/");
Config::setConfigFileName("m++.conf");
Config::setLogging(logging);
config = Config::parseConfig(&argc, argv);
logging->initialize();
config.get("DebugLevel", DebugLevel);
config.get("TimeLevel", TimeLevel);
Mpp::initialize(&argc, argv);
std::unique_ptr<MainProgram> mainProgram = std::make_unique<MainProgram>();
mainProgram->Initialize();
int rc = mainProgram->Run();
delete logging;
delete PPM;
clearQuadrature();
return rc;
}
......@@ -42,17 +42,30 @@ int MainProgram::Run() {
}
void MainProgram::PrintInfo() {
Vout(1) << endl
vout(1) << endl
<< "Main Info:" << endl
<< " Problem: " << problemName << endl
<< " Model: " << modelName << endl
<< " StochField: " << fieldName << endl
<< " Experiment: " << experimentName << endl;
if (experimentName == "MLMCExperiment")
Vout(1) << " epsilon: " << epsilon << endl;
Vout(1) << endl;
vout(1) << " epsilon: " << epsilon << endl;
vout(1) << endl;
}
//void NavierStokesData::printInfo() {
// mout << "Navier-Stokes Equations on an Unbounded Strip with Obstacle" << endl;
// problem->printInfo();
//
// mout.PrintInfo("Navier-Stokes", verbose,
// PrintInfoEntry("Approximation model", solutionModel),
// PrintInfoEntry("Defect bound model", defectModel),
// PrintInfoEntry("Norm bound model", normModel),
// PrintInfoEntry("Plot:", plotSolution),
// PrintInfoEntry("Verified computation", verified));
//}
void MainProgram::checkValues() {
if (verbose == -1) Exit("\nPlease set MainVerbose\n")
if (problemName.empty()) Exit("\nPlease set Problem\n")
......@@ -184,14 +197,14 @@ shared_ptr<Assemble> MainProgram::createAssemble() {
}
void MainProgram::headMLMCExperiment() {
Vout(1) << endl << "**********************************************************"
vout(1) << endl << "**********************************************************"
<< endl << "*** Run Multilevel Monte Carlo Method ***"
<< endl << "**********************************************************"
<< endl;
}
void MainProgram::headConvergenceTest() {
Vout(1) << endl << "**********************************************************"
vout(1) << endl << "**********************************************************"
<< endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
<< endl << "**********************************************************"
<< endl;
......
......@@ -40,15 +40,15 @@ void MultilevelMonteCarlo::clearMapMonteCarlo() {
}
void MultilevelMonteCarlo::Method() {
logger->StartMethod("Start MLMC Method", verbose);
mout.StartBlock("MLMC Method");
for (auto &mc : mapMonteCarlo)
mc.second->Method();
logger->EndMethod(verbose);
mout.EndBlock(verbose == 0);
}
void MultilevelMonteCarlo::Method(const double eps) {
logger->StartMethod("Start MLMC Method", verbose);
logger->InnerMsg("eps: " + to_string(eps), verbose);
mout.StartBlock("MLMC Method");
vout(1) << "eps=" << to_string(eps) << endl;
bool converged = false;
while (!converged) {
for (auto &mc : mapMonteCarlo)
......@@ -68,11 +68,12 @@ void MultilevelMonteCarlo::Method(const double eps) {
}
}
wrapUpResults(value, cost, levels, numsamples);
logger->EndMethod(verbose);
mout.EndBlock(verbose == 0);
}
void MultilevelMonteCarlo::Method(const vector<double> &epsLst) {
logger->StartMethod("Start MLMC Method over epsilon", verbose);
mout.StartBlock("Multiple MLMC Methods");
mout << "eps=[" << vec2str(epsLst) << "]" << endl;
for (auto &eps : epsLst) {
initializeMapMonteCarlo();
initializeValues();
......@@ -82,7 +83,7 @@ void MultilevelMonteCarlo::Method(const vector<double> &epsLst) {
valueOverEpsilon.push_back(value);
costOverEpsilon.push_back(cost);
}
logger->EndMethod(verbose);
mout.EndBlock(verbose == 0);
}
void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
......@@ -105,9 +106,7 @@ void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
res = linear_fit(levelList, costs);
gamma = res.first;
logger->InnerMsg("alpha: " + to_string(alpha) +
" beta: " + to_string(beta) +
" gamma: " + to_string(gamma), verbose);
vout(1) << "alpha=" << alpha << " beta=" << beta << " gamma=" << gamma << endl;
}
void MultilevelMonteCarlo::updateSampleAmount(const double &eps) {
......@@ -124,7 +123,7 @@ void MultilevelMonteCarlo::updateSampleAmount(const double &eps) {
string msg = "dM:";
for (auto &mc:mapMonteCarlo) msg += " " + to_string(mc.second->dM);
logger->InnerMsg(msg, verbose);
vout(1) << msg << endl;
}
bool MultilevelMonteCarlo::checkForNewLevel() {
......@@ -147,8 +146,7 @@ void MultilevelMonteCarlo::estimateNumericalError() {
numErr);
exponent--;
}
logger->InnerMsg("estimated numerical error: " + to_string(numErr), verbose);
vout(1) << "numErr=" << to_string(numErr) << endl;
}
void MultilevelMonteCarlo::estimateStatisticalError() {
......@@ -164,7 +162,7 @@ void MultilevelMonteCarlo::appendLevel(double eps) {
if (L > maxLevel) {
Exit ("Maximum level has been reached.")
} else {
logger->InnerMsg("new level: " + to_string(L), verbose);
vout(1) << "L_new=" << to_string(L) << endl;
mapMonteCarlo[L] = getMonteCarlo(L, 0, false);
mapMonteCarlo[L]->varY = mapMonteCarlo[L - 1]->varY / pow(2.0, beta);
mapMonteCarlo[L]->avgCost = mapMonteCarlo[L - 1]->avgCost * pow(2.0, gamma);
......@@ -186,7 +184,7 @@ void MultilevelMonteCarlo::wrapUpResults(double &val,
}
void MultilevelMonteCarlo::PrintInfo() {
Vout(1) << endl
vout(1) << endl
<< "MLMC Info:" << endl
<< " initLevels: " << vec2str(initLevels) << endl
<< " initSample: " << vec2str(initSampleAmount) << endl
......@@ -195,21 +193,21 @@ void MultilevelMonteCarlo::PrintInfo() {
}
void MultilevelMonteCarlo::ShowMCTableHead() {
Vout(1) << endl
vout(1) << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc]"
" V[Qf] kurtosis cost"
<< endl << line() << endl;
}
void MultilevelMonteCarlo::ShowResultsMLMCRunHead() {
Vout(1) << endl
vout(1) << endl
<< " eps value MLMC cost l"
" M"
<< endl << line() << endl;
}
void MultilevelMonteCarlo::ShowTableBottom() {
Vout(1) << line() << endl;
vout(1) << line() << endl;
}
string MultilevelMonteCarlo::line() {
......@@ -221,7 +219,7 @@ void MultilevelMonteCarlo::ShowMCTable() {
ShowMCTableHead();
for (auto &mc : mapMonteCarlo) {
int l = mc.first;
Vout(1) << setw(2) << l << setw(6) << mc.second->M
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
......@@ -233,7 +231,7 @@ void MultilevelMonteCarlo::ShowMCTable() {
void MultilevelMonteCarlo::ShowKurtosisWarning() {
for (auto &mc : mapMonteCarlo) {
if (mc.second->kurtosis > 100.0) {
Vout(1) << endl << "WARNING: kurtosis on level " << setw(2) << mc.first
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;
......@@ -243,7 +241,7 @@ void MultilevelMonteCarlo::ShowKurtosisWarning() {
void MultilevelMonteCarlo::ShowExponentResults() {
estimateExponents();
Vout(1) << endl << "alpha = " << setw(10) << alpha
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)"
......@@ -262,7 +260,7 @@ void MultilevelMonteCarlo::ShowResultsMLMCRun(double eps) {
ShowResultsMLMCRunHead();
string level_str = vec2str(levels);
string sample_str = vec2str(numsamples);
Vout(1) << 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();
}
......@@ -272,7 +270,7 @@ void MultilevelMonteCarlo::ShowResultsOverEpsilon(const vector<double> &epsLst)
int ctr = 0;
for (auto &eps: epsLst) {
string level_str, sample_str;
Vout(1) << setw(6) << eps << setw(12) << valueOverEpsilon[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)
......@@ -281,14 +279,14 @@ void MultilevelMonteCarlo::ShowResultsOverEpsilon(const vector<double> &epsLst)
level_str += to_string(levelsOverEpsilon[ctr][i]);
}
Vout(1) << 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]);
}
Vout(1) << setw(36) << sample_str << endl;
vout(1) << setw(36) << sample_str << endl;
ctr++;
}
ShowTableBottom();
......
#ifndef MULTILEVELMONTECARLO_HPP
#define MULTILEVELMONTECARLO_HPP
#include "utils/IndentedLogger.hpp"
#include "mc/MonteCarloElliptic.hpp"
#include "mc/MonteCarloTransport.hpp"
......@@ -11,8 +10,6 @@ public:
int verbose = -1;
int plotting = -1;
IndentedLogger *logger;
int maxLevel, pLevel;
vector<int> &initLevels, &initSampleAmount;
map<int, MonteCarlo *> mapMonteCarlo;
......@@ -44,9 +41,7 @@ public:
config.get("MLMCPlotting", plotting);
config.get("mcOnly", mcOnly);
config.get("MLMCVerbose", verbose);
logger = IndentedLogger::GetInstance();
maxLevel = meshes->Level();
pLevel = meshes->pLevel();
......
......@@ -2,7 +2,6 @@
#define TUTORIAL_DGTRANSPORT_HPP
#include "DGTAssemble.h"
#include "utils/IndentedLogger.hpp"
#include <memory>
......@@ -11,8 +10,6 @@ public:
// std::shared_ptr<IDiscretizationT<>> disc;
std::shared_ptr<StochasticTransportProblem> problem;
IndentedLogger *logger;
double flux_alpha = 1.0;
double diffusion = 0.0;
......@@ -22,7 +19,6 @@ public:
: DGTAssemble(dynamic_cast<DGDiscretization *>(disc.get())),
problem(move(problem)) {
config.get("flux_alpha", flux_alpha);
logger = IndentedLogger::GetInstance();
}
const char *Name() const override { return "DGTransportAssemble"; }
......@@ -234,17 +230,35 @@ public:
void PrintStep(double t, const Vector &u) const override {
std::pair<double, double> rate = InFlowOutFlowRate(u);
string msg = "Step(";
msg = msg.append(to_string(step)).append("), t(");
msg = msg.append(to_string(t)).append("): Energy(u) = ");
msg = msg.append(to_string(Energy(u))).append(" Mass(u) = ");
msg = msg.append(to_string(Mass(u)));
// config.get("precision")
// std::setw(precision) <<
vout(1) << "t(" << step << ")=" << std::setw(6) << t
<< " Energy(u)=" << std::setw(6) << Energy(u)
<< " Mass(u)=" << std::setw(6) << Mass(u);
if (abs(rate.first) > 1e-7)
msg = msg.append(" InFlowRate(u) = ").append(to_string(rate.first));
vout(1) << " IFR(u)=" << rate.first;
if (abs(rate.second) > 1e-7)
msg = msg.append(" OutFlowRate(u) = ").append(to_string(rate.second));
vout(1) << " OFR(u)=" << std::setw(6) << rate.second;
vout(4) << std::flush;
vout(3) << std::flush << std::beginl;
// if (verbose == 0) return;
// if (verbose == 1) {
// mout << "Step " << step << std::flush << std::beginl;
// } else if (verbose == 2) {
// mout << "lambda(" << step << ")= ";
// print(lambda);
// mout << std::flush << std::beginl;
// } else {
// mout << "lambda(" << step << ")= ";
// print(lambda);
// mout << endl;
// }
logger->LogMsgFlush(msg, verbose);
// vout(1) << endl;
}
void VtkPlotting(double t, const Vector &u) const override {
......@@ -320,7 +334,7 @@ public:
}
void PrintInfo(const Vector &u) const override {
Vout (1) << endl
vout (1) << endl
<< "Assemble Info:" << endl
<< " Name: " << Name() << endl
<< " Problem: " << problem->Name() << endl
......
......@@ -40,7 +40,7 @@ public:
const char *Name() const override { return "Lagrange"; }
void PrintInfo(const Vector &u) const override {
Vout (1) << endl
vout (1) << endl
<< "Assemble Info:" << endl
<< " Name: " << Name() << endl
<< " Problem: " << problem->Name() << endl
......
......@@ -2,20 +2,23 @@
void MonteCarlo::Method() {
logger->StartMethod("Start MC Method", verbose);
if (!onlyFineLevel)
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM), verbose);
else
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM)
+ " onlyFineLevel: True", verbose);
mout.StartBlock("MC Method");
vout(1) << "l=" << l << " dM=" << dM << " fine=" << onlyFineLevel << endl;
method();
updateAvg();
updateStatistic();
vout(1) << "|E[Y]|=" << avgY << " V[Y]=" << varY << endl;
mout.EndBlock(verbose == 0);
}
logger->InnerMsg("|E[Y]|: " + to_string(avgY) + " V[Y]: " + to_string(varY),
verbose);
logger->EndMethod(verbose);
void MonteCarlo::SolvePDE(int l, int m, bool fine,
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) << "Q=" << valueQ << " cost=" << cost << endl;
mout.EndBlock(verbose == 0);
}
void MonteCarlo::updateSums(double fineCost, double diffQ, double fineQ) {
......
......@@ -2,7 +2,6 @@
#define MLMC_MC_HPP
#include "stochastics/StochasticField.h"
#include "utils/IndentedLogger.hpp"
#include <utility>
#include <utils/MultilevelPlotter.hpp>
......@@ -12,7 +11,6 @@ class MonteCarlo {
public:
int plotting = 0;
int verbose = 0;
IndentedLogger *logger = nullptr;
int l = 0;
int pLevel = 0;
......@@ -45,18 +43,20 @@ public:
config.get("MCPlotting", plotting);
config.get("MCVerbose", verbose);
config.get("Functional", functional);
logger = IndentedLogger::GetInstance();
}
void Method();
void SolvePDE(int l, int m, bool fine,
double &valueQ, double &cost,
Vector &solution);
virtual void Initialize() = 0;
protected:
virtual void method() = 0;
virtual void solvePDE(int l, int m, bool fineLevel,
virtual void solvePDE(int l, int m, bool fine,
double &valueQ, double &cost,
Vector &solution) = 0;
......
......@@ -17,7 +17,7 @@ void MonteCarloElliptic::method() {
stochasticField->GetFineSample(l, fineInput);
assemble->problem->SetSample(&fineInput);
solvePDE(l, m, true, fineQ, fineCost, fineSolution);
SolvePDE(l, m, true, fineQ, fineCost, fineSolution);
if (plotting)
plotter->PlotVector("u", fineSolution,
......@@ -32,7 +32,7 @@ void MonteCarloElliptic::method() {
stochasticField->GetCoarseSample(l, fineInput, coarseInput);
assemble->problem->SetSample(&coarseInput);
solvePDE(l, 0, false,
SolvePDE(l, 0, false,
coarseQ, coarseCost, coarseSolution);
if (plotting)
......
......@@ -17,7 +17,7 @@ void MonteCarloTransport::method() {
stochasticField->GetFineSample(l, fineInput);
assemble->problem->SetSample(&fineInput);
solvePDE(l, m, true, fineQ, fineCost, fineSolution);
SolvePDE(l, m, true, fineQ, fineCost, fineSolution);
// if (plotting)
// plotter->PlotVector("u", fineSolution,
......@@ -32,7 +32,7 @@ void MonteCarloTransport::method() {
stochasticField->GetCoarseSample(l, fineInput, coarseInput);
assemble->problem->SetSample(&coarseInput);
solvePDE(l - 1, m, false,
SolvePDE(l - 1, m, false,
coarseQ, coarseCost, coarseSolution);
// if (plotting)
......@@ -45,9 +45,6 @@ void MonteCarloTransport::method() {
void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel,
double &valueQ, double &cost, Vector &solution) {
logger->StartMethod("Start Solving PDE", verbose);
logger->InnerMsg("l: " + to_string(l) + " m: " + to_string(m), verbose);
Preconditioner *pc = GetPC("PointBlockGaussSeidel");
Solver solver(pc, "GMRES");
......@@ -63,8 +60,6 @@ void MonteCarloTransport::solvePDE(int l, int m, bool fineLevel,
if (functional == "Energy") valueQ = assemble->Energy(solution);
if (functional == "Outflow")
valueQ = assemble->InFlowOutFlowRate(solution).second;
logger->EndMethod(verbose);
}
TimeSeries MonteCarloTransport::getTimeSeries(bool fineLevel) {
......
#ifndef M_COVARIANCEFUNCTION_H
#define M_COVARIANCEFUNCTION_H
#include "Config.hpp"
#include "Debug.h"
#include "utility/Config.hpp"
class CovarianceFunction {
......
......@@ -11,9 +11,9 @@ void HybridFluxGenerator::generateFineSample(Vector &fineField) {
assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, fineField);
if (plotting)
plotter->PlotVector("flux", *flux,
3, l, "CellData");
// if (plotting)
// plotter->PlotVector("flux", *flux,
// 3, l, "CellData");
}
void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
......@@ -26,9 +26,9 @@ void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
assemble->setFlux(*u, *flux);
assemble->SetNormalFlux(*u, coarseField);
if (plotting)
plotter->PlotVector("flux", *flux,
3, l - 1, "CellData");
// if (plotting)
// plotter->PlotVector("flux", *flux,