moved time information to problem, bug fixes and problem refactoring

mlmc/src/problems/StochasticTransportProblem.hpp

@@ -8,9 +8,25 @@ class IStochasticTransportProblem : public IStochasticProblem {
 private:
     bool rhs = true;
 
+    double CFL = 0.5;
+
+    double t0 = 0.0;
+
+    double T = 1.0;
+
 public:
     IStochasticTransportProblem(Meshes &meshes, GeneratorNames genNames = {})
-        : IStochasticProblem(meshes, genNames) {}
+        : IStochasticProblem(meshes, genNames) {
+        config.get("CFL", CFL);
+        config.get("t0", t0);
+        config.get("T", T);
+    }
+
+    double GetStartTime() { return t0; }
+
+    double GetEndTime() { return T; }
+
+    double GetCFL() { return CFL; }
 
     bool RHS() const { return rhs; }
 
@@ -135,13 +151,13 @@ public:
     Pollution2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
 
     Scalar Solution(double t, const Point &x) const override {
-        if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
-            return 1.0;
+        if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 4 * d))
+            return 32.0;
         return 0.0;
     }
 
     VectorField TransportFlux(const Point &x) const {
-        return VectorField(-1.0, 0.0, 0.0);
+        return VectorField(0.0, -1.0);
     }
 
     VectorField CellFlux(const cell &c, const Point &x) const override {
@@ -185,16 +201,15 @@ public:
 };
 
 class PollutionCosHat2D : public IStochasticTransportProblem {
-    double amplitude = 1.00;
-    double cc = 0.5;
+    double amplitude = 2.0 / Pi;
 public:
     explicit PollutionCosHat2D(Meshes &meshes) : IStochasticTransportProblem(meshes) {}
 
     double Solution(double t, const Point &x) const override {
-        Point midPoint = 5.0 * Point(cos(2. * Pi * t), sin(2. * Pi * t));
+        Point midPoint = 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t))
+            + Point(1.0, 1.0));
         double r = dist(midPoint, x);
-        if (r < 1 / cc)
-            return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
+        if (r < 1.0 / (2.0 * Pi)) return amplitude * pow(cos(r / 2.0), 2.0);
         return 0.0;
     }
 
@@ -233,13 +248,17 @@ public:
         return 1.0;
     }
 
+    VectorField CircleVectorField(const Point &x) const {
+        return 2 * Pi * VectorField(-x[1], x[0], 0.0);
+    }
+
     VectorField CellFlux(const cell &c, const Point &x) const override {
-        return vectorField;
+        return CircleVectorField(x);
     }
 
     Scalar FaceNormalFlux(const cell &c, int face,
                           const VectorField &N, const Point &x) const override {
-        return vectorField * N;
+        return CircleVectorField(x) * N;
     }
 
     string Name() const override { return "GaussHat2D"; }
@@ -247,38 +266,42 @@ public:
 
 class GaussHat2D : public IStochasticTransportProblem {
 private:
-    VectorField vectorField;
-    VectorField mu;
     Tensor sigma;
 
 public:
     GaussHat2D(Meshes &meshes) :
-        mu(VectorField(0.5, 0.5, 0.0)),
         sigma(Tensor()),
-        vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
         IStochasticTransportProblem(meshes) {
-        sigma[0][0] = 0.01;
+        sigma[0][0] = 0.005;
         sigma[0][1] = 0.0;
         sigma[1][0] = 0.0;
-        sigma[1][1] = 0.01;
-        sigma[2][2] = 1.0; // Todo might not be working for 2D Point
+        sigma[1][1] = 0.005;
+//        sigma[2][2] = 1.0; // Todo might not be working for 2D Point
+    }
+
+    Point mu(double t) const {
+        return 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t)) + Point(1.0, 1.0));
     }
 
     Scalar Solution(double t, const Point &x) const override {
         double factor = 1 / sqrt(pow(2 * Pi, 2) * sigma.det());
-        VectorField diff = (x - mu);
+        VectorField diff = (x - mu(t));
         VectorField temp = Invert(sigma) * diff;
         Scalar exponent = -0.5 * diff * temp;
         return factor * exp(exponent);
     }
 
+    VectorField CircleVectorField(const Point &x) const {
+        return 2 * Pi * VectorField(-x[1], x[0], 0.0);
+    }
+
     VectorField CellFlux(const cell &c, const Point &x) const override {
-        return vectorField;
+        return CircleVectorField(x);
     }
 
     Scalar FaceNormalFlux(const cell &c, int face,
                           const VectorField &N, const Point &x) const override {
-        return vectorField * N;
+        return CircleVectorField(x) * N;
     }
 
     string Name() const override { return "GaussHat2D"; }