worked on transport

mlmc/conf/elliptic.conf

-Distribution = RCB
-
-LinearSolver = GMRES
-
-Preconditioner = LIB_PS
-#Preconditioner = SGS
-#Preconditioner = GaussSeidel
-#Preconditioner = Jacobi
-
-LinearReduction = 1e-12
-LinearEpsilon = 1e-10
-LinearSteps = 3000
-
-NewtonSteps = 1
-NewtonLineSearchSteps = 0
-
-TimeLevel = -1
-LinearVerbose = -1
-BaseSolverVerbose = -1
-MuteLevel = -1
-NewtonVerbose = -1
-DebugLevel = -1

mlmc/conf/m++.conf

-loadconf = elliptic.conf
-#loadconf = transport.conf
 loadconf = mlmc.conf
\ No newline at end of file

mlmc/conf/mlmc.conf

-Model = LagrangeFEM
+#Model = LagrangeFEM
 #Model = MixedFEM
 #Model = HybridFEM
 #Model = DGFEM
 Model = DGTransport
 
 Overlap = dG1
+flux_alpha = 1  #Upwind: 1, Central: 0
 #sign = 1
 #penalty = 20
 
@@ -20,7 +21,7 @@ StochasticField = LogNormal
 Experiment = MLMCExperiment
 #Experiment = MLMCOverEpsilon
 
-initLevels = 3, 4, 5
+initLevels = 5, 6, 7
 initSampleAmount = 12, 6, 3
 
 #l_init = 3,  4,  5
@@ -32,7 +33,7 @@ epsilon_lst = 0.05, 0.01, 0.005
 
 maxLevel = 8
 plevel = 2
-degree = 1
+degree = 2
 
 functional = L2
 #functional = Energy
@@ -40,7 +41,7 @@ functional = L2
 #functional = Outflow
 #functional = Goal
 
-MCVerbose = 1
+MCVerbose = 2
 MLMCVerbose = 2
 MCPlotting = 3
 MLMCPlotting = 2
@@ -54,4 +55,44 @@ norm_p = 2
 lambda1 = 0.1
 lambda2 = 0.1
 alpha = 1.8
-evtol = 1e-10
\ No newline at end of file
+evtol = 1e-10
+
+# Other
+Distribution = RCB
+LinearSolver = GMRES
+Preconditioner = LIB_PS
+#Preconditioner = SGS
+#Preconditioner = GaussSeidel
+#Preconditioner = Jacobi
+
+LinearReduction = 1e-12
+LinearEpsilon = 1e-10
+LinearSteps = 3000
+
+NewtonSteps = 1
+NewtonLineSearchSteps = 0
+
+TimeLevel = -1
+LinearVerbose = -1
+BaseSolverVerbose = -1
+MuteLevel = -1
+NewtonVerbose = -1
+DebugLevel = -1
+
+
+#Transport
+TimeSeries = uniform
+plot_tStep = 1
+T = 1.6
+dt = 0.005
+Kmax = 250
+Keps = 1e-5
+
+#rkorder = 0  #Exponential integrator
+#rkorder = 1  #Explicit Euler
+#rkorder = 2  #Heun
+#rkorder = 2  #Runge
+rkorder = -2  #Implicit MP-rule
+#rkorder = 1000  #Polynomial Krylov-Arnoldi method
+
+gamma = 0.01
\ No newline at end of file

mlmc/src/MLMCMain.C

@@ -59,9 +59,9 @@ MatrixGraphs *MLMCMain::getSampleMatrixGraphs() {
 }
 
 StochasticField *MLMCMain::getStochasticField() {
-    if (problemName == "StochasticLaplace1D" || "StochasticLaplace2D")
+    if (problemName == "StochasticLaplace1D" ||  problemName == "StochasticLaplace2D")
         return new StochasticField(*meshes, "CirculantEmbedding");
-    if (problemName == "StochasticPollution1D" || "StochasticPollution2D")
+    if (problemName == "StochasticPollution1D" || problemName == "StochasticPollution2D")
         return new StochasticField(*meshes, "HybridFluxGenerator");
 }
 

mlmc/src/mc/MonteCarlo.h

@@ -77,6 +77,16 @@ protected:
         kurtosis = (avgY4 - 4.0 * avgY3 * avgY + 6.0 * avgY2 * avgY * avgY -
             3.0 * avgY * avgY * avgY * avgY) / varY / varY;
     }
+
+    const char *buildString(int m, string name, const string &suffix="") {
+        if (!suffix.empty())
+            name += "_" + suffix;
+        if (l != -1)
+            name += "_" + to_string(l);
+        if (m != -1)
+            name += "_" + to_string(m);
+        return const_cast<char *>(name.c_str());
+    }
 };
 
 #endif //MLMC_MC_HPP
\ No newline at end of file

mlmc/src/mc/MonteCarloElliptic.C

@@ -88,26 +88,12 @@ void MonteCarloElliptic::plotResults(int l, int m, const string &suffix,
                                      const Vector &u, const Vector &flux,
                                      const Vector &perm) {
 
-    const char *name = buildString(l, m, "u", suffix);
+    const char *name = buildString(m, "u", suffix);
     assemble->plotU(plot, u, name);
 
-    name = buildString(l, m, "flux", suffix);
+    name = buildString(m, "flux", suffix);
     assemble->plotFlux(plot, flux, name);
 
-    name = buildString(l, m, "perm", suffix);
+    name = buildString(m, "perm", suffix);
     assemble->plotPerm(plot, perm, name);
 }
-
-const char *MonteCarloElliptic::buildString(int l,
-                                            int m,
-                                            string name,
-                                            const string &suffix) {
-    if (!suffix.empty())
-        name += "_" + suffix;
-    if (l != -1)
-        name += "_" + to_string(l);
-    if (m != -1)
-        name += "_" + to_string(m);
-    return const_cast<char *>(name.c_str());
-}
-

mlmc/src/mc/MonteCarloElliptic.h

@@ -12,8 +12,6 @@ private:
     void solveElliptic(int l, double &valueQ, double &cost,
                        Vector &solution);
 
-    static const char *buildString(int l, int m, string name, const string &suffix);
-
     virtual void plotResults(int l, int m, const string &suffix,
                              const Vector &u, const Vector &flux,
                              const Vector &perm);

mlmc/src/mc/MonteCarloTransport.C

 #include "MonteCarloTransport.h"
 
+
 using namespace std;
 
 void MonteCarloTransport::Method() {
@@ -12,39 +13,61 @@ void MonteCarloTransport::Method() {
                              + " dM: " + to_string(dM)
                              + " onlyFineLevel: True");
 
-    MatrixGraphs fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3)));
-    Vector fineSolution(solMatrixGraphs[l - pLevel]);
     Vector fineFlux(fluxMatrixGraphs[l - pLevel]);
-    Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
     Vector coarseFlux(fluxMatrixGraphs[l - pLevel - 1]);
+    Vector coarseFluxUp(fineFlux);
+    Vector fineSolution(solMatrixGraphs[l - pLevel]);
+    Vector coarseSolution(solMatrixGraphs[l - pLevel - 1]);
+    Vector coarseSolutionUp(fineSolution);
 
-    fineFlux = 0.0;
+    fineFlux = 0.0, coarseFlux = 0.0, coarseFluxUp = 0.0;
+    fineSolution = 0.0, coarseSolution = 0.0, coarseSolutionUp = 0.0;
 
     double fineQ, coarseQ, fineCost, coarseCost;
     fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
 
+    assemble->plot = &plot;
+
     for (int m = M; m < M + dM; m++) {
         stochasticField->GetFineSample(l, fineFlux);
         assemble->problem->LoadNewSample(&fineFlux);
 
+        if (plotting > 1) {
+            const char *name = buildString(m, "flux");
+            plot.celldata(fineFlux, 3);
+            plot.vtk_cellvector(name);
+        }
+
         solveTransport(l, fineQ, fineCost, fineSolution);
 
         if (!onlyFineLevel) {
             stochasticField->GetCoarseSample(l - 1, fineFlux, coarseFlux);
             assemble->problem->LoadNewSample(&coarseFlux);
             solveTransport(l, fineQ, fineCost, fineSolution);
+        } else {
+            coarseQ = 0.0, coarseCost = 0.0, coarseSolutionUp = 0.0;
         }
+        updateSums(fineCost, fineQ - coarseQ, fineQ);
     }
+    updateAvg();
+    updateStatistic();
+
+    logger->logMSGV1("|E[Y]|: " + to_string(avgY)
+                         + " V[Y]: " + to_string(varY));
+    logger->endMethodMSG();
 }
 
-void MonteCarloTransport::solveTransport(int l,
-                                         double &valueQ,
-                                         double &cost,
-                                         Vector &u) {
+void MonteCarloTransport::solveTransport(int l, double &valueQ, double &cost,
+                                         Vector &solution) {
+    Date startSolving;
+
     Solver solver;
     Preconditioner *BJ = GetPC("PointBlockJacobi");
     Solver IM(BJ);
     DGTimeIntegrator timeIntegrator(solver, solver, solver, IM);
-//    assemble.PrintInfo(u);
-    timeIntegrator(*assemble, timeSeries, u);
+//    assemble.PrintInfo(solution);
+    timeIntegrator(*assemble, timeSeries, solution);
+
+    logger->logMSGV2("Solving on level " + to_string(l) + " took " +
+        to_string((startSolving - Date()).t) + " seconds");
 }
\ No newline at end of file

mlmc/src/mc/MonteCarloTransport.h

@@ -8,21 +8,27 @@
 
 class MonteCarloTransport : public MonteCarlo {
 private:
+    Plot plot;
+
     void solveTransport(int l, double &valueQ, double &cost, Vector &solution);
 
 public:
     DGTransportAssemble *assemble;
     TimeSeries timeSeries;
 
+    MatrixGraphs fluxMatrixGraphs;
     CellMatrixGraphs solMatrixGraphs;
 
+
     Transfer *transfer;
 
     MonteCarloTransport(int l, int dM, bool baseLevel, Meshes *meshes,
                         StochasticField *stochFields, DGTransportAssemble *assemble) :
         MonteCarlo(l, dM, baseLevel, meshes, stochFields),
         assemble(assemble),
-        solMatrixGraphs(*meshes, *assemble->disc) {
+        fluxMatrixGraphs(MatrixGraphs(*meshes, dof("cell", 3))),
+        solMatrixGraphs(*meshes, *assemble->disc),
+        plot(Plot((*meshes)[l], meshes->dim())) {
 
 //        transfer = new MatrixTransfer(*assemble);
 //        transfer->Construct((solMatrixGraphs)[l - pLevel], (solMatrixGraphs)[l - pLevel - 1]);