From ae9d1a422b2122935cfc50aa5a62b4fb158d1a3f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Tue, 9 Jun 2020 23:50:10 +0200
Subject: [PATCH] add problem dg4
---
dune/mmdg/problems/dgproblem4.hh | 166 +++++++++++++++++++++++++++++++
src/dg.cc | 18 +++-
src/parameterDG.ini | 4 +-
3 files changed, 186 insertions(+), 2 deletions(-)
create mode 100644 dune/mmdg/problems/dgproblem4.hh
diff --git a/dune/mmdg/problems/dgproblem4.hh b/dune/mmdg/problems/dgproblem4.hh
new file mode 100644
index 0000000..d00936a
--- /dev/null
+++ b/dune/mmdg/problems/dgproblem4.hh
@@ -0,0 +1,166 @@
+#ifndef __DUNE_MMDG_DGPROBLEM4_HH__
+#define __DUNE_MMDG_DGPROBLEM4_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+#include <cmath>
+
+template<class Coordinate, class Scalar = double>
+class DGProblem4 : public DGProblem<Coordinate, Scalar>
+{
+ public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = DGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ // static constexpr Scalar dmin_ = 0.05; //minimum aperture of the fracture
+ // static constexpr Scalar dmax_ = 0.2; //maximum aperture of the fracture
+
+ //require dim == 2
+ DGProblem4 (const Scalar dmin, const Scalar dmax)
+ : dmin_(dmin), dmax_(dmax) {}
+ //{
+ // if (dim != 2)
+ // {
+ // DUNE_THROW(Dune::Exception,
+ // "Problem dg4 is not implemented for dimension = " +
+ // std::to_string(dim) + "!");
+ // }
+ //}
+
+ //the exact solution at position pos
+ Scalar exactSolution (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ if ( pos[0] < 0.5 * ( 1.0 - d ) )
+ { //pos in Omega_1 (left of the fracture)
+ return 2.0 * ( pos[0] + d ) * exp(-d);
+ }
+
+ if ( pos[0] > 0.5 * ( 1.0 + d ) )
+ { //pos in Omega_2 (right of the fracture)
+ return 2.0 * pos[0] * exp(d);
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return ( 1.0 + d ) * exp( 2.0 * pos[0] - 1.0 );
+ }
+
+ //exact solution is implemented
+ bool hasExactSolution() const
+ {
+ return true;
+ };
+
+ //source term at position pos
+ Scalar q (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+ Scalar dPrime2 = aperturePrime(pos);
+ dPrime2 *= dPrime2;
+ const Scalar dPrimePrime = aperturePrimePrime(pos);
+
+
+ if ( pos[0] < 0.5 * ( 1.0 - d ) )
+ { //pos in Omega_1 (left of the fracture)
+ return -2.0 * exp(-d) / ( 1.0 - d ) * ( ( 1.0 - pos[0] - d )
+ * ( dPrimePrime + dPrime2 * d / ( 1.0 - d ) ) - dPrime2 );
+ }
+
+ if ( pos[0] > 0.5 * ( 1.0 + d ) )
+ { //pos in Omega_2 (right of the fracture)
+ return -2.0 * pos[0] * exp(d) / ( 1.0 + d )
+ * ( dPrimePrime + dPrime2 * d / ( 1.0 + d ) );
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return
+ -exp(2.0 * pos[0] - 1) * ( 4.0 + dPrimePrime );
+ }
+
+ //permeability tensor at position pos
+ Matrix K (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ Matrix permeability(0.0);
+
+ if ( pos[0] < 0.5 * ( 1.0 - d ) )
+ { //pos in Omega_1 (left of the fracture)
+ permeability[0][0] = 1.0;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = 1.0 / ( 1.0 - d );
+ }
+ }
+ else if ( pos[0] > 0.5 * ( 1.0 + d ) )
+ { //pos in Omega_2 (right of the fracture)
+ permeability[0][0] = 1.0;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = 1.0 / ( 1.0 + d );
+ }
+ }
+ else
+ { //pos in Omega_R (inside the fracture)
+ permeability[0][0] = 1.0 / ( 1.0 + d );
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = 1.0;
+ }
+ }
+
+ return permeability;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * q(x)
+ virtual int quadratureOrder_q () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * boundary(x) and K(x) * boundary(x)
+ virtual int quadratureOrderBoundary () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over an entry K[i][j] of the permeability tensor
+ virtual int quadratureOrder_K () const
+ {
+ return 10;
+ }
+
+ private:
+ //aperture of the fracture at position pos
+ const Scalar aperture (const Coordinate& pos) const
+ {
+ return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
+ ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ //derivative of aperture at position pos
+ const Scalar aperturePrime (const Coordinate& pos) const
+ {
+ return -4.0 * M_PI * ( dmax_ - dmin_ ) * std::sin(8.0 * M_PI * pos[1]);
+ }
+
+ //second derivative of aperture at position pos
+ const Scalar aperturePrimePrime (const Coordinate& pos) const
+ {
+ return -32.0 * M_PI * M_PI * ( dmax_ - dmin_ )
+ * std::cos(8.0 * M_PI * pos[1]);
+ }
+
+ const Scalar dmin_; //minimum aperture of the fracture
+ const Scalar dmax_; //maximum aperture of the fracture
+
+};
+
+#endif
diff --git a/src/dg.cc b/src/dg.cc
index 2dedbca..5304bf3 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -24,6 +24,7 @@
#include <dune/mmdg/problems/dgproblem1.hh>
#include <dune/mmdg/problems/dgproblem2.hh>
#include <dune/mmdg/problems/dgproblem3.hh>
+#include <dune/mmdg/problems/dgproblem4.hh>
int main(int argc, char** argv)
{
@@ -85,11 +86,26 @@ int main(int argc, char** argv)
"an infinite series; increase the cut-off if necessary!\n\n";
}
}
+ else if (pt["problem"] == "dg4")
+ {
+ if ( pt.hasKey("dmin") && pt.hasKey("dmax") )
+ {
+ const double dmin = std::stod(pt["dmin"]);
+ const double dmax = std::stod(pt["dmax"]);
+ problem = new DGProblem4<Coordinate>(dmin, dmax);
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg4 requires the parameters 'dmin' and 'dmax'. Check the "
+ "file parameterDG.ini and make sure that they are specified.");
+ }
+ }
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'dg1' to 'dg3')");
+ "parameterDG.ini (use 'dg1' to 'dg4')");
}
DG dg(gridView, mapper, *problem);
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 275babc..f6448b0 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,3 +1,5 @@
mu0 = 1000
-problem = dg2
+problem = dg4
fastEOC = true
+dmin = 1e-1
+dmax = 3e-1
--
GitLab