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

major refactoring and implementations of dg for elliptic and transport

parent 347756ff
......@@ -9,7 +9,7 @@ Preconditioner = LIB_PS
LinearReduction = 1e-12
LinearEpsilon = 1e-10
LinearSteps = 2000
LinearSteps = 3000
NewtonSteps = 1
NewtonLineSearchSteps = 0
......
Model = LinearFEM
#Model = LinearFEM
#Model = QuadraticFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = DGFEM
Model = FVFEM
#Model = LinearDGFEM
#Model = QuadraticDGFEM
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
......@@ -10,17 +12,17 @@ Problem = StochasticLaplace2D
StochasticField = LogNormal
#StochasticField = Gaussian
Experiment = ConvergenceTest
#Experiment = ConvergenceTest
#Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
Experiment = MLMCOverEpsilon
init = 2
l_init = 4 5
M_init = 1 1
init = 3
l_init = 3 4 5
M_init = 8 4 2
Lmax = 8
epsilon = 0.001
epsilon_lst = 0.01 0.005 0.001
epsilon = 0.01
epsilon_lst = 0.01 0.005 0.003
plevel = 2
functional = L2
......@@ -30,8 +32,5 @@ functional = L2
MCVerbose = 1
MLMCVerbose = 2
MCPlotting = 1
MCPlotting = 0
MLMCPlotting = 0
vtkplot = 1
gnuplot = 1
......@@ -2,12 +2,20 @@ set(MLMC_SRC
Main.C
Utils.C
MLMCMain.C
RNManager.C
MonteCarlo.C
StochasticFields.C
StochasticProblem.C
CirculantEmbedding.C
MultilevelMonteCarlo.C
mc/MonteCarloElliptic.C
mc/MonteCarloTransport.C
problem/StochasticEllipticProblem.h
problem/StochasticTransportProblem.C
assemble/EllipticAssemble.C
assemble/MixedEllipticAssemble.C
assemble/HybridEllipticAssemble.C
assemble/DGEllipticAssemble.C
stochfields/RNManager.C
stochfields/StochasticFields.C
stochfields/CirculantEmbedding.C
discgalerkin/DGDiscretization.C
discgalerkin/DGMatrixGraph.C
)
add_library(MLMC STATIC ${MLMC_SRC})
#include "MLMCMain.h"
void MLMCMain() {
double epsilon;
int init, plevel, Lmax;
vector<int> l_init, M_init;
vector<double> epsilon_lst;
string problemName, modelName, fieldName, experimentName;
using namespace std;
ReadConfig(Settings, "epsilon", epsilon);
ReadConfig(Settings, "plevel", plevel);
ReadConfig(Settings, "Lmax", Lmax);
ReadConfig(Settings, "init", init);
ReadConfig(Settings, "l_init", l_init, init);
ReadConfig(Settings, "M_init", M_init, init);
ReadConfig(Settings, "epsilon_lst", epsilon_lst, init);
ReadConfig(Settings, "Problem", problemName);
ReadConfig(Settings, "Model", modelName);
ReadConfig(Settings, "StochasticField", fieldName);
ReadConfig(Settings, "Experiment", experimentName);
void MLMCMain::initialize() {
meshes = getMeshes();
disc = getDiscretization();
stochFields = getStochasticFields();
problem = getStochasticProblem();
assemble = getAssemble();
graphs = getMatrixGraphs();
Meshes meshes = getMeshes(problemName, plevel, Lmax);
Discretization disc = getDiscretization(modelName, meshes);
StochasticFields stochasticFields = StochasticFields(fieldName, meshes);
StochasticProblem *problem = getStochasticProblem(stochasticFields, problemName);
Assemble *assemble = getAssemble(problemName, modelName, disc, problem);
MatrixGraphs graphs = getMatrixGraphs(modelName, meshes, disc);
mlmc = new MultilevelMonteCarlo(l_init, M_init, meshes,
stochFields, assemble, graphs);
}
MultilevelMonteCarlo mlmc(l_init, M_init, meshes,
stochasticFields, assemble, graphs);
Meshes *MLMCMain::getMeshes() {
if (problemName.find("1D") != string::npos)
return new Meshes("Line", plevel, Lmax);
if (problemName.find("2D") != string::npos)
return new Meshes("UnitSquare", plevel, Lmax);
Exit("Not implemented yet")
return nullptr;
}
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
mlmc.method();
mlmc.showMCTable();
mlmc.showKurtosisWarning();
mlmc.showExponentResults();
} else if (experimentName == "MLMCExperiment") {
headMLMCExperiment();
mlmc.method(epsilon);
mlmc.showMCTable();
mlmc.showResultsMLMCRun(epsilon);
mlmc.showKurtosisWarning();
mlmc.showExponentResults();
} else if (experimentName == "MLMCOverEpsilon") {
headMLMCOverEpsilon();
mlmc.method(epsilon_lst);
mlmc.showResultsOverEpsilon(epsilon_lst);
}
Discretization *MLMCMain::getDiscretization() {
if (modelName == "LinearFEM")
return new Discretization("linear", 1, meshes->dim());
if (modelName == "QuadraticFEM")
return new Discretization("quadratic", 1, meshes->dim());
if (modelName == "MixedFEM" || modelName == "HybridFEM")
return new Discretization("RT0_P0", 1, meshes->dim());
if (modelName == "FVFEM")
return new DGDiscretization(new DGDoF(0));
if (modelName == "LinearDGFEM")
return new DGDiscretization(new DGDoF(1));
if (modelName == "QuadraticDGFEM")
return new DGDiscretization(new DGDoF(2));
Exit("Not implemented yet")
return nullptr;
}
MatrixGraphs
getMatrixGraphs(const string &modelName, Meshes &meshes, Discretization &disc) {
if (modelName == "LinearFEM" || modelName == "QuadraticFEM")
return MatrixGraphs(meshes, disc);
if (modelName == "HybridFEM")
return MatrixGraphs(meshes, dof("face", 1));
if (modelName == "DGFEM") Exit("Not implemented yet")
StochasticFields *MLMCMain::getStochasticFields() {
return new StochasticFields(fieldName, *meshes);
}
StochasticProblem *MLMCMain::getStochasticProblem() {
if (problemName == "StochasticLaplace1D")
return new StochasticLaplace1D(*stochFields);
if (problemName == "StochasticLaplace2D")
return new StochasticLaplace2D(*stochFields);
Exit("Not implemented yet")
return nullptr;
}
Assemble *getAssemble(const string &problemName, const string &modelName,
Discretization &disc, StochasticProblem *problem) {
Assemble *MLMCMain::getAssemble() {
if (problemName.find("Laplace") != string::npos) {
if (modelName == "LinearFEM" || modelName == "QuadraticFEM")
return new EllipticAssemble(disc, *problem);
return new EllipticAssemble(disc, problem);
if (modelName == "MixedFEM")
return new MixedEllipticAssemble(disc, *problem);
return new MixedEllipticAssemble(disc, problem);
if (modelName == "HybridFEM")
return new HybridEllipticAssemble(disc, *problem);
if (modelName == "DGFEM") Exit("Not implemented yet")
return nullptr;
return new HybridEllipticAssemble(disc, problem);
if (modelName == "FVFEM" || modelName == "LinearDGFEM" ||
modelName == "QuadraticDGFEM")
return new DGEllipticAssemble(disc, problem);
} else if (problemName.find("Transport") != string::npos) {
if (modelName == "DGFEM") Exit("Not implemented yet")
return nullptr;
if (modelName == "FVFEM")
return nullptr;
if (modelName == "LinearDGFEM")
return nullptr;
if (modelName == "QuadraticDGFEM")
return nullptr;
}
Exit("Not implemented yet")
return nullptr;
}
Meshes getMeshes(const string &problemName, int plevel, int Lmax) {
if (problemName.find("1D") != string::npos)
return Meshes("Line", plevel, Lmax);
if (problemName.find("2D") != string::npos)
return Meshes("UnitSquare", plevel, Lmax);
}
Discretization getDiscretization(const string &modelName, Meshes &meshes) {
if (modelName == "LinearFEM")
return Discretization("linear", 1, meshes.dim());
if (modelName == "QuadraticFEM")
return Discretization("quadratic", 1, meshes.dim());
if (modelName == "MixedFEM")
return Discretization("RT0_P0", 1, meshes.dim());
MatrixGraphs *MLMCMain::getMatrixGraphs() {
if (modelName == "LinearFEM" || modelName == "QuadraticFEM" ||
modelName == "MixedFEM")
return new MatrixGraphs(*meshes, *disc);
if (modelName == "HybridFEM")
return Discretization("RT0_P0", 1, meshes.dim());
}
void headConvergenceTest() {
mout << endl << "**********************************************************"
<< endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
<< endl << "**********************************************************" << endl;
}
void headMLMCExperiment() {
mout << endl << "**********************************************************"
<< endl << "*** MultilevelMonteCarlo test run ***"
<< endl << "**********************************************************" << endl;
return new MatrixGraphs(*meshes, dof("face", 1));
if (modelName == "FVFEM" || modelName == "LinearDGFEM" ||
modelName == "QuadraticDGFEM")
return new CellMatrixGraphs(*meshes, *disc);
Exit("Not implemented yet")
return nullptr;
}
void headMLMCOverEpsilon() {
mout << endl << "**********************************************************"
<< endl << "*** MultilevelMonteCarlo over epsilon ***"
<< endl << "**********************************************************" << endl;
}
void MLMCMain::run() {
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
mlmc->method();
mlmc->showMCTable();
mlmc->showKurtosisWarning();
mlmc->showExponentResults();
} else if (experimentName == "MLMCExperiment") {
headMLMCExperiment();
mlmc->method(epsilon);
mlmc->showMCTable();
mlmc->showResultsMLMCRun(epsilon);
mlmc->showKurtosisWarning();
mlmc->showExponentResults();
} else if (experimentName == "MLMCOverEpsilon") {
headMLMCOverEpsilon();
mlmc->method(epsilon_lst);
mlmc->showResultsOverEpsilon(epsilon_lst);
}
}
\ No newline at end of file
#ifndef MLMC_MLMCMAIN_H
#define MLMC_MLMCMAIN_H
#include "Elliptic.h"
#include "Mixed.h"
#include "Hybrid.h"
#include "MonteCarlo.h"
#include "assemble/AssembleHeader.h"
#include "problem/StochasticEllipticProblem.h"
#include "problem/StochasticTransportProblem.h"
#include "MultilevelMonteCarlo.h"
#include "discgalerkin/DGDiscretization.h"
void MLMCMain();
class MLMCMain {
public:
double epsilon = 0.1;
int init = 3, plevel = 2, Lmax = 8;
vector<int> l_init = {3, 4, 5}, M_init = {8, 4, 2};
vector<double> epsilon_lst = {0.1, 0.05, 0.01};
string problemName = "", modelName = "";
string fieldName = "", experimentName = "";
Meshes getMeshes(const string &problemName, int plevel, int lmax);
Meshes *meshes = nullptr;
Discretization *disc = nullptr;
StochasticFields *stochFields = nullptr;
StochasticProblem *problem = nullptr;
Assemble *assemble = nullptr;
MatrixGraphs *graphs = nullptr;
Discretization getDiscretization(const string &modelName, Meshes &meshes);
MultilevelMonteCarlo *mlmc = nullptr;
Assemble *getAssemble(const string &problemName, const string &modelName,
Discretization &disc, StochasticProblem *problem);
MLMCMain() {
ReadConfig(Settings, "epsilon", epsilon);
ReadConfig(Settings, "plevel", plevel);
ReadConfig(Settings, "Lmax", Lmax);
ReadConfig(Settings, "init", init);
ReadConfig(Settings, "l_init", l_init, init);
ReadConfig(Settings, "M_init", M_init, init);
ReadConfig(Settings, "epsilon_lst", epsilon_lst, init);
ReadConfig(Settings, "Problem", problemName);
ReadConfig(Settings, "Model", modelName);
ReadConfig(Settings, "StochasticField", fieldName);
ReadConfig(Settings, "Experiment", experimentName);
MatrixGraphs getMatrixGraphs(const string &modelName, Meshes &meshes, Discretization &disc);
initialize();
}
void headConvergenceTest();
~MLMCMain() {
delete meshes, delete disc, delete stochFields;
delete problem, delete assemble, delete graphs;
delete mlmc;
}
void headMLMCExperiment();
void initialize();
void headMLMCOverEpsilon();
Meshes *getMeshes();
Discretization *getDiscretization();
StochasticFields *getStochasticFields();
StochasticProblem *getStochasticProblem();
Assemble *getAssemble();
MatrixGraphs *getMatrixGraphs();
void run();
static void headConvergenceTest() {
mout << endl << "**********************************************************"
<< endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
<< endl << "**********************************************************"
<< endl;
}
static void headMLMCExperiment() {
mout << endl << "**********************************************************"
<< endl << "*** MultilevelMonteCarlo test run ***"
<< endl << "**********************************************************"
<< endl;
}
static void headMLMCOverEpsilon() {
mout << endl << "**********************************************************"
<< endl << "*** MultilevelMonteCarlo over epsilon ***"
<< endl << "**********************************************************"
<< endl;
}
};
#endif //MLMC_MLMCMAIN_H
......@@ -10,6 +10,8 @@ int main(int argv, char **argc) {
SetSearchPath("../mlmc", argv, argc);
DPO dpo(&argv, &argc);
MLMCMain();
auto mlmcMain = new MLMCMain();
mlmcMain->run();
delete mlmcMain;
return 0;
}
......@@ -2,15 +2,33 @@
using namespace std;
MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool baseLevel) {
auto ellipticAssemble = dynamic_cast<EllipticAssemble *>(assemble);
auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
if (ellipticAssemble != nullptr)
return new MonteCarloElliptic(l, dM, baseLevel, meshes, stochFields,
graphs, ellipticAssemble);
if (transportAssemble != nullptr)
return new MonteCarloTransport(l, dM, baseLevel, meshes, stochFields,
graphs, transportAssemble);
return nullptr;
}
void MultilevelMonteCarlo::initializeMapMC() {
for (unsigned long i = 0; i < l_init.size(); i++) {
bool baseLevel = i == 0;
int l = l_init[i], M = M_init[i];
map_mc[l] = getMonteCarlo(l, M, baseLevel, meshes,
stochFields, assemble, graphs);
map_mc[l] = getMonteCarlo(l, M, baseLevel);
}
}
void MultilevelMonteCarlo::reinitializeMLMC() {
levels = {}, numsamples = {};
value = 0.0, cost = 0.0;
alpha = 0.0, beta = 0.0, gamma = 0.0;
}
void MultilevelMonteCarlo::method() {
Date Start;
vout (1) << "Start multilevel Monte Carlo method only for initial values" << endl;
......@@ -52,8 +70,7 @@ void MultilevelMonteCarlo::method(const double eps) {
Exit ("Maximum level has been reached.")
} else {
vout(2) << " Append level " << L << endl;
map_mc[L] = getMonteCarlo(L, 0, false, meshes,
stochFields, assemble, graphs);
map_mc[L] = getMonteCarlo(L, 0, false);
map_mc[L]->var_Y = map_mc[L - 1]->var_Y / pow(2.0, beta);
map_mc[L]->avg_cost = map_mc[L - 1]->avg_cost * pow(2.0, gamma);
updateNumSamples(eps);
......@@ -63,17 +80,30 @@ void MultilevelMonteCarlo::method(const double eps) {
}
}
}
wrapUpResults(value, cost, levels, numsamples);
Vector u((*graphs)[map_mc.rbegin()->first - pLevel]);
u = 0.0;
wrapUpResults(u, value, cost, levels, numsamples);
if (plotting >= 1) {
map_mc.rbegin()->second->plotSolution(u, -1, "");
}
if (plotting >= 2) {
}
vout (1) << "End after " << Date() - Start << endl;
}
void MultilevelMonteCarlo::method(const vector<double> &eps) {
void MultilevelMonteCarlo::method(const vector<double> &eps_lst) {
Date Start;
vout (1) << "Start multilevel Monte Carlo method for several epsilon" << endl;
vout (1) << "Start multilevel Monte Carlo method only for initial values" << endl;
for (auto &e : eps) {
for (auto &eps : eps_lst) {
initializeMapMC();
method(e);
reinitializeMLMC();
method(eps);
levels_vs_eps.push_back(levels);
numsamples_vs_eps.push_back(numsamples);
value_vs_eps.push_back(value);
cost_vs_eps.push_back(cost);
}
vout (1) << "End after " << Date() - Start << endl;
......@@ -91,6 +121,7 @@ void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
vars_Y.push_back(-log2(iter->second->var_Y));
costs.push_back(log2(iter->second->avg_cost));
}
vector<double> v = linear_fit(level_lst, avgs_Y);
alpha = v.front();
v = linear_fit(level_lst, vars_Y);
......@@ -128,24 +159,59 @@ double MultilevelMonteCarlo::estimateNumericalError() {
return maxErr;
}
void MultilevelMonteCarlo::wrapUpResults(double &val, double &mlmc_cost,
vector<int> &level_lst,
vector<int> &samples) {
void MultilevelMonteCarlo::wrapUpResults(Vector &u, double &val, double &mlmc_cost,
vector<int> &level_lst, vector<int> &samples) {
int L = map_mc.rbegin()->first;
for (auto &mc : map_mc) {
val += mc.second->avg_Y;
mlmc_cost += mc.second->sum_cost;
level_lst.push_back(mc.first);
samples.push_back(mc.second->M);
map<int, Vector *> map_u_temp;
for (int l = mc.first; l <= L; l++) {
map_u_temp[l] = new Vector((*graphs)[l - pLevel]);
*map_u_temp[l] = 0.0;
}
if (mc.first != L) {
for (int l = mc.first; l < L; l++) {
if (l == mc.first)
*map_u_temp[l] = mc.second->u_avg_diff;
map_mc[l + 1]->upscaleU(*map_u_temp[l], *map_u_temp[l + 1]);
if (l == L - 1)
u += *map_u_temp[l + 1];
}
} else
u += mc.second->u_avg_diff;
}
}
void MultilevelMonteCarlo::showMCTable() {
void MultilevelMonteCarlo::showMCTableHead() {
mout << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc] V[Qf] kurtosis cost"
<< endl
<< "--------------------------------------------------------------------------------"
<< 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;
}
void MultilevelMonteCarlo::showTableBottom() {
mout << line() << endl;
}
string MultilevelMonteCarlo::line() {
return "----------------------------------------"
"----------------------------------------";
}
void MultilevelMonteCarlo::showMCTable() {
showMCTableHead();
for (auto &mc : map_mc) {
int l = mc.first;
mout << setw(2) << l << setw(6) << mc.second->M
......@@ -154,9 +220,7 @@ void MultilevelMonteCarlo::showMCTable() {
<< setw(12) << mc.second->var_Q << setw(12) << mc.second->kurtosis
<< setw(12) << mc.second->avg_cost << endl;
}
mout
<< "--------------------------------------------------------------------------------"
<< endl;
showTableBottom();
}
void MultilevelMonteCarlo::showKurtosisWarning() {
......@@ -182,14 +246,8 @@ void MultilevelMonteCarlo::showExponentResults() {
}
void MultilevelMonteCarlo::showResultsMLMCRun(double eps) {
showResultsMLMCRunHead();
string level_str, sample_str;
mout << endl
<< " eps value MLMC cost l M"
<< endl
<< "--------------------------------------------------------------------------------"
<< endl;
mout << setw(6) << eps << setw(12) << value << setw(12) << cost;
for (unsigned long i = 0; i < levels.size(); i++) {
if (i != levels.size() - 1) level_str += to_string(levels[i]) + " ";
......@@ -202,28 +260,35 @@ void MultilevelMonteCarlo::showResultsMLMCRun(double eps) {
else sample_str += to_string(numsamples[i]);
}
mout << setw(36) << sample_str << endl;
showTableBottom();
}
void MultilevelMonteCarlo::showResultsOverEpsilon(const vector<double> &eps) {
string level_str, sample_str;
mout << endl
<< " eps value MLMC cost l M"
<< endl
<< "--------------------------------------------------------------------------------"
<< endl;
mout << setw(6) << eps << setw(12) << value << setw(12) << cost;
for (unsigned long i = 0; i < levels.size(); i++) {
if (i != levels.size() - 1) level_str += to_string(levels[i]) + " ";
else level_str += to_string(levels[i]);
}
void MultilevelMonteCarlo::showResultsOverEpsilon(const vector<double> &eps_lst) {
showResultsMLMCRunHead();