Commit 41a05672 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

worked on framework

parent 48e4820a
......@@ -4,23 +4,22 @@ Model = LinearFEM
#Model = HybridFEM
#Model = DGFEM
Problem = StochasticLaplace1D
#Problem = StochasticLaplace2D
#Problem = StochasticLaplace1D
Problem = StochasticLaplace2D
StochasticField = LogNormal
#StochasticField = Gaussian
#Experiment = ConvergenceTest
Experiment = MLMCExperiment
Experiment = ConvergenceTest
#Experiment = MLMCExperiment
#Experiment = MLMCOverEpsilon
#l_init = 3 4 5 6 7
l_init = 3 4 5
#M_init = 20 20 20 20 20
M_init = 40 20 10
init = 2
l_init = 4 5
M_init = 1 1
Lmax = 8
epsilon = 0.1
epsilon = 0.001
epsilon_lst = 0.01 0.005 0.001
plevel = 2
......@@ -29,12 +28,10 @@ functional = L2
#funcitonal = H1
#functional = Outflow
MCVerbose = 0
MLMCVerbose = 0
MCPlotting = 3
MCVerbose = 1
MLMCVerbose = 2
MCPlotting = 1
MLMCPlotting = 0
vtkplot = 1
gnuplot = 1
......@@ -7,6 +7,19 @@
using namespace std;
/*
* Todo:
* - Use Information from Mesh to be able to initialzie the right domain for
* circulant embedding
* - Remove hack in Stochastic fields for problem size
* -
*
*
*/
class CirculantEmbedding {
/* Circulant Embedding
......@@ -57,7 +70,7 @@ public:
int N_embedded;
int N;
double domain_start = 0; //TODO
double domain_start = 0;
double domain_end = 1;
explicit CirculantEmbedding1D(int N, int streamnum = 0, int nstreams = 1, int gtype = 0)
......@@ -91,7 +104,7 @@ public:
int N_embedded[2]{};
int N[2]{};
double domain1_start = 0; //TODO
double domain1_start = 0;
double domain1_end = 1;
double domain2_start = 0;
double domain2_end = 1;
......
......@@ -426,7 +426,7 @@ public:
const char *name, int l) {
plot.vertexdata(u);
plot.vtk_vertexdata(name);
plot.gnu_vertexdata(name);
// plot.gnu_vertexdata(name);
}
virtual void plotUex(Plot &plot, const Vector &u, Meshes &meshes,
......@@ -436,7 +436,7 @@ public:
setExactSolution(u_ex);
plot.vertexdata(u_ex);
plot.vtk_vertexdata(name);
plot.gnu_vertexdata(name);
// plot.gnu_vertexdata(name);
}
virtual void plotFlux(Plot &plot, const Vector &u, Meshes &meshes,
......@@ -447,7 +447,7 @@ public:
setFlux(u, flux);
plot.celldata(flux, 3);
plot.vtk_celldata(name);
plot.gnu_celldata(name);
// plot.gnu_celldata(name);
}
virtual void plotPermeability(Plot &plot, const Vector &u, Meshes &meshes,
......@@ -461,7 +461,7 @@ public:
}
plot.celldata(permeability);
plot.vtk_celldata(name);
plot.gnu_celldata(name);
// plot.gnu_celldata(name);
}
static const char *buildString(string name, int l, int m, const string &suffix) {
......
......@@ -2,7 +2,7 @@
void MLMCMain() {
double epsilon;
int plevel, Lmax;
int init, plevel, Lmax;
vector<int> l_init, M_init;
vector<double> epsilon_lst;
string problemName, modelName, fieldName, experimentName;
......@@ -10,9 +10,10 @@ void MLMCMain() {
ReadConfig(Settings, "epsilon", epsilon);
ReadConfig(Settings, "plevel", plevel);
ReadConfig(Settings, "Lmax", Lmax);
ReadConfig(Settings, "l_init", l_init, 3);
ReadConfig(Settings, "M_init", M_init, 3);
ReadConfig(Settings, "epsilon_lst", epsilon_lst, 3);
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);
......@@ -44,6 +45,7 @@ void MLMCMain() {
} else if (experimentName == "MLMCOverEpsilon") {
headMLMCOverEpsilon();
mlmc.method(epsilon_lst);
mlmc.showResultsOverEpsilon(epsilon_lst);
}
}
......
......@@ -31,28 +31,41 @@ void MonteCarlo::updateStatistic() {
3.0 * avg_Y * avg_Y * avg_Y * avg_Y) / var_Y / var_Y;
}
void MonteCarlo::upscaleU() {
// TransferMatrix transferMatrix(TransferMatrixGraph(graphs[l - plevel],
// graphs[l - plevel - 1]));
// assemble.AssembleTransfer(transferMatrix);
}
void MonteCarlo::updateUSum() {
}
void MonteCarlo::updateUAvg() {
}
class MonteCarloElliptic : public MonteCarlo {
public:
EllipticAssemble &assemble;
Transfer *transfer;
MonteCarloElliptic(int l, int dM, bool baseLevel, Meshes &meshes,
StochasticFields &stochFields, MatrixGraphs &graphs,
EllipticAssemble &assemble) : assemble(assemble),
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs) {}
EllipticAssemble &assemble) :
assemble(assemble),
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs) {
transfer = new MatrixTransfer(assemble);
transfer->Construct(graphs[l - plevel], graphs[l - plevel - 1]);
}
~MonteCarloElliptic() {
transfer->Destruct();
delete transfer;
}
void upscaleU(const Vector &uc, Vector &uc_up) override {
uc_up = (*transfer) * uc;
}
void updateUSum(const Vector &uc_up, const Vector &uf) override {
u_sum += uf;
u_sum_diff += uf - uc_up;
}
void updateUAvg() override {
Scalar M_inverted = 1.0 / M;
u_avg = M_inverted * u_sum;
u_avg_diff = M_inverted * u_sum_diff;
}
void method() override {
Date Start;
......@@ -77,50 +90,58 @@ public:
}
}
double Qf, Qc;
double costf, costc;
double Qf, Qc, costf, costc;;
Qf = 0.0, Qc = 0.0, costf = 0.0, costc = 0.0;
u_f = 0.0, u_c = 0.0;
Vector uf(u_sum), uc_up(u_sum), uc(graphs[l - plevel - 1]);
uf = 0.0, uc = 0.0, uc_up = 0.0;
Date Start_solving;
vout(2) << " Solving PDE on level " << l;
solvePDE(Qf, costf, u_f);
solvePDE(Qf, costf, uf);
vout(2) << " took " << Start_solving - Date() << endl;
if (!baseLevel) {
Date Start_solving_c;
vout(2) << " Solving PDE on level " << l - 1;
stochFields.setFineField(false);
solvePDE(Qc, costc, u_c);
solvePDE(Qc, costc, uc);
stochFields.setFineField(true);
vout(2) << " took " << Start_solving_c - Date() << endl;
} else {
Qc = 0.0;
costc = 0.0;
u_c = 0.0;
upscaleU();
uc = 0.0;
}
double sums[7] = {costf, Qf - Qc, pow((Qf - Qc), 2),
pow((Qf - Qc), 3), pow((Qf - Qc), 4),
Qf, pow(Qf, 2)};
updateSum(sums);
updateUSum();
upscaleU(uc, uc_up);
updateUSum(uc_up, uf);
if (plotting >= 2)
assemble.plotSolution(plot_f, u_f, meshes, l, m);
assemble.plotSolution(plot, uf, meshes, l, m);
if (plotting == 3 && !baseLevel) {
stochFields.setFineField(false);
assemble.plotSolution(plot_c, u_c, meshes, l - 1, m, "ip");
assemble.plotSolution(plot, uc_up, meshes, l, m, "up");
stochFields.setFineField(true);
}
}
updateAvg();
updateUAvg();
updateStatistic();
if (plotting >= 1)
// assemble.plotSolution(plot_f, u_avg, meshes, l, -1, "avg");
if (plotting >= 1) {
mout << OUT(u_avg_diff) << endl;
assemble.plotSolution(plot, u_avg_diff, meshes, l,
-1, "avg_diff");
assemble.plotSolution(plot, u_avg, meshes, l,
-1, "avg");
}
vout (1) << " End after " << Date() - Start << ": |E(Y)|=" << avg_Y
<< " V(Y)=" << var_Y << endl;
......@@ -149,29 +170,7 @@ MonteCarlo *getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
if (ellipticAssemble != nullptr)
return new MonteCarloElliptic(l, dM, baseLevel, meshes,
stochFields, graphs, (*ellipticAssemble));
if(transportAssemble != nullptr)
if (transportAssemble != nullptr)
return nullptr;
return nullptr;
}
void estimateExponents(const map<int, MonteCarlo *> &df_mc,
double &alpha, double &beta, double &gamma,
bool excludeBaseLevel) {
vector<double> levels, avgs_Y, vars_Y, costs;
auto iter = df_mc.begin();
if (excludeBaseLevel)
iter++;
for (; iter != df_mc.end(); iter++) {
levels.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));
}
vector<double> v = linear_fit(levels, avgs_Y);
alpha = v.front();
v = linear_fit(levels, vars_Y);
beta = v.front();
v = linear_fit(levels, costs);
gamma = v.front();
}
......@@ -46,9 +46,9 @@ public:
StochasticFields &stochFields;
MatrixGraphs &graphs;
Vector u_sum, u_avg, u_f, u_c;
Vector u_sum, u_avg, u_sum_diff, u_avg_diff;
Plot plot_c, plot_f;
Plot plot;
int M = 0, dM = 0;
......@@ -64,20 +64,20 @@ public:
MonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
StochasticFields &stochFields, MatrixGraphs &graphs) :
l(l), dM(dM), baseLevel(baseLevel), meshes(meshes), plevel(meshes.pLevel()),
l(l), dM(dM), baseLevel(baseLevel),
meshes(meshes), plevel(meshes.pLevel()),
stochFields(stochFields), graphs(graphs),
u_sum(graphs[l - plevel]), u_avg(u_sum),
u_f(u_sum), u_c(graphs[l - plevel - 1]),
plot_c(meshes[l - 1], meshes.dim(), 2),
plot_f(meshes[l], meshes.dim(), 2) {
u_sum_diff(u_sum), u_avg_diff(u_sum),
plot(meshes[l], meshes.dim(), 2) {
ReadConfig(Settings, "MCVerbose", verbose);
ReadConfig(Settings, "MCPlotting", plotting);
ReadConfig(Settings, "functional", functional);
}
u_sum = 0.0, u_avg = 0.0, u_sum_diff = 0.0, u_avg_diff = 0.0;
virtual ~MonteCarlo() = default;;
}
void updateSum(const double *sums);
......@@ -85,24 +85,20 @@ public:
void updateStatistic();
void upscaleU();
virtual void upscaleU(const Vector &uc, Vector &uc_up) = 0;
void updateUSum();
virtual void updateUSum(const Vector &uc_up, const Vector &uf) = 0;
void updateUAvg();
virtual void updateUAvg() = 0;
virtual void method() = 0;
virtual void solvePDE(double &Q, double &cost, Vector &u) {};
};
MonteCarlo *
getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
StochasticFields &stochFields, Assemble *assemble,
MatrixGraphs &graphs);
MonteCarlo *getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
StochasticFields &stochFields, Assemble *assemble,
MatrixGraphs &graphs);
void estimateExponents(const map<int, MonteCarlo *> &df_mc,
double &alpha, double &beta, double &gamma,
bool excludeBaseLevel = true);
#endif //MLMC_MC_HPP
\ No newline at end of file
......@@ -2,6 +2,15 @@
using namespace std;
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);
}
}
void MultilevelMonteCarlo::method() {
Date Start;
vout (1) << "Start multilevel Monte Carlo method only for initial values" << endl;
......@@ -22,12 +31,11 @@ void MultilevelMonteCarlo::method(const double eps) {
if (mc.second->dM > 0)
mc.second->method();
estimateExponents(map_mc, alpha, beta, gamma);
estimateExponents();
vout (2) << " Estimated exponents: " << "alpha=" << alpha
<< " beta=" << beta << " gamma=" << gamma << endl;
updateNumSamples(eps);
bool check_for_new_level = true;
for (auto &mc : map_mc) {
if (mc.second->dM > 0) {
......@@ -45,8 +53,7 @@ void MultilevelMonteCarlo::method(const double eps) {
} else {
vout(2) << " Append level " << L << endl;
map_mc[L] = getMonteCarlo(L, 0, false, meshes,
stochFields,
assemble, graphs);
stochFields, assemble, graphs);
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);
......@@ -56,26 +63,42 @@ void MultilevelMonteCarlo::method(const double eps) {
}
}
}
for (auto &mc : map_mc) {
value += mc.second->avg_Y;
cost += mc.second->sum_cost;
levels.push_back(mc.first);
numsamples.push_back(mc.second->M);
}
wrapUpResults(value, cost, levels, numsamples);
vout (1) << "End after " << Date() - Start << endl;
}
void MultilevelMonteCarlo::method(const vector<double> &eps) {
Date Start;
vout (1) << "Start multilevel Monte Carlo method for several eps" << endl;
vout (1) << "Start multilevel Monte Carlo method for several epsilon" << endl;
for (auto &e : eps)
for (auto &e : eps) {
initializeMapMC();
method(e);
}
vout (1) << "End after " << Date() - Start << endl;
}
void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
vector<double> level_lst, avgs_Y, vars_Y, costs;
auto iter = map_mc.begin();
if (excludeBaseLevel)
iter++;
for (; iter != map_mc.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));
}
vector<double> v = linear_fit(level_lst, avgs_Y);
alpha = v.front();
v = linear_fit(level_lst, vars_Y);
beta = v.front();
v = linear_fit(level_lst, costs);
gamma = v.front();
}
void MultilevelMonteCarlo::updateNumSamples(const double &eps) {
int M_optimal;
double factor = 0.0;
......@@ -84,6 +107,7 @@ void MultilevelMonteCarlo::updateNumSamples(const double &eps) {
for (auto &mc : map_mc) {
M_optimal = (int) (ceil(2 * pow(eps, -2) * factor *
sqrt(mc.second->var_Y / mc.second->avg_cost)));
if (M_optimal == 1) M_optimal++; // Hack
map_mc[mc.first]->dM = M_optimal - map_mc[mc.first]->M;
}
vout(2) << " Missing samples:";
......@@ -104,6 +128,17 @@ double MultilevelMonteCarlo::estimateNumericalError() {
return maxErr;
}
void MultilevelMonteCarlo::wrapUpResults(double &val, double &mlmc_cost,
vector<int> &level_lst,
vector<int> &samples) {
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);
}
}
void MultilevelMonteCarlo::showMCTable() {
mout << endl
<< " l M E[Qf-Qc] E[Qf] V[Qf-Qc] V[Qf] kurtosis cost"
......@@ -136,7 +171,7 @@ void MultilevelMonteCarlo::showKurtosisWarning() {
}
void MultilevelMonteCarlo::showExponentResults() {
estimateExponents(map_mc, alpha, beta, gamma, true);
estimateExponents();
mout << endl << "alpha = " << setw(8) << alpha
<< " (exponent for the weak convergence of the FEM)"
<< endl << " beta = " << setw(8) << beta
......@@ -156,7 +191,7 @@ void MultilevelMonteCarlo::showResultsMLMCRun(double eps) {
<< endl;
mout << setw(6) << eps << setw(12) << value << setw(12) << cost;
for (unsigned long i = 0 ; i < levels.size(); i++) {
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]);
}
......@@ -169,5 +204,26 @@ void MultilevelMonteCarlo::showResultsMLMCRun(double eps) {
mout << setw(36) << sample_str << endl;
}
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]);
}
mout << setw(14) << level_str;
for (unsigned long i = 0; i < numsamples.size(); i++) {
if (i != numsamples.size() - 1) sample_str += to_string(numsamples[i]) + " ";
else sample_str += to_string(numsamples[i]);
}
mout << setw(36) << sample_str << endl;
}
......@@ -10,6 +10,7 @@ public:
int verbose = 0, plotting = 0;
int Lmax;
vector<int> &l_init, &M_init;
map<int, MonteCarlo *> map_mc;
Meshes &meshes;
......@@ -21,38 +22,41 @@ public:
double value = 0.0, cost = 0.0;
double alpha = 0.0, beta = 0.0, gamma = 0.0;
MultilevelMonteCarlo(vector<int> l_init, vector<int> M_init, Meshes &meshes,
StochasticFields &stochFields, Assemble *assemble,
MatrixGraphs &graphs) :
stochFields(stochFields), assemble(assemble), graphs(graphs),
meshes(meshes), Lmax(meshes.Level()) {
MultilevelMonteCarlo(vector<int> &l_init, vector<int> &M_init,
Meshes &meshes, StochasticFields &stochFields,
Assemble *assemble, MatrixGraphs &graphs) :
l_init(l_init), M_init(M_init), stochFields(stochFields), assemble(assemble),
graphs(graphs), meshes(meshes), Lmax(meshes.Level()) {
ReadConfig(Settings, "MLMCVerbose", verbose);
ReadConfig(Settings, "MLMCPlotting", plotting);
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);
}
initializeMapMC();
}
void initializeMapMC();
~MultilevelMonteCarlo() {
for (auto &iter: map_mc) {
delete iter.second;
}
}
void estimateExponents(bool excludeBaseLevel = true);
void updateNumSamples(const double &eps);
double estimateNumericalError();
void wrapUpResults(double &val, double &mlmc_cost,
vector<int> &level_lst, vector<int> &samples);
void method();
void method(double eps);
void method(const vector<double>& eps);
void method(const vector<double> &eps);
void showMCTable();
......@@ -62,7 +66,7 @@ public:
void showResultsMLMCRun(double eps);
void showResultsOverEpsilon();
void showResultsOverEpsilon(const vector<double> &eps);