Commit 1be73b07 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

stochastic fields on new data structure + refactoring

parent e6596357
......@@ -74,11 +74,11 @@ target_link_libraries(MLMC-M++ MLMC sprng SRC LIB_PS ${SUPERLU} blas lapack fftw
#---------------------------------------------------------------------------------------#
# Test Executables
add_executable(TestCirculantEmbedding mlmc/tests/TestCirculantEmbedding.C)
add_executable(TestRNManager mlmc/tests/TestRNManager.C)
#add_executable(TestCirculantEmbedding mlmc/tests/TestCirculantEmbedding.C)
#add_executable(TestRNManager mlmc/tests/TestRNManager.C)
# Linking
target_link_libraries(TestCirculantEmbedding MultilevelMonteCarlo sprng SRC SRC_Solver
blas lapack ${SUPERLU} fftw3 m ${GTEST_LIB})
target_link_libraries(TestRNManager MultilevelMonteCarlo sprng fftw3 m ${GTEST_LIB})
#target_link_libraries(TestCirculantEmbedding MLMC sprng SRC LIB_PS ${SUPERLU}
# blas lapack fftw3 m ${GTEST_LIB})
#target_link_libraries(TestRNManager MLMC sprng fftw3 m ${GTEST_LIB})
#---------------------------------------------------------------------------------------#
loadconf = elliptic.conf
#loadconf = transport.conf
loadconf = mlmc.conf
\ No newline at end of file
loadconf = stochfield.conf
\ No newline at end of file
Model = LinearFEM
#Model = QuadraticFEM
Model = LagrangeFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = FVFEM # Don't forget to swith on Overlap, sign and penalty
#Model = LinearDGFEM # Don't forget to swith on Overlap, sign and penalty
#Model = QuadraticDGFEM # Don't forget to swith on Overlap, sign and penalty
#Model = DGFEM
#Overlap = dG1
#sign = 1
......@@ -18,24 +15,24 @@ Problem = StochasticLaplace2D
StochasticField = LogNormal
#StochasticField = Gaussian
Experiment = ConvergenceTest
#Experiment = ConvergenceTest
#Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
Experiment = MLMCOverEpsilon
size_init = 4
l_init = 3 4 5 6
M_init = 3 3 3 3
initLevels = 3, 4, 5
initSampleAmount = 12, 6, 3
size_init = 3
l_init = 3 4 5
M_init = 45 15 5
#l_init = 3, 4, 5
#M_init = 45, 15, 5
epsilon = 0.01
size_lst = 2
epsilon_lst = 0.05 0.01
mcOnly = false
epsilon_lst = 0.05, 0.01, 0.005
maxLevel = 8
pLevel = 2
degree = 1
Lmax = 9
plevel = 2
functional = L2
#functional = Energy
#funcitonal = H1
......@@ -53,7 +50,7 @@ enablePython = 1
mean = 1.0
sigma = 1.0
norm_p = 2
lambda1 = 0.15
lambda2 = 0.15
lambda1 = 0.1
lambda2 = 0.1
alpha = 1.8
evtol = 1e-10
\ No newline at end of file
......@@ -9,9 +9,8 @@ set(MLMC_SRC
assemble/MixedEllipticAssemble.C
assemble/HybridEllipticAssemble.C
assemble/DGEllipticAssemble.C
stochfields/RNManager.C
stochfields/StochasticFields.C
stochfields/CirculantEmbedding.C
)
stochastics/RNManager.C
stochastics/CirculantEmbedding.C
problem/StochasticEllipticProblem.h stochastics/SampleGenerator.h stochastics/HybridFluxCalculator.C stochastics/HybridFluxCalculator.h)
add_library(MLMC STATIC ${MLMC_SRC})
#include "MLMCMain.h"
using namespace std;
void MLMCMain::initialize() {
meshes = getMeshes();
disc = getDiscretization();
stochFields = getStochasticFields();
// fieldSample = getFieldSampleVector();
// stochasticField = getStochasticField();
stochasticField = new StochasticField(meshes);
problem = getStochasticProblem();
assemble = getAssemble();
graphs = getMatrixGraphs();
mlmc = new MultilevelMonteCarlo(l_init, M_init, meshes,
stochFields, assemble, graphs);
mlmc = new MultilevelMonteCarlo(initLevels, initSampleAmount, meshes,
stochasticField, assemble, graphs);
}
Meshes *MLMCMain::getMeshes() {
if (problemName.find("1D") != string::npos)
return new Meshes("Line", plevel, Lmax);
return new Meshes("Line", pLevel, maxLevel);
if (problemName.find("2D") != string::npos)
return new Meshes("UnitSquare", plevel, Lmax);
return new Meshes("UnitSquare", pLevel, maxLevel);
Exit("\nMesh name not found\n")
}
Discretization *MLMCMain::getDiscretization(string _modelName_) {
if (_modelName_.empty()) _modelName_ = modelName;
if (_modelName_ == "LinearFEM")
Discretization *MLMCMain::getDiscretization() {
if (modelName == "LagrangeFEM") {
if (degree == 1)
return new Discretization("linear", 1, meshes->dim());
if (_modelName_ == "QuadraticFEM")
return new Discretization("quadratic", 1, meshes->dim());
if (_modelName_ == "MixedFEM" || _modelName_ == "HybridFEM")
if (degree == 2)
return new Discretization("serendipity", 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));
if (modelName == "DGFEM")
return new DGDiscretization(new DGDoF(degree));
if (modelName == "EGFEM")
return new EGDiscretization(new EGDoF(degree));
if (modelName == "LinearWC")
return nullptr;
if (modelName == "QuadraticWC")
return nullptr;
if (modelName == "LinearDPG")
return nullptr;
if (modelName == "QuadraticDPG")
return nullptr;
if (modelName == "DGTransport")
return new DGDiscretization(new DGDoF(degree));
Exit("\nDiscretization not implemented yet\n")
}
StochasticFields *MLMCMain::getStochasticFields() {
return new StochasticFields(fieldName, *meshes);
}
StochasticProblem *MLMCMain::getStochasticProblem(string _problemName_) {
if (_problemName_.empty()) _problemName_ = problemName;
if (_problemName_ == "StochasticLaplace1D")
return new StochasticLaplace1D(*stochFields);
if (_problemName_ == "StochasticLaplace2D")
return new StochasticLaplace2D(*stochFields);
if (_problemName_ == "StochasticTransport1D")
return nullptr;
if (_problemName_ == "StochasticTransport2D") {
// Discretization *_disc_ = getDiscretization("HybridFEM");
// StochasticProblem *_problem_ = getStochasticProblem(
// "StochasticLaplace2D");
//
// Assemble *hybridAssemble = getAssemble("HybridFEM",
// "StochasticLaplace2D",
// _disc_, _problem_);
// Vector normalFlux(graphs->fine());
// return new StochasticPollution2D(hybridAssemble, normalFlux, stochFields);
Vector *MLMCMain::getFieldSampleVector() {
if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D") {
MatrixGraphs cellMatrixGraphs(*meshes, dof("cell", 1));
return new Vector(cellMatrixGraphs.fine());
}
if (problemName == "StochasticPollution1D" || "StochasticPollution2D") {
MatrixGraphs cellMatrixGraphs(*meshes, dof("cell", 3));
return new Vector(cellMatrixGraphs.fine());
}
}
StochasticField *MLMCMain::getStochasticField() {
if (problemName == "StochasticLaplace1D")
return new StochasticField(meshes);
if (problemName == "StochasticLaplace1D")
return new StochasticField(meshes);
if (problemName == "StochasticPollution1D")
return new StochasticVectorField(meshes);
if (problemName == "StochasticPollution2D")
return new StochasticField(meshes);
}
StochasticProblem *MLMCMain::getStochasticProblem() {
if (problemName == "StochasticLaplace1D")
return new StochasticLaplace1D();
if (problemName == "StochasticLaplace2D")
return new StochasticLaplace2D();
if (problemName == "StochasticPollution1D")
return new StochasticPollution1D();
if (problemName == "StochasticPollution2D")
return new StochasticPollution2D();
Exit("\nStochastic problem not implemented yet\n")
return nullptr;
}
Assemble *MLMCMain::getAssemble(string _modelName_, string _problemName_,
Discretization *_disc_,
StochasticProblem *_problem_) {
if (_problemName_.empty()) _problemName_ = problemName;
if (_modelName_.empty()) _modelName_ = modelName;
if (_disc_ == nullptr) _disc_ = disc;
if (_problem_ == nullptr) _problem_ = problem;
if (_problemName_.find("Laplace") != string::npos) {
if (_modelName_ == "LinearFEM" || _modelName_ == "QuadraticFEM")
return new EllipticAssemble(_disc_, _problem_);
if (_modelName_ == "MixedFEM")
return new MixedEllipticAssemble(_disc_, _problem_);
if (_modelName_ == "HybridFEM")
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_ == "FVFEM")
return nullptr;
if (_modelName_ == "LinearDGFEM")
return nullptr;
if (_modelName_ == "QuadraticDGFEM")
return nullptr;
Assemble *MLMCMain::getAssemble() {
if (problemName.find("Laplace") != string::npos) {
auto ellipticProblem = dynamic_cast<StochasticEllipticProblem *>(problem);
if (modelName == "LagrangeFEM")
return new EllipticAssemble(disc, ellipticProblem);
if (modelName == "MixedFEM")
return new MixedEllipticAssemble(disc, ellipticProblem);
if (modelName == "HybridFEM")
return new HybridEllipticAssemble(disc, ellipticProblem);
if (modelName == "DGFEM")
return new DGEllipticAssemble(disc, ellipticProblem);
} else if (problemName.find("Pollution") != string::npos) {
auto transportProblem = dynamic_cast<StochasticTransportProblem *>(problem);
Plot plot(meshes->fine(), meshes->dim(), 2);
return new DGTransportAssemble(disc, transportProblem, &plot);
}
Exit("\nAssembling for this problem not implemented yet\n")
return nullptr;
}
MatrixGraphs *MLMCMain::getMatrixGraphs(string _modelName_,
Meshes *_meshes_,
Discretization *_disc_) {
if (_modelName_.empty()) _modelName_ = modelName;
if (_meshes_ == nullptr) _meshes_ = meshes;
if (_disc_ == nullptr) _disc_ = disc;
if (_modelName_ == "LinearFEM" || _modelName_ == "QuadraticFEM" ||
_modelName_ == "MixedFEM")
return new MatrixGraphs(*_meshes_, *_disc_);
if (_modelName_ == "HybridFEM")
return new MatrixGraphs(*_meshes_, dof("face", 1));
if (_modelName_ == "FVFEM" || _modelName_ == "LinearDGFEM" ||
_modelName_ == "QuadraticDGFEM")
return new CellMatrixGraphs(*_meshes_, *_disc_);
Exit("\nMatrix graph Not implemented yet\n")
return nullptr;
MatrixGraphs *MLMCMain::getMatrixGraphs() {
if (modelName == "LagrangeFEM" || modelName == "MixedFEM")
return new MatrixGraphs(*meshes, *disc);
if (modelName == "HybridFEM")
return new MatrixGraphs(*meshes, dof("face", 1));
if (modelName == "DGFEM")
return new CellMatrixGraphs(*meshes, *disc);
Exit("\nMatrix graph not implemented yet\n")
}
void MLMCMain::run() {
printParameters();
//Todo ConvergenceTest bool mcOnly
if (experimentName == "ConvergenceTest") {
headConvergenceTest();
mlmc->method();
mlmc->showMCTable();
mlmc->showKurtosisWarning();
mlmc->showExponentResults();
mlmc->Method();
mlmc->ShowMCTable();
mlmc->ShowKurtosisWarning();
mlmc->ShowExponentResults();
if (enablePython == 1) {
system("python3 ../tools/plot_statistics.py ConvergenceTest");
}
} else if (experimentName == "MLMCExperiment") {
headMLMCExperiment();
mlmc->method(epsilon);
mlmc->showMCTable();
mlmc->showResultsMLMCRun(epsilon);
mlmc->showKurtosisWarning();
mlmc->showExponentResults();
mlmc->Method(epsilon);
mlmc->ShowMCTable();
mlmc->ShowResultsMLMCRun(epsilon);
mlmc->ShowKurtosisWarning();
mlmc->ShowExponentResults();
if (enablePython == 1 && PPM->master()) {
system("python3 ../tools/plot_statistics.py MLMCExperiment");
system("python3 ../tools/plot_mlmc.py MLMCExperiment");
}
} else if (experimentName == "MLMCOverEpsilon") {
headMLMCOverEpsilon();
mlmc->method(epsilon_lst);
mlmc->showMCTable();
mlmc->showResultsOverEpsilon(epsilon_lst);
mlmc->Method(epsilonList);
mlmc->ShowMCTable();
mlmc->ShowResultsOverEpsilon(epsilonList);
if (enablePython == 1 && PPM->master()) {
system("python3 ../tools/plot_statistics.py MLMCOverEpsilon");
system("python3 ../tools/plot_mlmc.py MLMCOverEpsilon");
......@@ -156,13 +159,13 @@ void MLMCMain::printParameters() {
mout << "Model \t=\t" << modelName << endl;
mout << "Experiment\t=\t" << experimentName << endl;
mout << "StochField\t=\t" << fieldName << endl;
mout << "l_init \t=\t" << vec2str(l_init) << endl;
mout << "M_init \t=\t" << vec2str(M_init) << endl;
mout << "initLevels\t=\t" << vec2str(initLevels) << endl;
mout << "initSample\t=\t" << vec2str(initSampleAmount) << endl;
if (experimentName == "MLMCExperiment")
mout << "eps \t=\t" << epsilon << endl;
if (experimentName == "MLMCOverEpsilon")
mout << "eps_lst \t=\t" << vec2str(epsilon_lst) << endl;
mout << "Lmax \t=\t" << Lmax << endl;
mout << "plevel \t=\t" << plevel << endl;
mout << "eps_lst \t=\t" << vec2str(epsilonList) << endl;
mout << "maxLevel \t=\t" << maxLevel << endl;
mout << "pLevel \t=\t" << pLevel << endl;
mout << "enablePy \t=\t" << enablePython << endl;
}
......@@ -5,47 +5,55 @@
#include "problem/StochasticEllipticProblem.h"
#include "problem/StochasticTransportProblem.h"
#include "MultilevelMonteCarlo.h"
#include "discgalerkin/DGDiscretization.h"
#include "discretization/DGDiscretization.h"
#include "stochastics/StochasticField.h"
class MLMCMain {
public:
int enablePython = 1;
double epsilon = 0.1;
int size_init = 3, plevel = 2, Lmax = 8, size_lst = 3;
vector<int> l_init = {3, 4, 5}, M_init = {12, 6, 3};
vector<double> epsilon_lst = {0.1, 0.05, 0.01};
string problemName = "", modelName = "";
string fieldName = "", experimentName = "";
int enablePython = 1;
double epsilon = 0.1;
int pLevel = 0, maxLevel = 8, degree = 0;
vector<int> initLevels = {3, 4, 5}, initSampleAmount = {12, 6, 3};
vector<double> epsilonList = {0.1, 0.05, 0.01};
Meshes *meshes = nullptr;
Discretization *disc = nullptr;
StochasticFields *stochFields = nullptr;
StochasticField *stochasticField = nullptr;
StochasticProblem *problem = nullptr;
Assemble *assemble = nullptr;
MatrixGraphs *graphs = nullptr;
Plot *plot = nullptr;
MultilevelMonteCarlo *mlmc = nullptr;
MLMCMain() {
ReadConfig(Settings, "enablePython", enablePython);
ReadConfig(Settings, "epsilon", epsilon);
ReadConfig(Settings, "plevel", plevel);
ReadConfig(Settings, "Lmax", Lmax);
ReadConfig(Settings, "size_init", size_init);
ReadConfig(Settings, "l_init", l_init, size_init);
ReadConfig(Settings, "M_init", M_init, size_init);
ReadConfig(Settings, "size_lst", size_lst);
ReadConfig(Settings, "epsilon_lst", epsilon_lst, size_lst);
ReadConfig(Settings, "Problem", problemName);
ReadConfig(Settings, "Model", modelName);
ReadConfig(Settings, "StochasticField", fieldName);
ReadConfig(Settings, "Experiment", experimentName);
config.get("enablePython", enablePython);
config.get("epsilon", epsilon);
config.get("pLevel", pLevel);
config.get("maxLevel", maxLevel);
config.get("degree", degree);
config.get("initLevels", initLevels);
config.get("initSampleAmount", initSampleAmount);
config.get("epsilonList", epsilonList);
config.get("Problem", problemName);
config.get("Model", modelName);
config.get("StochasticField", fieldName);
config.get("Experiment", experimentName);
if (initSampleAmount.size() != initLevels.size()) Exit(
"\ninitLevels and initSampleAmount have to be of the same size.\n")
if (maxLevel < initLevels.back()) Exit(
"\nLast element of initLevels has to be smaller or equal to maxLevel.\n")
initialize();
}
~MLMCMain() {
delete meshes, delete disc, delete stochFields;
delete meshes, delete disc, delete stochasticField;
delete problem, delete assemble, delete graphs;
delete mlmc;
}
......@@ -54,19 +62,17 @@ public:
Meshes *getMeshes();
Discretization *getDiscretization(string _modelName_ = "");
Discretization *getDiscretization();
Vector *getFieldSampleVector();
StochasticFields *getStochasticFields();
StochasticField *getStochasticField();
StochasticProblem *getStochasticProblem(string _problemName_ = "");
StochasticProblem *getStochasticProblem();
Assemble *getAssemble(string _modelName_ = "", string _problemName_ = "",
Discretization *_disc_ = nullptr,
StochasticProblem *_problem_ = nullptr);
Assemble *getAssemble();
MatrixGraphs *getMatrixGraphs(string _modelName_ = "",
Meshes *_meshes_ = nullptr,
Discretization *_disc_ = nullptr);
MatrixGraphs *getMatrixGraphs();
void run();
......
#include "MultilevelMonteCarlo.h"
using namespace std;
MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool baseLevel) {
MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool onlyFineLevel) {
auto ellipticAssemble = dynamic_cast<EllipticAssemble *>(assemble);
// auto transportAssemble = dynamic_cast<DGTransportAssemble *>(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 new MonteCarloElliptic(l, dM, onlyFineLevel, meshes,
stochasticField, ellipticAssemble);
if (transportAssemble != nullptr)
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes,
stochasticField, transportAssemble);
return nullptr;
}
void MultilevelMonteCarlo::initializeMapMC() {
clearMapMC();
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);
void MultilevelMonteCarlo::initializeMapMonteCarlo() {
clearMapMonteCarlo();
for (unsigned long i = 0; i < initLevels.size(); i++) {
bool onlyFineLevel = (i == 0) || (mcOnly);
int l = initLevels[i], M = initSampleAmount[i];
mapMonteCarlo[l] = getMonteCarlo(l, M, onlyFineLevel);
}
}
void MultilevelMonteCarlo::initializeValuesMLMC() {
void MultilevelMonteCarlo::initializeValues() {
levels = {}, numsamples = {};
value = 0.0, cost = 0.0;
alpha = 0.0, beta = 0.0, gamma = 0.0;
}
void MultilevelMonteCarlo::clearMapMC() {
for (auto &iter: map_mc) {
void MultilevelMonteCarlo::clearMapMonteCarlo() {
for (auto &iter: mapMonteCarlo) {
delete iter.second;
}
map_mc.clear();
mapMonteCarlo.clear();
}
void MultilevelMonteCarlo::method() {
void MultilevelMonteCarlo::Method() {
logger->startMethodMSG();
for (auto &mc : map_mc)
mc.second->method();
for (auto &mc : mapMonteCarlo)
mc.second->Method();
logger->endMethodMSG();
}
void MultilevelMonteCarlo::method(const double eps) {
void MultilevelMonteCarlo::Method(const double eps) {
logger->startMethodMSG();
logger->logMSGV1("eps: " + to_string(eps));
bool converged = false;
while (!converged) {
for (auto &mc : map_mc)
for (auto &mc : mapMonteCarlo)
if (mc.second->dM > 0)
mc.second->method();
mc.second->Method();
estimateExponents();
updateNumSamples(eps);
updateSampleAmount(eps);
if (checkForNewLevel()) {
if (estimateNumericalError() > eps / sqrt(2.0))
appendLevel(eps);
......@@ -61,23 +61,23 @@ void MultilevelMonteCarlo::method(const double eps) {
converged = true;
}
}
Vector u((*graphs)[map_mc.rbegin()->first - pLevel]);
Vector u((*graphs)[mapMonteCarlo.rbegin()->first - pLevel]);
u = 0.0;
wrapUpResults(u, value, cost, levels, numsamples);
logger->endMethodMSG();
}
void MultilevelMonteCarlo::method(const vector<double> &eps_lst) {
void MultilevelMonteCarlo::Method(const vector<double> &eps_lst) {
logger->startMethodMSG();
// logger->increase_indent();
for (auto &eps : eps_lst) {
initializeMapMC();
initializeValuesMLMC();
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);
initializeMapMonteCarlo();
initializeValues();
Method(eps);
levelsOverEpsilon.push_back(levels);
sampleAmountOverEpsilon.push_back(numsamples);
valueOverEpsilon.push_back(value);
costOverEpsilon.push_back(cost);
}
// logger->decrease_indent();
logger->endMethodMSG();
......@@ -85,15 +85,15 @@ void MultilevelMonteCarlo::method(const vector<double> &eps_lst) {
void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
vector<double> level_lst, avgs_Y, vars_Y, costs;
auto iter = map_mc.begin();
auto iter = mapMonteCarlo.begin();
if (excludeBaseLevel)
iter++;
for (; iter != map_mc.end(); iter++) {
for (; iter != mapMonteCarlo.end(); iter++) {
level_lst.push_back((double) iter->first);
avgs_Y.push_back(-log2(iter->second->avg_Y));
vars_Y.push_back(-log2(iter->second->var_Y));
costs.push_back(log2(iter->second->avg_cost));
avgs_Y.push_back(-log2(iter->second->avgY));
vars_Y.push_back(-log2(iter->second->varY));