added error functions

mlmc/src/MultilevelMonteCarlo.C

@@ -28,6 +28,7 @@ void MultilevelMonteCarlo::initializeValues() {
     levels = {}, numsamples = {};
     value = 0.0, cost = 0.0;
     alpha = 0.0, beta = 0.0, gamma = 0.0;
+    numErr = 0.0, statErr = 0.0, totalErr = 0.0;
 }
 
 void MultilevelMonteCarlo::clearMapMonteCarlo() {
@@ -55,10 +56,14 @@ void MultilevelMonteCarlo::Method(const double eps) {
         estimateExponents();
         updateSampleAmount(eps);
         if (checkForNewLevel()) {
-            if (estimateNumericalError() > eps / sqrt(2.0))
+            estimateNumericalError();
+            if (numErr > eps) {
                 appendLevel(eps);
-            else
+            } else {
+                estimateStatisticalError();
+                totalErr = statErr + numErr;
                 converged = true;
+            }
         }
     }
     wrapUpResults(value, cost, levels, numsamples);
@@ -132,19 +137,25 @@ bool MultilevelMonteCarlo::checkForNewLevel() {
     return newLevel;
 }
 
-double MultilevelMonteCarlo::estimateNumericalError() {
-    double maxErr = 0.0;
+void MultilevelMonteCarlo::estimateNumericalError() {
+    numErr = 0.0;
     unsigned long exponent = mapMonteCarlo.size() - 2;
     for (auto &mc : mapMonteCarlo) {
         if (mc.first == mapMonteCarlo.begin()->first) continue;
-        maxErr = max(mc.second->avgY / (pow(2.0, alpha) - 1) /
-                         (pow(2.0, exponent)),
-                     maxErr);
+        numErr = max(mc.second->avgY / (pow(2.0, alpha) - 1) / (pow(2.0, exponent)),
+                     numErr);
         exponent--;
     }
 
-    logger->InnerMsg("err: " + to_string(maxErr), verbose);
-    return maxErr;
+    logger->InnerMsg("estimated numerical error: " + to_string(numErr), verbose);
+}
+
+void MultilevelMonteCarlo::estimateStatisticalError() {
+    statErr = 0.0;
+    for (auto &mc : mapMonteCarlo) {
+        statErr += mc.second->varY / mc.second->M;
+    }
+    statErr = sqrt(statErr);
 }
 
 void MultilevelMonteCarlo::appendLevel(double eps) {
@@ -231,12 +242,18 @@ void MultilevelMonteCarlo::ShowKurtosisWarning() {
 
 void MultilevelMonteCarlo::ShowExponentResults() {
     estimateExponents();
-    mout << endl << "alpha = " << setw(8) << alpha
-         << "  (exponent for the weak convergence of the FEM)"
-         << endl << " beta = " << setw(8) << beta
-         << "  (exponent for the decay of the variance)"
-         << endl << "gamma = " << setw(8) << gamma
-         << "  (exponent for the growth of the cost)"
+    mout << 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)"
+         << endl << "gamma   = " << setw(10) << gamma
+         << "   (exponent for the growth of the cost)"
+         << endl << "numErr  = " << setw(10) << numErr
+         << "   (estimated numerical error)"
+         << endl << "statErr = " << setw(10) << statErr
+         << "   (estimated statistical error)"
+         << endl << "error   = " << setw(10) << totalErr
+         << "   (estimated total error)"
          << endl;
 }
 

mlmc/src/MultilevelMonteCarlo.h

@@ -27,7 +27,9 @@ private:
 
     void updateSampleAmount(const double &eps);
 
-    double estimateNumericalError();
+    void estimateNumericalError();
+
+    void estimateStatisticalError();
 
     void appendLevel(double eps);
 
@@ -55,6 +57,7 @@ public:
     vector<double> valueOverEpsilon = {}, costOverEpsilon = {};
 
     double alpha = 0.0, beta = 0.0, gamma = 0.0;
+    double numErr = 0.0, statErr = 0.0, totalErr = 0.0;
 
     MultilevelMonteCarlo(vector<int> &initLevels,
                          vector<int> &initSampleAmount,

mlmc/src/mc/MonteCarloElliptic.C

@@ -22,8 +22,6 @@ void MonteCarloElliptic::method() {
         assemble->plot = finePlot;
         assemble->problem->LoadNewSample(&finePermeability);
 
-        mout << finePermeability << endl;
-
         solvePDE(l, m, true, fineQ, fineCost, fineSolution);
 
         if (plotting) assemble->PlotResults(fineSolution);

mlmc/src/mc/MonteCarloTransport.C

@@ -71,7 +71,7 @@ void MonteCarloTransport::solvePDE(int l,
     if (functional == "Energy") valueQ = assemble->Energy(solution);
     if (functional == "Outflow")
         valueQ = assemble->InFlowOutFlowRate(solution).second;
-    mout << endl;
+
     logger->EndMethod(verbose);
 }