Merge branch '21-transport-tests-and-notebook' into 'feature'

added transport test to testpdesolver

Closes #21

See merge request !44

mluq/src/problems/StochasticTransportProblem.cpp

@@ -3,30 +3,33 @@
 
 IStochasticTransportProblem *
 CreateStochasticTransportProblem(const string &problemName, const Meshes &meshes) {
-    if (problemName == "StochasticPollution1D")
-        return new StochasticPollution1D(meshes);
-    if (problemName == "Pollution1D")
-        return new Pollution1D(meshes);
+  if (problemName == "StochasticPollution1D")
+    return new StochasticPollution1D(meshes);
+  if (problemName == "Pollution1D")
+    return new Pollution1D(meshes);
 
-    if (problemName == "StochasticPollutionCosHat1D")
-        return new StochasticPollutionCosHat1D(meshes);
+  if (problemName == "StochasticPollutionCosHat1D")
+    return new StochasticPollutionCosHat1D(meshes);
 //    if (problemName == "PollutionCosHat1D")
 //        return new PollutionCosHat1D(meshes);
 
-    if (problemName == "StochasticPollution2D")
-        return new StochasticPollution2D(meshes);
-    if (problemName == "Pollution2D")
-        return new Pollution2D(meshes);
+  if (problemName == "StochasticPollution2D")
+    return new StochasticPollution2D(meshes);
+  if (problemName == "Pollution2D")
+    return new Pollution2D(meshes);
 
-    if (problemName == "StochasticPollutionCosHat2D")
-        return new StochasticPollutionCosHat2D(meshes);
-    if (problemName == "PollutionCosHat2D")
-        return new PollutionCosHat2D(meshes);
+  if (problemName == "StochasticPollutionCosHat2D")
+    return new StochasticPollutionCosHat2D(meshes);
+  if (problemName == "PollutionCosHat2D")
+    return new PollutionCosHat2D(meshes);
 
-    if (problemName == "StochasticGaussHat2D")
-        return new StochasticGaussHat2D(meshes);
-    if (problemName == "GaussHat2D")
-        return new GaussHat2D(meshes);
+  if (problemName == "StochasticGaussHat2D")
+    return new StochasticGaussHat2D(meshes);
+  if (problemName == "GaussHat2D")
+    return new GaussHat2D(meshes);
 
-    Exit(problemName + " not found")
+  if (problemName == "CosHat2D")
+    return new CosHat2D(meshes);
+
+  Exit(problemName + " not found")
 }

mluq/src/problems/StochasticTransportProblem.hpp

@@ -23,7 +23,7 @@ public:
     config.get("T", T);
   }
 
-  double GetStepSize(int h) const { return CFL * h; }
+  double GetStepSize(double h) const { return CFL * h; }
 
   double GetStartTime() const { return t0; }
 
@@ -276,6 +276,39 @@ public:
   string Name() const override { return "PollutionCosHat2D"; }
 };
 
+class CosHat2D : public IStochasticTransportProblem {
+  double amplitude = 1.0;
+  double frequency = 8.0;
+public:
+  explicit CosHat2D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  double Solution(double t, const Point &x) const override {
+    Point midPoint = 0.25 * (Point(cos(2.0 * Pi * t), sin(2.0 * Pi * t)))
+      + Point(0.5, 0.5);
+    double r = dist(midPoint, x);
+    if (r < (Pi / 2.0) / frequency)
+      return amplitude * pow(cos(frequency * r), 2.0);
+    return 0.0;
+  }
+
+  VectorField CircleVectorField(const Point &x) const {
+    return 2 * Pi * VectorField(0.5 - x[1], x[0] - 0.5, 0.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return CircleVectorField(x);
+  };
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return CircleVectorField(x) * N;
+  };
+
+  string Name() const override { return "CosHat2D"; }
+};
+
 class StochasticGaussHat2D : public IStochasticTransportProblem {
 private:
   double mu = 0.0;

mluq/tests/pdesolvers/TestPDESolver.cpp

 #include "TestPDESolver.hpp"
 
 
+/*
+ * Todo : Tests for CommSplit
+ */
+
 PDESOLVERTEST_TESTRUN(TestLaplace1D)
 
 PDESOLVERTEST_TESTRUN(TestLaplace2D)
 
 //PDESOLVERTEST_TESTRUN(TestTransport1D)
 
-//PDESOLVERTEST_TESTRUN(TestTransport2D)
+PDESOLVERTEST_TESTRUN(TestTransport2D)
+
 
 int main(int argc, char **argv) {
     return MppTest(

mluq/tests/pdesolvers/TestPDESolver.hpp

@@ -10,159 +10,144 @@
 constexpr double PDESOLVER_TEST_TOLERANCE = 1e-3;
 
 struct TestParams {
-    std::string model;
-    std::string problem;
-    std::string quantity;
-    double Q;
-
-    std::string GetOverlap() const {
-        if (model.find("DG") != -1) return "dG1";
-        else return "NoOverlap";
-    }
-
-    int GetDegree() const {
-        if (model.find("Laplace") != -1) return 1;
-        else return 2;
-    }
+  std::string model;
+  std::string problem;
+  std::string quantity;
+  double Q;
+
+  std::string GetOverlap() const {
+    if (model.find("DG") != -1) return "dG1";
+    else return "NoOverlap";
+  }
+
+  int GetDegree() const {
+    if (model.find("Laplace") != -1) return 1;
+    else return 2;
+  }
 };
 
 std::ostream &operator<<(std::ostream &s, const TestParams &testParams) {
-    return s << "{Model: " << testParams.model
-             << ", Problem: " << testParams.problem
-             << ", Quantity: " << testParams.quantity << "}" << endl;
+  return s << "{Model: " << testParams.model
+           << ", Problem: " << testParams.problem
+           << ", Quantity: " << testParams.quantity << "}" << endl;
 }
 
 class TestPDESolver : public TestWithParam<TestParams> {
 protected:
-    SampleID id;
-
-    Meshes *meshes;
-
-    PDESolver *pdeSolver;
-
-    TestPDESolver(const std::string &meshName, int commSplit = 0) {
-        std::filesystem::remove_all("data/vtk/");
-        std::filesystem::create_directory("data/vtk/");
-
-        id = SampleID(4, 0, false);
-
-        meshes = MeshesCreator(meshName).
-            WithOverlap(GetParam().GetOverlap()).
-            WithPLevel(id.cLevel).
-            WithLevel(id.fLevel).
-            WithCommSplit(commSplit).
-            WithDistribute("RCB").
-            Create();
-
-        pdeSolver = PDESolverCreator().
-            WithDegree(GetParam().GetDegree()).
-            WithQuantity(GetParam().quantity).
-            WithProblem(GetParam().problem).
-            WithModel(GetParam().model).
-            Create(*meshes);
-    }
-
-    void TearDown() override {
-        Plotting::InstanceM().Clear();
-        delete pdeSolver;
-        delete meshes;
-    }
+  SampleID id;
+
+  std::unique_ptr<Meshes> meshes;
+
+  std::unique_ptr<PDESolver> pdeSolver;
+
+  TestPDESolver(const std::string &meshName, int commSplit = 0) {
+    std::filesystem::remove_all("data/vtk/");
+    std::filesystem::create_directory("data/vtk/");
+
+    id = SampleID(4, 0, false);
+
+    meshes = MeshesCreator(meshName).
+      WithOverlap(GetParam().GetOverlap()).
+      WithPLevel(id.cLevel).
+      WithLevel(id.fLevel).
+      WithCommSplit(commSplit).
+      WithDistribute("RCB").
+      CreateUnique();
+
+    pdeSolver = PDESolverCreator().
+      WithDegree(GetParam().GetDegree()).
+      WithQuantity(GetParam().quantity).
+      WithProblem(GetParam().problem).
+      WithModel(GetParam().model).
+      CreateUnique(*meshes);
+  }
+
+  void TearDown() override {
+    PPM->Barrier(0);
+    Plotting::InstanceM().Clear();
+  }
 };
 
 class TestLaplace1D : public TestPDESolver {
 public:
-    TestLaplace1D() : TestPDESolver("Interval") {};
+  TestLaplace1D() : TestPDESolver("Interval") {};
 };
 
 INSTANTIATE_TEST_SUITE_P(TestPDESolver, TestLaplace1D, Values(
-    TestParams{"LagrangeElliptic", "Laplace1D", "L2Error", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace1D", "EnergyError", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace1D", "FluxError", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace1D", "Inflow", -1.0},
-    TestParams{"LagrangeElliptic", "Laplace1D", "Outflow", 1.0},
-
-    TestParams{"MixedElliptic", "Laplace1D", "L2CellAverageError", 0.0},
-    TestParams{"MixedElliptic", "Laplace1D", "EnergyError", 0.0},
-    TestParams{"MixedElliptic", "Laplace1D", "FluxError", 0.0},
-    TestParams{"MixedElliptic", "Laplace1D", "Inflow", -1.0},
-    TestParams{"MixedElliptic", "Laplace1D", "Outflow", 1.0},
-
-    TestParams{"HybridElliptic", "Laplace1D", "L2CellAverageError", 0.0},
-    TestParams{"HybridElliptic", "Laplace1D", "Inflow", -1.0},
-    TestParams{"HybridElliptic", "Laplace1D", "Outflow", 1.0},
-    TestParams{"HybridElliptic", "Laplace1D", "EnergyError", 0.0},
-    TestParams{"HybridElliptic", "Laplace1D", "FluxError", 0.0},
-
-    TestParams{"DGElliptic", "Laplace1D", "L2Error", 0.0},
-    TestParams{"DGElliptic", "Laplace1D", "EnergyError", 0.0},
-    TestParams{"DGElliptic", "Laplace1D", "FluxError", 0.0},
-    TestParams{"DGElliptic", "Laplace1D", "Inflow", -1.0},
-    TestParams{"DGElliptic", "Laplace1D", "Outflow", 1.0}
+  TestParams{"LagrangeElliptic", "Laplace1D", "L2Error", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace1D", "EnergyError", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace1D", "FluxError", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace1D", "Inflow", -1.0},
+  TestParams{"LagrangeElliptic", "Laplace1D", "Outflow", 1.0},
+
+  TestParams{"MixedElliptic", "Laplace1D", "L2CellAverageError", 0.0},
+  TestParams{"MixedElliptic", "Laplace1D", "EnergyError", 0.0},
+  TestParams{"MixedElliptic", "Laplace1D", "FluxError", 0.0},
+  TestParams{"MixedElliptic", "Laplace1D", "Inflow", -1.0},
+  TestParams{"MixedElliptic", "Laplace1D", "Outflow", 1.0},
+
+  TestParams{"HybridElliptic", "Laplace1D", "L2CellAverageError", 0.0},
+  TestParams{"HybridElliptic", "Laplace1D", "Inflow", -1.0},
+  TestParams{"HybridElliptic", "Laplace1D", "Outflow", 1.0},
+  TestParams{"HybridElliptic", "Laplace1D", "EnergyError", 0.0},
+  TestParams{"HybridElliptic", "Laplace1D", "FluxError", 0.0},
+
+  TestParams{"DGElliptic", "Laplace1D", "L2Error", 0.0},
+  TestParams{"DGElliptic", "Laplace1D", "EnergyError", 0.0},
+  TestParams{"DGElliptic", "Laplace1D", "FluxError", 0.0},
+  TestParams{"DGElliptic", "Laplace1D", "Inflow", -1.0},
+  TestParams{"DGElliptic", "Laplace1D", "Outflow", 1.0}
 ));
 
 class TestLaplace2D : public TestPDESolver {
 public:
-    TestLaplace2D() : TestPDESolver("Square") {};
+  TestLaplace2D() : TestPDESolver("Square") {};
 };
 
 INSTANTIATE_TEST_SUITE_P(TestPDESolver, TestLaplace2D, Values(
-    TestParams{"LagrangeElliptic", "Laplace2D", "L2Error", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace2D", "EnergyError", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace2D", "FluxError", 0.0},
-    TestParams{"LagrangeElliptic", "Laplace2D", "Inflow", -1.0},
-    TestParams{"LagrangeElliptic", "Laplace2D", "Outflow", 1.0},
-
-    TestParams{"MixedElliptic", "Laplace2D", "L2CellAverageError", 0.0},
-    TestParams{"MixedElliptic", "Laplace2D", "EnergyError", 0.0},
-    TestParams{"MixedElliptic", "Laplace2D", "FluxError", 0.0},
-    TestParams{"MixedElliptic", "Laplace2D", "Inflow", -1.0},
-    TestParams{"MixedElliptic", "Laplace2D", "Outflow", 1.0},
-
-    TestParams{"HybridElliptic", "Laplace2D", "L2CellAverageError", 0.0},
-    TestParams{"HybridElliptic", "Laplace2D", "EnergyError", 0.0},
-    TestParams{"HybridElliptic", "Laplace2D", "FluxError", 0.0},
-    TestParams{"HybridElliptic", "Laplace2D", "Inflow", -1.0},
-    TestParams{"HybridElliptic", "Laplace2D", "Outflow", 1.0},
-
-    TestParams{"DGElliptic", "Laplace2D", "L2Error", 0.0},
-    TestParams{"DGElliptic", "Laplace2D", "EnergyError", 0.0},
-    TestParams{"DGElliptic", "Laplace2D", "FluxError", 0.0},
-    TestParams{"DGElliptic", "Laplace2D", "Inflow", -1.0},
-    TestParams{"DGElliptic", "Laplace2D", "Outflow", 1.0}
+  TestParams{"LagrangeElliptic", "Laplace2D", "L2Error", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace2D", "EnergyError", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace2D", "FluxError", 0.0},
+  TestParams{"LagrangeElliptic", "Laplace2D", "Inflow", -1.0},
+  TestParams{"LagrangeElliptic", "Laplace2D", "Outflow", 1.0},
+
+  TestParams{"MixedElliptic", "Laplace2D", "L2CellAverageError", 0.0},
+  TestParams{"MixedElliptic", "Laplace2D", "EnergyError", 0.0},
+  TestParams{"MixedElliptic", "Laplace2D", "FluxError", 0.0},
+  TestParams{"MixedElliptic", "Laplace2D", "Inflow", -1.0},
+  TestParams{"MixedElliptic", "Laplace2D", "Outflow", 1.0},
+
+  TestParams{"HybridElliptic", "Laplace2D", "L2CellAverageError", 0.0},
+  TestParams{"HybridElliptic", "Laplace2D", "EnergyError", 0.0},
+  TestParams{"HybridElliptic", "Laplace2D", "FluxError", 0.0},
+  TestParams{"HybridElliptic", "Laplace2D", "Inflow", -1.0},
+  TestParams{"HybridElliptic", "Laplace2D", "Outflow", 1.0},
+
+  TestParams{"DGElliptic", "Laplace2D", "L2Error", 0.0},
+  TestParams{"DGElliptic", "Laplace2D", "EnergyError", 0.0},
+  TestParams{"DGElliptic", "Laplace2D", "FluxError", 0.0},
+  TestParams{"DGElliptic", "Laplace2D", "Inflow", -1.0},
+  TestParams{"DGElliptic", "Laplace2D", "Outflow", 1.0}
 ));
 
-//class TestTransport1D : public TestPDESolver {
-//public:
-//    TestTransport1D() : TestPDESolver("Interval") {};
-//};
-//
-//INSTANTIATE_TEST_SUITE_P(TestPDESolver, TestTransport1D, Values(
-//    TestParams{"DGTransport", "Pollution1D", "Mass", 1.0}
-//));
-//
-//TEST_P(TestTransport1D, TestRun) {
-//    SampleSolution solution(pdeSolver->MGraphs(), id);
-//    pdeSolver->Run(solution);
-//    EXPECT_NEAR(solution.Q, GetParam().Q, PDESOLVER_TEST_TOLERANCE);
-//}
+class TestTransport1D : public TestPDESolver {
+public:
+    TestTransport1D() : TestPDESolver("Interval") {};
+};
+
+INSTANTIATE_TEST_SUITE_P(TestPDESolver, TestTransport1D, Values(
+    TestParams{"DGTransport", "CosHat1D", "Mass", 1.0}
+));
 
 class TestTransport2D : public TestPDESolver {
 public:
-    TestTransport2D() : TestPDESolver("Square") {};
+  TestTransport2D() : TestPDESolver("Square") {};
 };
 
 INSTANTIATE_TEST_SUITE_P(TestPDESolver, TestTransport2D, Values(
-//    TestParams{"DGTransport", "Pollution2D", "Mass", 1.0}
-//    TestParams{"DGTransport", "PollutionCosHat2D", "Mass", 1.0},
-    TestParams{"DGTransport", "GaussHat2D", "Mass", 1.0}
+  TestParams{"DGTransport", "CosHat2D", "Mass", 0.036012725}
 ));
 
-//TEST_P(TestTransport2D, TestInitialValue) {
-////    SampleSolution solution(pdeSolver->MGraphs(), id);
-////    pdeSolver->Run(solution);
-////    EXPECT_NEAR(solution.Q, GetParam().Q, PDESOLVER_TEST_TOLERANCE);
-//}
-
 #define PDESOLVERTEST_TESTRUN(TestClass)\
 \
 TEST_P(TestClass, TestRun) {\