Major refactoring in MLMC Framework

CMakeLists.txt

@@ -75,8 +75,8 @@ target_link_libraries(MLMC-M++ MLMC sprng SRC LIB_PS ${SUPERLU} blas lapack fftw
 #add_executable(TestRNManager mlmc/tests/TestRNManager.C)
 
 # Linking
-#target_link_libraries(TestCirculantEmbedding MLMC sprng SRC SRC_Solver
+#target_link_libraries(TestCirculantEmbedding MultilevelMonteCarlo sprng SRC SRC_Solver
 #       blas lapack ${SUPERLU} fftw3 m ${GTEST_LIB})
-#target_link_libraries(TestRNManager MLMC sprng fftw3 m ${GTEST_LIB})
+#target_link_libraries(TestRNManager MultilevelMonteCarlo sprng fftw3 m ${GTEST_LIB})
 #set(MPI_COMMAND mpirun -np 4 )
 #---------------------------------------------------------------------------------------#

mlmc/conf/mlmc.conf

 #----- ConvergenceTest -----#
-Experiment = ConvergenceTest
-SubroutineEstimator = MCElliptic
-l_mc = 3
-L_mc = 7
-M = 50
+#Experiment = ConvergenceTest
+#l_mc = 3
+#L_mc = 7
+#M = 10
 
 #----- MLMCExperiment -----#
-#Experiment = MLMCExperiment
-#SubroutineEstimator = MCElliptic
-#epsilon = 0.5
-#l_init = 3 4 5
-#M_init = 40 20 10
-#Lmax = 9
+Experiment = MLMCExperiment
+epsilon = 0.01
+l_init = 3 4 5
+M_init = 40 20 10
+Lmax = 8
 
 #----- MLMCExperimentOverEpsilon-----#
 #Experiment = MLMCExperimentOverEpsilon
-#SubroutineEstimator = MCElliptic
+#SubroutineEstimator = MonteCarloElliptic
 #epsilon = 0.5 0.3 0.1
 #l_init = 3 4 5
 #M_init = 40 20 5
@@ -31,4 +29,8 @@ functional = L2
 MCVerbose = 0
 MLMCVerbose = 0
 MCPlotting = 0
+MLMCPlotting = 0
+
+vtkplot = 1
+
 

mlmc/src/CMakeLists.txt

 set(MLMC_SRC
         Main.C
-        MLMC.C
-        MC.C
         Utils.C
+        MLMCMain.C
         RNManager.C
-        ConvergenceTest.C
-        MLMCExperiment.C
-        MCElliptic.C
+        MonteCarlo.C
+        StochasticFields.C
+        StochasticProblem.C
         CirculantEmbedding.C
-        PlottingUtils.C
+        MultilevelMonteCarlo.C
         )
 
 add_library(MLMC STATIC ${MLMC_SRC})

mlmc/src/CirculantEmbedding.h

@@ -4,9 +4,6 @@
 #include "RNManager.h"
 #include "CovarianceFunction.h"
 #include "Utils.h"
-#include <math.h>
-#include <complex>
-#include <vector>
 
 using namespace std;
 

mlmc/src/ConvergenceTest.C

-#include "m++.h"
-#include "Utils.h"
-#include "MCProblem.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;
-
-    mout << endl << " l    E(Qf-Qc)       E(Qf)    V(Qf-Qc)       V(Qf)    kurtosis        cost"
-         << endl << "--------------------------------------------------------------------------" << endl;
-
-    map<int, MC *> map_mc;
-    for (int l = l_mc; l <= L_mc; l++) {
-        bool base_l = l == l_mc;
-        map_mc[l] = GetSubroutineEstimator(l, M, base_l);
-        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;
-
-    if (map_mc[L_mc]->kurtosis > 100.0) {
-        mout << endl << "WARNING: kurtosis on finest level = " << map_mc[L_mc - 1]->kurtosis
-             << endl << "indicates MLMC correction dominated by a few rare paths!" << endl;
-    }
-
-    double alpha = 0.0, beta = 0.0, gamma = 0.0;
-    estimate_exponents(map_mc, alpha, beta, gamma, true);
-    mout << endl << "alpha = " << setw(8) << alpha << "  (exponent for MLMC weak convergence)"
-         << endl << " beta = " << setw(8) << beta << "  (exponent for MLMC variance)"
-         << endl << "gamma = " << setw(8) << gamma << "  (exponent for MLMC cost)" << endl;
-}

mlmc/src/Elliptic.h

+#ifndef _LAPLACE_H_
+#define _LAPLACE_H_
+
+#include "m++.h"
+#include "StochasticProblem.h"
+
+class EllipticAssemble : public Assemble {
+public:
+    const Discretization &disc;
+    const StochasticProblem &problem;
+
+    EllipticAssemble(const Discretization &disc, const StochasticProblem &problem)
+            : disc(disc), problem(problem) {}
+
+    void PrintInfo(const Vector &u) override {
+        Assemble::PrintInfo(u);
+        Vout(1) << " Problem:        " << problem.Name() << endl
+                << " Discretization: " << disc.DiscName() << endl
+                << endl;
+    }
+
+    virtual Element *getElement(const cell &c, const Vector &u) const {
+        return new ScalarElement(disc, u, c);
+    }
+
+    virtual FaceElement *getFaceElement(const cell &c, const Vector &u, int face) const {
+        return new ScalarFaceElement(disc, u, c, face);
+    }
+
+    virtual Scalar getValue(const cell &c, const Vector &u, Element *elem, int q) const {
+        auto scalarElem = dynamic_cast<ScalarElement *>(elem);
+        return scalarElem->Value(q, u);
+    }
+
+    virtual Scalar getFaceValue(const cell &c, const Vector &u,
+                                FaceElement *faceElem, int q, int face) const {
+        auto scalarFaceElem = dynamic_cast<ScalarFaceElement *>(faceElem);
+        return scalarFaceElem->Value(q, u);
+    }
+
+    virtual VectorField getFlux(const cell &c, const Vector &u,
+                                Element *elem, int q) const {
+        auto scalarElem = dynamic_cast<ScalarElement *>(elem);
+        Tensor K = problem.Permeability(c.Subdomain(), elem->QPoint(q));
+        return K * scalarElem->Derivative(q, u);
+    }
+
+    virtual VectorField getFaceFlux(const cell &c, const Vector &u,
+                                    FaceElement *faceElem, int q, int face) const {
+        auto scalarFaceElem = dynamic_cast<ScalarFaceElement *>(faceElem);
+        Tensor K = problem.Permeability(c.Subdomain(), scalarFaceElem->QPoint(q));
+        return K * scalarFaceElem->Derivative(q, u);
+    }
+
+    virtual VectorField getGradient(const cell &c, const Vector &u,
+                                    Element *elem, int q) const {
+        auto scalarElem = dynamic_cast<ScalarElement *>(elem);
+        return scalarElem->Derivative(q, u);
+    }
+
+    const char *Name() const override { return "Lagrange"; }
+
+    void Dirichlet(const cell &c, Vector &u) const override {
+        RowBndValues u_c(u, c);
+        if (!u_c.onBnd()) return;
+        Element *elem = getElement(c, u);
+        for (int face = 0; face < c.Faces(); ++face) {
+            if (u_c.bc(face) == 1) {
+                for (int j = 0; j < disc.NodalPointsOnFace(c, face); ++j) {
+                    int k = disc.NodalPointOnFace(c, face, j);
+                    u_c(k) = problem.Solution(elem->NodalPoint(k));
+                    u_c.D(k) = true;
+                }
+            }
+        }
+        delete elem;
+    }
+
+    void Residual(const cell &c, const Vector &u, Vector &r) const override {
+        ScalarElement elem(disc, u, c);
+        RowBndValues r_c(r, c);
+        for (int q = 0; q < elem.nQ(); ++q) {
+            double w = elem.QWeight(q);
+            const Point &z = elem.QPoint(q);
+            Tensor K = problem.Permeability(c.Subdomain(), z);
+            Scalar f = problem.Load(c.Subdomain(), z);
+            VectorField DU = elem.Derivative(q, u);
+            VectorField K_DU = K * DU;
+            for (int i = 0; i < elem.size(); ++i) {
+                Scalar U_i = elem.Value(q, i);
+                VectorField DU_i = elem.Derivative(q, i);
+                r_c(i) += w * (K_DU * DU_i - f * U_i);
+            }
+        }
+        if (!r_c.onBnd()) return;
+        for (int f = 0; f < c.Faces(); ++f) {
+            int bnd = r_c.bc(f);
+            if (bnd != 2) continue;
+            ScalarFaceElement felem(disc, u, c, f);
+            RowValues r_f(r, felem);
+            for (int q = 0; q < felem.nQ(); ++q) {
+                double w = felem.QWeight(q);
+                Scalar F = problem.Flux(bnd, felem.QPoint(q)) * felem.QNormal(q);
+                for (int i = 0; i < felem.size(); ++i) {
+                    Scalar U_i = felem.Value(q, i);
+                    r_f(i) -= w * U_i * F;
+                }
+            }
+        }
+    }
+
+    void Jacobi(const cell &c, const Vector &u, Matrix &A) const override {
+        ScalarElement E(disc, u, c);
+        RowEntries A_c(A, E);
+        for (int q = 0; q < E.nQ(); ++q) {
+            double w = E.QWeight(q);
+            Tensor K = problem.Permeability(c.Subdomain(), E.QPoint(q));
+            for (unsigned int i = 0; i < E.size(); ++i) {
+                VectorField Dv = K * E.Derivative(q, i);
+                for (unsigned int j = 0; j < E.size(); ++j) {
+                    VectorField Dw = E.Derivative(q, j);
+                    A_c(i, j) += w * Dv * Dw;
+                }
+            }
+        }
+    }
+
+    double Energy(const Vector &u) const override {
+        double energy = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            for (int q = 0; q < elem->nQ(); ++q) {
+                double w = elem->QWeight(q);
+                VectorField DU = getGradient(c, u, elem, q);
+                VectorField K_DU = getFlux(c, u, elem, q);
+                energy += w * K_DU * DU;
+            }
+            delete elem;
+        }
+        return sqrt(PPM->Sum(energy));
+    }
+
+    virtual double EnergyError(const Vector &u) const {
+        double err = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            for (int q = 0; q < elem->nQ(); ++q) {
+                double w = elem->QWeight(q);
+                Tensor IK = Invert(problem.Permeability(c.Subdomain(), elem->QPoint(q)));
+                VectorField diff = (getFlux(c, u, elem, q) -
+                                    problem.Flux(c.Subdomain(), elem->QPoint(q)));
+                err += (IK * diff) * diff;
+            }
+            delete elem;
+        }
+        return sqrt(PPM->Sum(err));
+    }
+
+    virtual double L2(const Vector &u) const {
+        double l2 = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            for (int q = 0; q < elem->nQ(); ++q) {
+                double w = elem->QWeight(q);
+                Scalar U = getValue(c, u, elem, q);
+                l2 += w * U * U;
+            }
+            delete elem;
+        }
+        return sqrt(PPM->Sum(l2));
+    }
+
+    virtual double L2Error(const Vector &u) const {
+        double err = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            for (int q = 0; q < elem->nQ(); ++q) {
+                double w = elem->QWeight(q);
+                Scalar U = getValue(c, u, elem, q);
+                Scalar Sol = problem.Solution(elem->QPoint(q));
+                err += w * (U - Sol) * (U - Sol);
+            }
+            delete elem;
+        }
+        return sqrt(PPM->Sum(err));
+    }
+
+    virtual double L2CellAverageError(const Vector &u) const {
+        double err = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            double w = 0;
+            Scalar U = 0;
+            Scalar Sol = 0;
+            for (int q = 0; q < elem->nQ(); ++q) {
+                w += elem->QWeight(q);
+                U += getValue(c, u, elem, q);
+                Sol += problem.Solution(elem->QPoint(q));
+            }
+            err += w * (U - Sol) * (U - Sol);
+            delete elem;
+        }
+        return sqrt(PPM->Sum(err));
+    }
+
+    virtual double MaxError(const Vector &u) const {
+        double err = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            for (int q = 0; q < elem->nQ(); ++q) {
+                Scalar U = getValue(c, u, elem, q);
+                err = max(err, abs(U - problem.Solution(elem->QPoint(q))));
+            }
+            delete elem;
+        }
+        return PPM->Max(err);
+    }
+
+    virtual double FluxError(const Vector &u) const {
+        double flux_error = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            BFParts bnd(u.GetMesh(), c);
+            if (!bnd.onBnd()) continue;
+            for (int face = 0; face < c.Faces(); ++face) {
+                if (bnd[face] == 2) {
+                    auto *faceElem = getFaceElement(c, u, face);
+                    for (int q = 0; q < faceElem->nQ(); ++q) {
+                        double w = faceElem->QWeight(q);
+                        VectorField G = problem.Flux(bnd[face], faceElem->QPoint(q));
+                        Scalar F = (getFaceFlux(c, u, faceElem, q, face) - G) *
+                                   faceElem->QNormal(q);
+                        flux_error += w * F * F;
+                    }
+                    delete faceElem;
+                } else if (bnd[face] == 0) {
+                    auto *faceElem = getFaceElement(c, u, face);
+                    for (int q = 0; q < faceElem->nQ(); ++q) {
+                        double w = faceElem->QWeight(q);
+                        Scalar F = getFaceFlux(c, u, faceElem, q, face)
+                                   * faceElem->QNormal(q);
+                        flux_error += w * F * F;
+                    }
+                    delete faceElem;
+                }
+            }
+        }
+        return sqrt(PPM->Sum(flux_error));
+    }
+
+    virtual double FaceError(const Vector &u) const {
+        double face_error = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            for (int face = 0; face < c.Faces(); ++face) {
+                auto *faceElem = getFaceElement(c, u, face);
+                for (int q = 0; q < faceElem->nQ(); ++q) {
+                    double w = faceElem->QWeight(q);
+                    Scalar U = getFaceValue(c, u, faceElem, q, face);
+                    Scalar Sol = problem.Solution(faceElem->QPoint(q));
+                    face_error += w * (U - Sol) * (U - Sol);
+                }
+                delete faceElem;
+            }
+        }
+        return sqrt(PPM->Sum(face_error));
+    }
+
+    virtual pair<double, double> getInflowOutflow(const Vector &u) const {
+        double inflow = 0;
+        double outflow = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            BFParts bnd(u.GetMesh(), c);
+            if (!bnd.onBnd()) continue;
+            for (int face = 0; face < c.Faces(); ++face) {
+                if (bnd[face] < 0) continue;
+                auto *faceElem = getFaceElement(c, u, face);
+                for (int q = 0; q < faceElem->nQ(); ++q) {
+                    double w = faceElem->QWeight(q);
+                    Scalar F = getFaceFlux(c, u, faceElem, q, face) *
+                               faceElem->QNormal(q);
+                    if (F > 0) outflow += w * F;
+                    else inflow += w * F;
+                }
+                delete faceElem;
+            }
+        }
+        return {PPM->Sum(inflow), PPM->Sum(outflow)};
+    }
+
+    virtual pair<double, double> getPrescribedInflowOutflow(const Vector &u) const {
+        double inflow = 0;
+        double outflow = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            BFParts bnd(u.GetMesh(), c);
+            if (!bnd.onBnd()) continue;
+            for (int face = 0; face < c.Faces(); ++face) {
+                if (bnd[face] < 2) continue;
+                auto *faceElem = getFaceElement(c, u, face);
+                for (int q = 0; q < faceElem->nQ(); ++q) {
+                    double w = faceElem->QWeight(q);
+                    Scalar F = problem.Flux(bnd[face], faceElem->QPoint(q)) *
+                               faceElem->QNormal(q);
+                    if (F > 0) outflow += w * F;
+                    else inflow += w * F;
+                }
+                delete faceElem;
+            }
+        }
+        return {PPM->Sum(inflow), PPM->Sum(outflow)};
+    }
+
+    virtual pair<double, double> getOutflowLeftRight(const Vector &u) const {
+        double outflowLeft = 0;
+        double outflowRight = 0;
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            BFParts bnd(u.GetMesh(), c);
+            if (!bnd.onBnd()) continue;
+            for (int face = 0; face < c.Faces(); ++face) {
+                if (bnd[face] != 1) continue;
+                auto *faceElem = getFaceElement(c, u, face);
+                for (int q = 0; q < faceElem->nQ(); ++q) {
+                    double w = faceElem->QWeight(q);
+                    Scalar F = getFaceFlux(c, u, faceElem, q, face) *
+                               faceElem->QNormal(q);
+                    if (F > 0 && c()[0] <= 1) outflowLeft += w * F;
+                    else if (F > 0 && c()[0] >= 2) outflowRight += w * F;
+                }
+                delete faceElem;
+            }
+        }
+        return {PPM->Sum(outflowLeft), PPM->Sum(outflowRight)};
+    }
+
+    virtual void setExactSolution(Vector &u) const {
+        for (row r = u.rows(); r != u.rows_end(); ++r)
+            u(r, 0) = problem.Solution(r());
+    }
+
+    virtual void setFlux(const Vector &u, Vector &flux) {
+        for (cell c = u.cells(); c != u.cells_end(); ++c) {
+            auto *elem = getElement(c, u);
+            VectorField F = zero;
+            double area = 0;
+            for (int q = 0; q < elem->nQ(); ++q) {
+                double w = elem->QWeight(q);
+                VectorField Q = getFlux(c, u, elem, q);
+                F += w * Q;
+                area += w;
+            }
+            delete elem;
+            F *= (1 / area);
+            for (int d = 0; d < 3; ++d)
+                flux(c(), d) = F[d];
+        }
+    }
+
+    virtual void plotU(Plot &plot, const Vector &u, Meshes &meshes) {
+        plot.vertexdata(u);
+        plot.vtk_vertexdata("u");
+    }
+
+    virtual void plotUex(Plot &plot, const Vector &u, Meshes &meshes) {
+        if (!problem.SetSolution()) return;
+        Vector u_ex(u);
+        setExactSolution(u_ex);
+        plot.vertexdata(u_ex);
+        plot.vtk_vertexdata("u_ex");
+    }
+
+    virtual void plotFlux(Plot &plot, const Vector &u, Meshes &meshes) {
+        dof fluxdof("cell", 3);
+        MatrixGraphs fluxMGraph(meshes, fluxdof);
+        Vector flux(fluxMGraph.fine());
+        setFlux(u, flux);
+        plot.celldata(flux, 3);
+        plot.vtk_cellvector("flux");
+        for (row r = flux.rows(); r != flux.rows_end(); ++r) {
+            VectorField F(flux(r, 0), flux(r, 1), flux(r, 2));
+            flux(r, 0) = norm(F);
+        }
+        plot.celldata(flux, 1);
+        plot.vtk_celldata("norm_flux");
+    }
+
+    virtual void plotPermeability(Plot &plot, const Vector &u, Meshes &meshes) {
+        dof permdof("cell", 3);
+        MatrixGraphs permMGraph(meshes, permdof);
+        Vector permeability(permMGraph.fine());
+        for (cell c = permeability.cells(); c != permeability.cells_end(); ++c) {
+            Tensor K = problem.Permeability(c.Subdomain(), c());
+            permeability(c(), 0) = K[0][0];
+        }
+        plot.celldata(permeability);
+        plot.vtk_celldata("permeability");
+    }
+
+    void AssembleTransfer(TransferMatrix &TM) const override {
+        TM = 0;
+        const matrixgraph &cg = TM.CoarseMatrixGraph();
+        for (cell c = cg.cells(); c != cg.cells_end(); ++c) {
+            for (int cf = 0; cf < c.Corners(); ++cf)
+                TM(c.Corner(cf), c.Corner(cf))[0] = 1;
+            for (int ce = 0; ce < c.Edges(); ++ce)
+                for (int cf = 0; cf < 2; ++cf)
+                    TM(c.EdgeCorner(ce, cf), c.Edge(ce))[0] = 0.5;
+            if (c.Corners() == 4)
+                for (int cf = 0; cf < c.Corners(); ++cf)
+                    TM(c.Corner(cf), c())[0] = 0.25;
+        }
+    }
+
+    void printSolutionInfo(const Vector &u) {
+        pair<double, double> prescribedFlows = getPrescribedInflowOutflow(u);
+        pair<double, double> calculatedFlows = getInflowOutflow(u);
+
+        mout << endl << endl;
+        mout << "Prescribed Inflow:   " << setw(3) << prescribedFlows.first
+             << "  Calculated Inflow:   " << setw(3) << calculatedFlows.first << endl
+             << "Prescribed Outflow:  " << setw(3) << prescribedFlows.second
+             << "  Calculated Outflow:  " << setw(3) << calculatedFlows.second
+             << endl << endl;
+
+        mout << "Flux Loss           "
+             << calculatedFlows.first + calculatedFlows.second << endl;
+        mout << "Flux Error          " << FluxError(u)
+             << endl << endl;
+
+        if (problem.Name() == "Rock Problem") {
+            pair<double, double> outflowLeftRight = getOutflowLeftRight(u);
+
+            mout << "LeftOutflow:        " << outflowLeftRight.first << endl
+                 << "RightOutflow:       " << outflowLeftRight.second
+                 << endl << endl;
+        }
+
+        if (problem.SetSolution()) {
+            mout << "Exact solution available...       " << endl;
+            mout << "Supremum Error                    " << MaxError(u) << endl;
+            mout << "Energy Error                      " << EnergyError(u) << endl;
+            mout << "L2 Error                          " << L2Error(u) << endl;
+            mout << "L2 Error Average                  " << L2CellAverageError(u);
+            mout << endl << endl;
+        }
+    }
+
+    void plotSolution(const Vector &u, int l) {
+        Meshes meshes(l);
+        Plot plot(meshes[l], meshes.dim(), 2);
+        if (plot.vtk_plot() == 1) {
+            plotU(plot, u, meshes);
+            plotUex(plot, u, meshes);
+            plotFlux(plot, u, meshes);
+            plotPermeability(plot, u, meshes);
+            plot.vtk_load();
+        }
+    }
+};
+
+#endif

mlmc/src/MC.C

-#include "MC.h"
-
-using namespace std;