change problem mmdg6 to allow arbitrary xi

dune/mmdg/problems/mmdgproblem6.hh

@@ -12,7 +12,10 @@ public:
   static constexpr int dim = Coordinate::dimension;
 
   //constructor
-  MMDGProblem6 (const Scalar d0) : d0_(d0) {}
+  MMDGProblem6 (const Scalar d0) : d0_(d0), xi_(1.0) {}
+
+  //constructor
+  MMDGProblem6 (const Scalar d0, const Scalar xi) : d0_(d0), xi_(xi) {}
 
   //the exact bulk solution at position pos
   Scalar exactSolution (const Coordinate& pos) const
@@ -26,7 +29,11 @@ public:
   //the exact solution on the interface at position pos
   Scalar exactInterfaceSolution (const Coordinate& pos) const
   {
-    return std::cos(4.0 * pos[1]);
+    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
+    constexpr Scalar sinhsqrt2 = std::sinh(sqrt(2.0));
+
+    return std::cos(4.0 * pos[1]) * (coshsqrt2 * coshsqrt2
+      - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
   }
 
   //indicates whether an exact solution is implemented for the problem
@@ -49,8 +56,12 @@ public:
   {
     const Scalar cos4y = std::cos(4.0 * pos[1]);
     constexpr Scalar sqrt2 = 2.0 * sqrt(2.0);
+    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
+    constexpr Scalar sinhsqrt2 = std::sinh(sqrt(2.0));
+    const Scalar hyperbolicFactor =
+      (coshsqrt2 * coshsqrt2 - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
 
-    return -48.0 * d0_ * d0_ * std::exp(-8.0 * pos[1]) * ( cos4y
+    return -48.0 * d0_ * d0_ * std::exp(-8.0 * pos[1]) * hyperbolicFactor * ( cos4y
       + 3.0 * std::sin(4.0 * pos[1]) ) - sqrt2 * std::sinh(sqrt2) * cos4y;
   }
 
@@ -117,6 +128,7 @@ public:
 
 private:
   Scalar d0_; //prefactor of the aperture
+  Scalar xi_; //coupling constant
 
 };
 

src/mmdg.cc

@@ -95,9 +95,9 @@ int main(int argc, char** argv)
     }
     else if (pt["problem"] == "mmdg6")
     {
-      problem = new MMDGProblem6<Coordinate>(aperture);
+      xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
+      problem = new MMDGProblem6<Coordinate>(aperture, xi);
       gridType = "mmdg2_5e-2";
-      xi = 1.0;
     }
     else
     {