new structure in main

mlmc/src/MLMCMain.C

 #include "MLMCMain.h"
 
-void ConvergenceTest() {
-    int l_mc, L_mc, M;
-    ReadConfig(Settings, "l_mc", l_mc);
-    ReadConfig(Settings, "L_mc", L_mc);
-    ReadConfig(Settings, "M", M);
-
-    mout << endl << "**********************************************************"
-         << endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
-         << endl << "**********************************************************" << endl;
+void MLMCMain() {
+    double epsilon;
+    int plevel, Lmax;
+    vector<int> l_init, M_init;
+    vector<double> epsilon_lst;
+    string problemName, modelName, fieldName, experimentName;
 
-    mout << endl
-         << " l    E(Qf-Qc)       E(Qf)    V(Qf-Qc)       V(Qf)    kurtosis        cost"
-         << endl
-         << "--------------------------------------------------------------------------"
-         << endl;
+    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, "Problem", problemName);
+    ReadConfig(Settings, "Model", modelName);
+    ReadConfig(Settings, "StochasticField", fieldName);
+    ReadConfig(Settings, "Experiment", experimentName);
 
-    Meshes meshes(L_mc);
-    Discretization disc(1, meshes.dim());
-    StochasticFields stochFields(meshes);
-    StochasticProblem *problem = getStochasticProblem(stochFields,
-                                                      "StochLaplace2D");
-    EllipticAssemble ellipticAssemble(disc, *problem);
-    MatrixGraphs graphs(meshes, disc);
-    Vector u(graphs.fine());
-    u = 0;
+    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);
 
-    map<int, MonteCarlo *> map_mc;
-    for (int l = l_mc; l <= L_mc; l++) {
-        bool base_l = l == l_mc;
-        map_mc[l] = getMonteCarlo(l, M, base_l, stochFields, ellipticAssemble,
-                                  graphs,
-                                  "MonteCarloElliptic");
-        map_mc[l]->method();
-        mout << setw(2) << l << setw(12) << map_mc[l]->avg_Y
-             << setw(12) << map_mc[l]->avg_Q << setw(12) << map_mc[l]->var_Y
-             << setw(12) << map_mc[l]->var_Q << setw(12) << map_mc[l]->kurtosis
-             << setw(12) << map_mc[l]->avg_cost << endl;
-    }
-    mout << "--------------------------------------------------------------------------"
-         << endl;
+    MultilevelMonteCarlo mlmc(l_init, M_init, meshes,
+                              stochasticFields, assemble, graphs);
 
-    if (map_mc[L_mc]->kurtosis > 100.0) {
-        mout << endl << "WARNING: kurtosis on finest level = "
-             << map_mc[L_mc - 1]->kurtosis
-             << endl
-             << "indicates MultilevelMonteCarlo correction dominated by a few rare paths!"
-             << endl;
+    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);
     }
+}
 
-    double alpha = 0.0, beta = 0.0, gamma = 0.0;
-    estimateExponents(map_mc, alpha, beta, gamma, true);
-    mout << endl << "alpha = " << setw(8) << alpha
-         << "  (exponent for MultilevelMonteCarlo weak convergence)"
-         << endl << " beta = " << setw(8) << beta
-         << "  (exponent for MultilevelMonteCarlo variance)"
-         << endl << "gamma = " << setw(8) << gamma
-         << "  (exponent for MultilevelMonteCarlo cost)"
-         << endl;
+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")
 }
 
-void MLMCExperiment() {
-    double epsilon;
-    int Lmax;
-    vector<int> l_init, M_init;
-    ReadConfig(Settings, "epsilon", epsilon);
-    ReadConfig(Settings, "Lmax", Lmax);
-    ReadConfig(Settings, "l_init", l_init, 3);
-    ReadConfig(Settings, "M_init", M_init, 3);
+Assemble *getAssemble(const string &problemName, const string &modelName,
+                      Discretization &disc, StochasticProblem *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 == "DGFEM") Exit("Not implemented yet")
+        return nullptr;
+    } else if (problemName.find("Transport") != string::npos) {
+        if (modelName == "DGFEM") Exit("Not implemented yet")
+        return nullptr;
+    }
+}
 
-    mout << endl << "*********************************************************"
-         << endl << "***           MultilevelMonteCarlo test run           ***"
-         << endl << "*********************************************************"
-         << endl << endl;
+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);
+}
 
-    Meshes meshes(Lmax);
-    Discretization disc(1, meshes.dim());
-    StochasticFields stochFields(meshes);
-    StochasticProblem *problem = getStochasticProblem(stochFields,
-                                                      "StochLaplace2D");
-    EllipticAssemble ellipticAssemble(disc, *problem);
-    MatrixGraphs graphs(meshes, disc);
+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());
+    if (modelName == "HybridFEM")
+        return Discretization("RT0_P0", 1, meshes.dim());
+}
 
-    MultilevelMonteCarlo mlmc(l_init, M_init, Lmax, stochFields, ellipticAssemble,
-                              graphs);
-    mlmc.method(epsilon);
+void headConvergenceTest() {
+    mout << endl << "**********************************************************"
+         << endl << "*** Convergence tests, kurtosis, telescoping sum check ***"
+         << endl << "**********************************************************" << endl;
+}
 
-    mout << endl << "Value Q: " << mlmc.value
-         << endl << "Cost:    " << mlmc.cost
-         << endl << "Levels:  ";
-    for (auto &iter : mlmc.levels) mout << iter << " ";
-    mout << endl << "Samples: ";
-    for (auto &iter : mlmc.numsamples) mout << iter << " ";
-    mout << endl << endl;
+void headMLMCExperiment() {
+    mout << endl << "**********************************************************"
+         << endl << "***           MultilevelMonteCarlo test run            ***"
+         << endl << "**********************************************************" << endl;
 }
 
-//void MLMCExperimentOverEpsilon() {
-//    vector<double> epsilons;
-//    ReadConfig(Settings, "epsilon", epsilons, 3);
-//
-//    mout << endl << "*****************************"
-//         << endl << "*** MultilevelMonteCarlo complexity tests ***"
-//         << endl << "*****************************" << endl;
-//
-//    mout << endl << "  eps       value   MultilevelMonteCarlo cost             l M"
-//         << endl << "-----------------------------------------------------------------"
-//         << endl;
-//
-//    for (double &eps : epsilons) {
-//        MultilevelMonteCarlo mlmc = MultilevelMonteCarlo();
-//        mlmc.method(eps);
-//
-//        mout << setw(5) << eps << setw(12) << mlmc.value
-//             << setw(12) << mlmc.cost;
-//        string level_str, sample_str;
-//        for (auto &iter : mlmc.levels) level_str += to_string(iter) + " ";
-//        mout << setw(14) << level_str;
-//        for (auto &iter : mlmc.numsamples) sample_str += to_string(iter) + " ";
-//        mout << sample_str << endl;
-//    }
-//    mout << "-----------------------------------------------------------------" << endl;
-//
-//}
\ No newline at end of file
+void headMLMCOverEpsilon() {
+    mout << endl << "**********************************************************"
+         << endl << "***         MultilevelMonteCarlo over epsilon          ***"
+         << endl << "**********************************************************" << endl;
+}

mlmc/src/MLMCMain.h

@@ -2,13 +2,26 @@
 #define MLMC_MLMCMAIN_H
 
 #include "Elliptic.h"
+#include "Mixed.h"
+#include "Hybrid.h"
 #include "MonteCarlo.h"
 #include "MultilevelMonteCarlo.h"
 
-void ConvergenceTest();
+void MLMCMain();
 
-void MLMCExperiment();
+Meshes getMeshes(const string &problemName, int plevel, int lmax);
 
-void MLMCExperimentOverEpsilon();
+Discretization getDiscretization(const string &modelName, Meshes &meshes);
+
+Assemble *getAssemble(const string &problemName, const string &modelName,
+                      Discretization &disc, StochasticProblem *problem);
+
+MatrixGraphs getMatrixGraphs(const string &modelName, Meshes &meshes, Discretization &disc);
+
+void headConvergenceTest();
+
+void headMLMCExperiment();
+
+void headMLMCOverEpsilon();
 
 #endif //MLMC_MLMCMAIN_H

mlmc/src/Main.C

@@ -10,12 +10,6 @@ int main(int argv, char **argc) {
     SetSearchPath("../mlmc", argv, argc);
     DPO dpo(&argv, &argc);
 
-    string Experiment;
-    ReadConfig(Settings, "Experiment", Experiment);
-
-    if (Experiment == "ConvergenceTest") ConvergenceTest();
-    else if (Experiment == "MLMCExperiment") MLMCExperiment();
-//    else if (Experiment == "MLMCExperimentOverEpsilon") MLMCExperimentOverEpsilon();
-
+    MLMCMain();
     return 0;
 }