diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index ff022d9742290c56a2fc674a861b4860ac9ca611..e45f08dbf889e9955da3cc97e45a0025c5ab230c 100644
--- a/dune/mmdg/problems/coupleddgproblem.hh
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -3,7 +3,7 @@
#include <dune/mmdg/problems/dgproblem.hh>
-//abstract class providing an interface for a coupled flow problem in a
+//abstract class that provides an interface for a coupled flow problem in a
//fractured porous medium that is to be solved with a coupled discontinuous
//Galerkin scheme
template<class Coordinate, class Scalar = double>
@@ -38,12 +38,13 @@ public:
return Coordinate(0.0);
}
+ //interface boundary condition at position pos
virtual Scalar interfaceBoundary(const Coordinate& pos) const
{
return exactInterfaceSolution(pos);
}
- //tangential permeability per aperture tensor of the interface at position pos
+ //tangential permeability per aperture of the interface at position pos
virtual Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
@@ -68,7 +69,7 @@ public:
//over x * qInterface(x)
virtual int quadratureOrder_qInterface () const
{
- return 10;//1;
+ return 1;
}
//returns the recommended quadrature order to compute an integral
@@ -77,7 +78,7 @@ public:
//and d^2 * interfaceBoundary(x) * Kpar
virtual int quadratureOrderInterfaceBoundary () const
{
- return 10;//1;
+ return 1;
}
//returns the recommended quadrature order to compute integrals
@@ -87,14 +88,14 @@ public:
//and x * d^2 * Kpar
virtual int quadratureOrder_Kparallel () const
{
- return 10;//2;
+ return 2;
}
//returns the recommended quadrature order to compute an integral
//over x * Kperp(x)
virtual int quadratureOrder_Kperp () const
{
- return 10;//1;
+ return 1;
}
};
diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh
index 38ed3077b35279cf29c0d64ba1e5fc92c46d603d..4cd466905c396d72fb8ee6de7a088ecef07d414e 100644
--- a/dune/mmdg/problems/dgproblem.hh
+++ b/dune/mmdg/problems/dgproblem.hh
@@ -1,8 +1,8 @@
#ifndef __DUNE_MMDG_DGPROBLEM_HH__
#define __DUNE_MMDG_DGPROBLEM_HH__
-//abstract class providing an interface for a problem that is to be solved with
-//a discontinuous Galerkin scheme
+//abstract class that provides an interface for a problem that is to be solved
+//with a discontinuous Galerkin scheme
template<class Coordinate, class Scalar = double>
class DGProblem
{
diff --git a/dune/mmdg/problems/dgproblem1.hh b/dune/mmdg/problems/dgproblem1.hh
index 39f67481129a76e9217813072d0bf0cf96c54d10..52ca8ea29375e169269508989738f4ac05a0152d 100644
--- a/dune/mmdg/problems/dgproblem1.hh
+++ b/dune/mmdg/problems/dgproblem1.hh
@@ -11,9 +11,9 @@ class DGProblem1 : public DGProblem<Coordinate, Scalar>
static constexpr int dim = Coordinate::dimension;
//the exact solution at position pos
- //1d: f(x) = 0.5*x^2,
- //2d: f(x,y) = 0.5*x^2*y^2,
- //3d: f(x,y,z) = 0.5*x^2*y^2*z^2
+ //1d: f(x) = 0.5 * x^2,
+ //2d: f(x,y) = 0.5 * x^2 * y^2,
+ //3d: f(x,y,z) = 0.5 * x^2 * y^2 * z^2
Scalar exactSolution (const Coordinate& pos) const
{
Scalar solution(0.5);
@@ -46,7 +46,7 @@ class DGProblem1 : public DGProblem<Coordinate, Scalar>
}
if constexpr (dim == 3)
- { //q(x,y,z) = -x^2*y^2 - x^2*z^2 - y^2*z^2
+ { //q(x,y,z) = -x^2 * y^2 - x^2 * z^2 - y^2 * z^2
Scalar pos1_2 = pos[1] * pos[1];
Scalar pos2_2 = pos[2] * pos[2];
return -pos[0] * pos[0] * (pos1_2 + pos2_2) - pos1_2 * pos2_2;
diff --git a/dune/mmdg/problems/dgproblem2.hh b/dune/mmdg/problems/dgproblem2.hh
index 5d8240452bc414b30a6f886b73bdc64d0068709e..96415cd3611b1817fb211f33f852d975ba3bbfa8 100644
--- a/dune/mmdg/problems/dgproblem2.hh
+++ b/dune/mmdg/problems/dgproblem2.hh
@@ -3,6 +3,7 @@
#include <dune/mmdg/problems/dgproblem.hh>
+//a test problem with anisotropic permeability but without fracture
template<class Coordinate, class Scalar = double>
class DGProblem2 : public DGProblem<Coordinate, Scalar>
{
@@ -17,7 +18,6 @@ class DGProblem2 : public DGProblem<Coordinate, Scalar>
//3d: p(x,y,z) = x^2 + y^2 + z^2
Scalar exactSolution (const Coordinate& pos) const
{
- //return pos * Coordinate(1.0);
return pos * pos;
}
diff --git a/dune/mmdg/problems/dgproblem3.hh b/dune/mmdg/problems/dgproblem3.hh
index 21659d847274b1984e4a126d6f1446e68b896894..1ba4e5d6816f178fa835590b8e6bfd847d95e138 100644
--- a/dune/mmdg/problems/dgproblem3.hh
+++ b/dune/mmdg/problems/dgproblem3.hh
@@ -2,11 +2,9 @@
#define __DUNE_MMDG_DGPROBLEM3_HH__
#include <cmath>
-#include <string>
-
#include <dune/mmdg/problems/dgproblem.hh>
-//Poisson problem with right side q = -1
+//Poisson problem with right side q = -1 and homogeneous Dirichlet boundary
template<class Coordinate, class Scalar = double>
class DGProblem3 : public DGProblem<Coordinate, Scalar>
{
@@ -76,13 +74,13 @@ class DGProblem3 : public DGProblem<Coordinate, Scalar>
//source term
Scalar q (const Coordinate& pos) const
{
- return Scalar(-1.0);
+ return -1.0;
}
//boundary condition at position pos
- virtual Scalar boundary (const Coordinate& pos) const
+ Scalar boundary (const Coordinate& pos) const
{
- return Scalar(0.0);
+ return 0.0;
}
};
diff --git a/dune/mmdg/problems/dgproblem4.hh b/dune/mmdg/problems/dgproblem4.hh
index 7ea27cbd877ee3567c832b4eb1c86697da6047b4..07b8dd0aea89b1e81d52ece3c0d719337f8bb365 100644
--- a/dune/mmdg/problems/dgproblem4.hh
+++ b/dune/mmdg/problems/dgproblem4.hh
@@ -4,6 +4,8 @@
#include <dune/mmdg/problems/dgproblem.hh>
#include <cmath>
+//a problem with cosine-shaped fracture and continuous pressure and Darcy
+//velocity
template<class Coordinate, class Scalar = double>
class DGProblem4 : public DGProblem<Coordinate, Scalar>
{
@@ -12,7 +14,7 @@ class DGProblem4 : public DGProblem<Coordinate, Scalar>
using Base = DGProblem<Coordinate, Scalar>;
using Matrix = typename Base::Matrix;
- //require dim != 1
+ //constructor, require dim != 1
DGProblem4 (const Scalar dmin, const Scalar dmax) : dmin_(dmin), dmax_(dmax)
{
if (dim == 1)
diff --git a/dune/mmdg/problems/dgproblem5.hh b/dune/mmdg/problems/dgproblem5.hh
index 880f6b9af0effefc16c784282acfc0c5b33469f8..678c0e5283d87b8e6db6558ffc6362a31c65a397 100644
--- a/dune/mmdg/problems/dgproblem5.hh
+++ b/dune/mmdg/problems/dgproblem5.hh
@@ -4,6 +4,11 @@
#include <dune/mmdg/problems/dgproblem.hh>
#include <cmath>
+//a problem with a cosine-shaped fracture and continuous pressure
+//NOTE: This problem does not show convergence as the normal Darcy velocity is
+// only continuous with respect to the normal of the hyperplane
+// Gamma = { x_1 = 0.5 } but not with respect to the normal of the actual
+// interfaces between the bulk and fracture domains
template<class Coordinate, class Scalar = double>
class DGProblem5 : public DGProblem<Coordinate, Scalar>
{
@@ -16,17 +21,16 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
static constexpr Scalar K_par = 3.0; //tangential fracture permeability
static constexpr Scalar K_perp = 0.5; //normal fracture permeability
- //constructor
- DGProblem5 (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) + "!");
- // }
- // }
+ //constructor, require dim != 1
+ DGProblem5 (const Scalar dmin, const Scalar dmax) : dmin_(dmin), dmax_(dmax)
+ {
+ if (dim == 1)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg5 is not implemented for dimension = " +
+ std::to_string(dim) + "!");
+ }
+ }
//the exact solution at position pos
Scalar exactSolution (const Coordinate& pos) const
@@ -120,19 +124,6 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
( 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
};
diff --git a/dune/mmdg/problems/mmdgproblem1.hh b/dune/mmdg/problems/mmdgproblem1.hh
index e1caa9b773a8869f508b410b0286097b6235eee2..d581459eceb434ce520f374cb86d6f64c1021f32 100644
--- a/dune/mmdg/problems/mmdgproblem1.hh
+++ b/dune/mmdg/problems/mmdgproblem1.hh
@@ -4,7 +4,12 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem
+//a coupled dg problem wiht non-constant aperture
+//NOTE: To deal with the non-constant aperture, an adaption of the coupled DG
+// scheme implemented in mmdg.hh is used to solve this problem. Thus, the
+// interface source term and the tangential permeability per aperture
+// inside the fracture are multiplied with the aperture function d in
+// the following problem definition
template<class Coordinate, class Scalar = double>
class MMDGProblem1 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -14,7 +19,7 @@ public:
using Matrix = typename Base::Matrix;
//constructor
- MMDGProblem1 (const Scalar d) : d_(d) {}
+ MMDGProblem1 (const Scalar dmax) : d_ ( dmax / 1.01 ) {}
//the exact bulk solution at position pos
Scalar exactSolution (const Coordinate& pos) const
@@ -28,7 +33,7 @@ public:
return 1.0;
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -38,55 +43,76 @@ public:
Scalar aperture (const Coordinate& pos) const
{
return d_ * (pos[1] + 0.01);
- return d_;
}
//interface source term at position pos
+ //NOTE: contains an extra factor d
Scalar qInterface (const Coordinate& pos) const
{
- const Scalar de = aperture(pos);
- const Scalar dPrime2 = gradAperture(pos) * gradAperture(pos);
- const Scalar dPrimePrime = 0.0;
+ const Scalar d = aperture(pos);
+ const Coordinate gradD = gradAperture(pos);
- return (-de * dPrimePrime - dPrime2) * de;
-
- return Scalar(0.0);
+ // q = -d * d'' - (d')^2 with extra factor d and d'' = 0
+ return -d * (gradD * gradD);
}
//tangential gradient of the aperture d of the fracture at position pos
Coordinate gradAperture (const Coordinate& pos) const
{
- Coordinate gradD(0.0);
- gradD[1] = d_;
- return gradD;
-
- return Coordinate(0.0);
+ Coordinate gradD(0.0);
+ gradD[1] = d_;
+ return gradD;
}
//tangential permeability per aperture tensor of the interface at position pos
+ //NOTE: contains an extra factor d
Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
for (int i = 0; i < dim; i++)
{
- permeability[i][i] = 1.0 * aperture(pos);
+ permeability[i][i] = aperture(pos);
}
- //return identity matrix
return permeability;
}
-
//returns the recommended quadrature order to compute an integral
//over x * boundary(x)
int quadratureOrderBoundary () const
{
- return 10;//2;
+ return 2;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * qInterface(x)
+ int quadratureOrder_qInterface () const
+ {
+ return 2;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
+ int quadratureOrderInterfaceBoundary () const
+ {
+ return 3;
+ }
+
+ //returns the recommended quadrature order to compute integrals
+ //over x^2 * (Kpar * grad(d) * grad(d))
+ //and x^2 * d^2
+ //and x^2 * d * (Kpar * grad(d))
+ //and x * d^2 * Kpar
+ int quadratureOrder_Kparallel () const
+ {
+ return 4;
}
private:
- Scalar d_; //aperture
+ Scalar d_; //precfactor of the aperture function
};
#endif
diff --git a/dune/mmdg/problems/mmdgproblem2.hh b/dune/mmdg/problems/mmdgproblem2.hh
index efe1f3f4e531caa7e1663cf412e5b6fa5fa79549..61624636efc2f67d55c37648f5c237f12aff14a1 100644
--- a/dune/mmdg/problems/mmdgproblem2.hh
+++ b/dune/mmdg/problems/mmdgproblem2.hh
@@ -4,8 +4,8 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem for dim = 2,
-//taken from [Antonietti et al. (2019): Example 6.3]
+//a coupled dg problem with constant aperture
+//taken with adaptions from [Antonietti et al. (2019): Example 6.3]
template<class Coordinate, class Scalar = double>
class MMDGProblem2 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -21,18 +21,18 @@ public:
Scalar exactSolution (const Coordinate& pos) const
{
Scalar solution =
- (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
- return solution * std::cos(M_PI * pos[1]);
+ (pos[0] < 0.5) ? sin(4 * pos[0]) : cos(4 * pos[0]);
+ return solution * cos(M_PI * pos[1]);
}
//the exact solution on the interface at position pos
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
- constexpr Scalar sourcefactor = 0.75 * (std::cos(2.0) + std::sin(2.0));
- return sourcefactor * std::cos(M_PI * pos[1]);
+ constexpr Scalar solFactor = 0.75 * ( cos(2.0) + sin(2.0) );
+ return solFactor * cos(M_PI * pos[1]);
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -43,29 +43,29 @@ public:
{
constexpr Scalar sourcefactor = 16.0 + M_PI * M_PI;
Scalar source =
- (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
- return sourcefactor * source * std::cos(M_PI * pos[1]);
+ (pos[0] < 0.5) ? sin(4 * pos[0]) : cos(4 * pos[0]);
+ return sourcefactor * source * cos(M_PI * pos[1]);
}
//interface source term at position pos
Scalar qInterface (const Coordinate& pos) const
{
- constexpr Scalar sourcefactor = std::cos(2.0) + std::sin(2.0);
- return sourcefactor * (4.0 + 0.75 * d_ * d_ * M_PI * M_PI)
- * std::cos(M_PI * pos[1]);
+ constexpr Scalar sourcefactor = cos(2.0) + sin(2.0);
+ return sourcefactor * ( 4.0 + 0.75 * d_ * d_ * M_PI * M_PI )
+ * cos(M_PI * pos[1]);
}
- //aperture d of the fracture at position pos
+ //constant aperture d of the fracture
Scalar aperture (const Coordinate& pos) const
{
return d_;
}
- //permeability per aperture of the fracture per aperture in normal direction
+ //permeability per aperture of the fracture per aperture in normal direction
//at position pos
Scalar Kperp (const Coordinate& pos) const
{
- return Scalar(2.0);
+ return 2.0;
}
//returns the recommended quadrature order to compute an integral
diff --git a/dune/mmdg/problems/mmdgproblem3.hh b/dune/mmdg/problems/mmdgproblem3.hh
index cf7f896a9aa6b61d19dc18a5cecf7d20c9b1b0fa..441aebc2d4ff3ae9b36b3b82a2ba26a47a9da7f2 100644
--- a/dune/mmdg/problems/mmdgproblem3.hh
+++ b/dune/mmdg/problems/mmdgproblem3.hh
@@ -4,8 +4,8 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem for dim = 2,
-//taken from [Antonietti et al. (2019): Example 6.1]
+//a coupled dg problem with constant aperture,
+//taken with adaptions from [Antonietti et al. (2019): Example 6.1]
template<class Coordinate, class Scalar = double>
class MMDGProblem3 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -21,18 +21,18 @@ public:
Scalar exactSolution (const Coordinate& pos) const
{
Scalar xPlusy = pos[0] + pos[1];
- Scalar solution = std::exp(xPlusy);
+ Scalar solution = exp(xPlusy);
return (xPlusy < 1.0) ? solution :
- 0.5 * solution + std::exp(1.0) * (0.5 + 1.5 / sqrt(2.0) * d_);
+ 0.5 * solution + exp(1.0) * (0.5 + 1.5 / sqrt(2.0) * d_);
}
//the exact solution on the interface at position pos
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
- return std::exp(1.0) * ( 1.0 + 1.0 / sqrt(2.0) * d_ );
+ return exp(1.0) * ( 1.0 + 1.0 / sqrt(2.0) * d_ );
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -42,18 +42,18 @@ public:
Scalar q (const Coordinate& pos) const
{
Scalar xPlusy = pos[0] + pos[1];
- Scalar source = -std::exp(xPlusy);
+ Scalar source = -exp(xPlusy);
return (xPlusy < 1.0) ? 2.0 * source : source;
}
//interface source term at position pos
Scalar qInterface (const Coordinate& pos) const
{
- constexpr Scalar iSource = std::exp(1.0) / sqrt(2.0);
+ constexpr Scalar iSource = exp(1.0) / sqrt(2.0);
return iSource;
}
- //aperture d of the fracture at position pos
+ //constant aperture d of the fracture
Scalar aperture (const Coordinate& pos) const
{
return d_;
diff --git a/dune/mmdg/problems/mmdgproblem4.hh b/dune/mmdg/problems/mmdgproblem4.hh
index 347d4096782730bf0c66a586636b4a8878b6d19f..c8adbfff494997ba68a2cad962f7da6bb05dcb07 100644
--- a/dune/mmdg/problems/mmdgproblem4.hh
+++ b/dune/mmdg/problems/mmdgproblem4.hh
@@ -4,18 +4,24 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem
+//a coupled dg problem wiht non-constant aperture
+//NOTE: To deal with the non-constant aperture, an adaption of the coupled DG
+// scheme implemented in mmdg.hh is used to solve this problem. Thus, the
+// interface source term and the tangential permeability per aperture
+// inside the fracture are multiplied with the aperture function d in
+// the following problem definition
template<class Coordinate, class Scalar = double>
class MMDGProblem4 : public CoupledDGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
- static constexpr Scalar jump = 1.0;
+ static constexpr Scalar jump = 1.0; //jump of the pressure at the interface
+
using Base = CoupledDGProblem<Coordinate, Scalar>;
using Matrix = typename Base::Matrix;
//constructor
- MMDGProblem4 (const Scalar d) : d_(d) {}
+ MMDGProblem4 (const Scalar dmax) : d_ (dmax / 1.01) {}
//the exact bulk solution at position pos
Scalar exactSolution (const Coordinate& pos) const
@@ -36,53 +42,53 @@ public:
return 0.25 - pos[1] * pos[1] + 0.5 * jump;
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
};
//interface source term at position pos
+ //NOTE: contains an extra factor d
Scalar qInterface (const Coordinate& pos) const
{
- // return 2.0 * d_ * d_;
- const Scalar de = aperture(pos);
+ const Scalar d = aperture(pos);
const Scalar dPrime = sqrt(gradAperture(pos) * gradAperture(pos));
- const Scalar dPrimePrime = 0.;
+ const Scalar dPrimePrime = 2.0 * d_;
- return (d_ * d_ * ( 10. * pos[1] * 0.01 + 9. * pos[1] * pos[1] + 2. * 0.01 * 0.01 - 0.5 * (jump + 0.5))) * de;
+ return d * ( -( dPrime * dPrime + d * dPrimePrime ) *
+ exactInterfaceSolution(pos) + 2.0 * d * ( 3.0 * pos[1] * dPrime + d ) );
}
//aperture d of the fracture at position pos
Scalar aperture (const Coordinate& pos) const
{
- return d_ * (pos[1] + 0.01);;
+ return d_ * ( pos[1] * pos[1] + 0.01 );
}
//tangential gradient of the aperture d of the fracture at position pos
Coordinate gradAperture (const Coordinate& pos) const
{
Coordinate gradD(0.0);
- gradD[1] = d_;
+ gradD[1] = 2.0 * d_ * pos[1];
return gradD;
-
- return Coordinate(0.0);
}
//permeability of the fracture in normal direction at position pos
Scalar Kperp (const Coordinate& pos) const
{
- return Scalar(1.0) / jump;
+ return 1.0 / jump;
}
//tangential permeability tensor of the interface at position pos
+ //NOTE: contains an extra factor d
Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
for (int i = 0; i < dim; i++)
{
- permeability[i][i] = 1.0 * aperture(pos);//d_;
+ permeability[i][i] = aperture(pos);
}
return permeability;
@@ -92,7 +98,14 @@ public:
//over x * boundary(x)
int quadratureOrderBoundary () const
{
- return 10;//3;
+ return 3;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * qInterface(x)
+ int quadratureOrder_qInterface () const
+ {
+ return 7;
}
//returns the recommended quadrature order to compute an integral
@@ -101,7 +114,17 @@ public:
//and d^2 * interfaceBoundary(x) * Kpar
int quadratureOrderInterfaceBoundary () const
{
- return 10;//3;
+ return 8;
+ }
+
+ //returns the recommended quadrature order to compute integrals
+ //over x^2 * (Kpar * grad(d) * grad(d))
+ //and x^2 * d^2
+ //and x^2 * d * (Kpar * grad(d))
+ //and x * d^2 * Kpar
+ int quadratureOrder_Kparallel () const
+ {
+ return 7;
}
private:
diff --git a/dune/mmdg/problems/mmdgproblem5.hh b/dune/mmdg/problems/mmdgproblem5.hh
index ee614898bdfdaadf48f5e3b1e97aa5986101799d..7b276a0e37a5864b802c7a67252cdaccdc590eb9 100644
--- a/dune/mmdg/problems/mmdgproblem5.hh
+++ b/dune/mmdg/problems/mmdgproblem5.hh
@@ -4,7 +4,7 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem
+//a coupled dg problem with constant aperture
template<class Coordinate, class Scalar = double>
class MMDGProblem5 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -29,7 +29,7 @@ public:
return 0.5 / (1.0 + pos[1] * pos[1]);
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -51,7 +51,7 @@ public:
return 0.5 * d_ * d_ - 2.0 * sqrt(2.0) / (1.0 + pos[1] * pos[1]);
}
- //aperture d of the fracture at position pos
+ //constant aperture d of the fracture
Scalar aperture (const Coordinate& pos) const
{
return d_;
diff --git a/dune/mmdg/problems/mmdgproblem6.hh b/dune/mmdg/problems/mmdgproblem6.hh
index ee0320ff268b4b0f3b7e72d58e47194a6bc91bc3..90d267497bc2170f5499ba26a83921a6380b22d4 100644
--- a/dune/mmdg/problems/mmdgproblem6.hh
+++ b/dune/mmdg/problems/mmdgproblem6.hh
@@ -4,7 +4,13 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem
+//a coupled dg problem with non-constant exponentially-shaped aperture
+//with exact solution for arbitrary coupling constant xi > 0.5
+//NOTE: To deal with the non-constant aperture, an adaption of the coupled DG
+// scheme implemented in mmdg.hh is used to solve this problem. Thus, the
+// interface source term and the tangential permeability per aperture
+// inside the fracture are multiplied with the aperture function d in
+// the following problem definition
template<class Coordinate, class Scalar = double>
class MMDGProblem6 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -13,11 +19,11 @@ public:
using Base = CoupledDGProblem<Coordinate, Scalar>;
using Matrix = typename Base::Matrix;
- //constructor
- MMDGProblem6 (const Scalar d0) : d0_(d0), xi_(1.0) {}
+ //constructor with default coupling parameter xi = 1.0
+ MMDGProblem6 (const Scalar dmax) : dmax_(dmax), xi_(1.0) {}
- //constructor
- MMDGProblem6 (const Scalar d0, const Scalar xi) : d0_(d0), xi_(xi) {}
+ //constructor with arbitrary coupling parameter xi
+ MMDGProblem6 (const Scalar dmax, const Scalar xi) : dmax_(dmax), xi_(xi) {}
//the exact bulk solution at position pos
Scalar exactSolution (const Coordinate& pos) const
@@ -25,20 +31,20 @@ public:
Scalar xArg = (pos[0] < 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
: sqrt(2.0) * (2.0 * pos[0] + 1.0);
- return std::cos(4.0 * pos[1]) * std::cosh(xArg);
+ return cos(4.0 * pos[1]) * cosh(xArg);
}
//the exact solution on the interface at position pos
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
- constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
- constexpr Scalar sinhsqrt2 = std::sinh(sqrt(2.0));
+ constexpr Scalar coshsqrt2 = cosh(sqrt(2.0));
+ constexpr Scalar sinhsqrt2 = sinh(sqrt(2.0));
- return std::cos(4.0 * pos[1]) * (coshsqrt2 * coshsqrt2
- - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
+ return cos(4.0 * pos[1]) * ( coshsqrt2 * coshsqrt2
+ - ( 2 * xi_ - 1 ) * sinhsqrt2 * sinhsqrt2 );
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -50,34 +56,36 @@ public:
Scalar xArg = (pos[0] < 0.5) ? sqrt(2.0) * (2.0 * pos[0] - 1.0)
: sqrt(2.0) * (2.0 * pos[0] + 1.0);
- return 8.0 * std::cos(4.0 * pos[1]) * std::cosh(xArg);
+ return 8.0 * cos(4.0 * pos[1]) * cosh(xArg);
}
//interface source term at position pos
+ //NOTE: contains an extra factor d
Scalar qInterface (const Coordinate& pos) const
{
- const Scalar cos4y = std::cos(4.0 * pos[1]);
+ const Scalar cos4y = 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));
+ constexpr Scalar coshsqrt2 = cosh(sqrt(2.0));
+ constexpr Scalar sinhsqrt2 = sinh(sqrt(2.0));
const Scalar hyperbolicFactor =
(coshsqrt2 * coshsqrt2 - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
- return -16.0 * aperture(pos) * d0_ * d0_ * std::exp(-8.0 * pos[1]) * hyperbolicFactor * (
- cos4y + 3.0 * std::sin(4.0 * pos[1]) ) - sqrt2 * std::sinh(sqrt2) * cos4y;
+ return -16.0 * aperture(pos) * dmax_ * dmax_ * exp(-8.0 * pos[1]) *
+ hyperbolicFactor * ( cos4y + 3.0 * sin(4.0 * pos[1]) )
+ - sqrt2 * sinh(sqrt2) * cos4y;
}
//aperture d of the fracture at position pos
Scalar aperture (const Coordinate& pos) const
{
- return d0_ * std::exp(-4.0 * pos[1]);
+ return dmax_ * exp(-4.0 * pos[1]);
}
//tangential gradient of the aperture d of the fracture at position pos
Coordinate gradAperture (const Coordinate& pos) const
{
Coordinate gradD(0.0);
- gradD[1] = -4.0 * d0_ * std::exp(-4.0 * pos[1]);
+ gradD[1] = -4.0 * dmax_ * exp(-4.0 * pos[1]);
return gradD;
}
@@ -85,17 +93,18 @@ public:
//at position pos
Scalar Kperp (const Coordinate& pos) const
{
- return sqrt(2.0) / std::tanh(sqrt(2.0));
+ return sqrt(2.0) / tanh(sqrt(2.0));
}
//tangential permeability tensor of the interface at position pos
+ //NOTE: contains an extra factor d
Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
for (int i = 0; i < dim; i++)
{
- permeability[i][i] = 1.0 * aperture(pos);
+ permeability[i][i] = aperture(pos);
}
return permeability;
@@ -142,7 +151,7 @@ public:
}
private:
- Scalar d0_; //prefactor of the aperture
+ Scalar dmax_; //maximum aperture
Scalar xi_; //coupling constant
};
diff --git a/dune/mmdg/problems/mmdgproblem7.hh b/dune/mmdg/problems/mmdgproblem7.hh
index 19b89847227aafb6952812a6a01f2a9c6ef3a3a3..b0a606327bb2765ea3858f5a0943d9f8886731f6 100644
--- a/dune/mmdg/problems/mmdgproblem7.hh
+++ b/dune/mmdg/problems/mmdgproblem7.hh
@@ -4,7 +4,16 @@
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
-//a coupled dg problem
+//coupled dg problem with non constant cosine-shaped fracture,
+//this is the averaged problem corresponding to the fully resolved problem dg4
+//NOTE: The "exact" solution of the averaged problem that is implemented below
+// does not solve the reduced problem. Thus, there is no convergence for
+// this problem.
+//NOTE: To deal with the non-constant aperture, an adaption of the coupled DG
+// scheme implemented in mmdg.hh is used to solve this problem. Thus, the
+// interface source term and the tangential permeability per aperture
+// inside the fracture are multiplied with the aperture function d in
+// the following problem definition
template<class Coordinate, class Scalar = double>
class MMDGProblem7 : public CoupledDGProblem<Coordinate, Scalar>
{
@@ -32,14 +41,15 @@ public:
}
//the exact solution on the interface at position pos
+ //NOTE: this is the averaged pressure form problem dg4 but not an exact
+ // solution of the system with reduced fracture
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
-
const Scalar d = aperture(pos);
- return ( 1.0 + d ) / d * std::sinh(d);
+ return ( 1.0 + d ) / d * sinh(d);
}
- //indicates whether an exact solution is implemented for the problem
+ //exact solution is implemented
bool hasExactSolution () const
{
return true;
@@ -61,27 +71,27 @@ public:
//pos in Omega_2 (right of the fracture)
return -2.0 * pos[0] * exp(d) / ( 1.0 + d )
* ( dPrimePrime + dPrime2 * d / ( 1.0 + d ) );
-
}
//interface source term at position pos
+ //NOTE: contains an extra factor d
Scalar qInterface (const Coordinate& pos) const
{
- return -( 4.0 + aperturePrimePrime(pos) ) * std::sinh(aperture(pos)) * aperture(pos);
+ return -( 4.0 + aperturePrimePrime(pos) ) * sinh(aperture(pos))
+ * aperture(pos);
}
//aperture d of the fracture at position pos
Scalar aperture (const Coordinate& pos) const
{
- return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
- ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+ return dmin_ + 0.5 * ( dmax_ - dmin_ ) * ( 1.0 + cos(8.0 * M_PI * pos[1]) );
}
//tangential gradient of the aperture d of the fracture at position pos
Coordinate gradAperture (const Coordinate& pos) const
{
Coordinate gradD(0.0);
- gradD[1] = -4.0 * M_PI * ( dmax_ - dmin_ ) * std::sin(8.0 * M_PI * pos[1]);
+ gradD[1] = -4.0 * M_PI * ( dmax_ - dmin_ ) * sin(8.0 * M_PI * pos[1]);
return gradD;
}
@@ -114,6 +124,13 @@ public:
return permeability;
}
+ //returns the recommended quadrature order to compute an integral
+ //over an entry K[i][j] of the permeability tensor
+ int quadratureOrder_K () const
+ {
+ return 10;
+ }
+
//permeability per aperture of the fracture in normal direction
//at position pos
Scalar Kperp (const Coordinate& pos) const
@@ -124,13 +141,14 @@ public:
}
//tangential permeability tensor of the interface at position pos
+ //NOTE: contains an extra factor d
Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
for (int i = 0; i < dim; i++)
{
- permeability[i][i] = 1.0;// / aperture(pos);
+ permeability[i][i] = 1.0;
}
return permeability;
@@ -176,12 +194,18 @@ public:
return 10;
}
+ //returns the recommended quadrature order to compute an integral
+ //over x * Kperp(x)
+ int quadratureOrder_Kperp () const
+ {
+ return 10;
+ }
+
private:
//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]);
+ return -32.0 * M_PI * M_PI * ( dmax_ - dmin_ ) * cos(8.0 * M_PI * pos[1]);
}
const Scalar dmin_; //minimum aperture of the fracture
diff --git a/src/mmdgAnalysis.py b/src/mmdgAnalysis.py
index d403120c698a3b4b61c189ecc732de6c9392dcee..5ca18824b79ef9b8b4f346856e272cbde8946408 100644
--- a/src/mmdgAnalysis.py
+++ b/src/mmdgAnalysis.py
@@ -57,7 +57,7 @@ def readData():
# determines the problem type and the parameter "fastEOC" from the ini file
def getProblemType ():
- # storage for problem type (1 to 6)
+ # storage for problem type (1 to 7)
problemType = None
# boolean that determines whether to exclude computationally expensive grids
diff --git a/src/plotDG.py b/src/plotDG.py
index dbbed59eac9b5503802f0b16e6d93c287bfe3bbb..5bfaf15154960f4c933c3db701a9ee26d8fda828 100644
--- a/src/plotDG.py
+++ b/src/plotDG.py
@@ -3,6 +3,22 @@ import matplotlib.pyplot as plt
from mpltools.annotation import slope_marker
import matplotlib.pylab as pylab
+# determines the problem type from the ini file
+def getProblemType ():
+ # storage for problem type (1 to 5)
+ problemType = None
+
+ # read file content and get problem type
+ with open("parameterDG.ini", "r") as file:
+ lines = file.read().splitlines()
+
+ for i in range(len(lines)):
+ if (lines[i].startswith("problem")):
+ problemType = int(lines[i].rstrip()[-1])
+
+ return problemType
+
+
params = {'legend.fontsize': 'x-large',
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
@@ -10,17 +26,24 @@ params = {'legend.fontsize': 'x-large',
'ytick.labelsize':'x-large'}
pylab.rcParams.update(params)
-data_1d = np.loadtxt("plots/dgAnalysis_1d")
+problemType = getProblemType()
+
+data_1d = None
data_2d = np.loadtxt("plots/dgAnalysis_2d")
data_3d = np.loadtxt("plots/dgAnalysis_3d")
-plt.loglog(np.reciprocal(data_1d[:,2]), data_1d[:,0], 'g^-', label='$n=1$')
+if problemType in [1,2,3]:
+ data_1d = np.loadtxt("plots/dgAnalysis_1d")
+ plt.loglog(np.reciprocal(data_1d[:,2]), data_1d[:,0], 'g^-', label='$n=1$')
+
plt.loglog(np.reciprocal(data_2d[:,2]), data_2d[:,0], 'bo-', label='$n=2$')
plt.loglog(np.reciprocal(data_3d[:,2]), data_3d[:,0], 'rs-', label='$n=3$')
+slopeSize = 0.15 if problemType in [1,2,3] else 0.075
x = 1.0 / data_3d[3,2]
-y = 10 ** (0.2 * np.log10(data_3d[3,0]) + 0.8 * np.log10(data_3d[4,0]))
-slope_marker((x, y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
+y = 10 ** (0.1 * np.log10(data_3d[3,0]) + 0.9 * np.log10(data_3d[4,0]))
+
+slope_marker((x, y), -2, size_frac = slopeSize, invert = True, pad_frac = 0.12,\
poly_kwargs={'ec': 'black', 'fill': False}, \
text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})
plt.xlabel('$h^{-1}$')
diff --git a/src/plotMMDG.py b/src/plotMMDG.py
index 40a6e68899afebaf702e069f05d59e7dd0b0dd9a..f4f67c35698160d55c07560ab47e894cf5cf205f 100644
--- a/src/plotMMDG.py
+++ b/src/plotMMDG.py
@@ -3,6 +3,23 @@ import matplotlib.pyplot as plt
from mpltools.annotation import slope_marker
import matplotlib.pylab as pylab
+
+# determines the problem type from the ini file
+def getProblemType ():
+ # storage for problem type (1 to 7)
+ problemType = None
+
+ # read file content and get problem type
+ with open("parameterMMDG.ini", "r") as file:
+ lines = file.read().splitlines()
+
+ for i in range(len(lines)):
+ if (lines[i].startswith("problem")):
+ problemType = int(lines[i].rstrip()[-1])
+
+ return problemType
+
+
params = {'legend.fontsize': 'x-large',
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
@@ -13,6 +30,8 @@ pylab.rcParams.update(params)
data_2d = np.loadtxt("plots/mmdgAnalysis_2d")
data_3d = np.loadtxt("plots/mmdgAnalysis_3d")
+problemType = getProblemType()
+
# plot total error
plt.figure()
plt.loglog(np.reciprocal(np.maximum(data_2d[:,3], data_2d[:,4])), data_2d[:,2],\
@@ -20,7 +39,7 @@ plt.loglog(np.reciprocal(np.maximum(data_2d[:,3], data_2d[:,4])), data_2d[:,2],\
plt.loglog(np.reciprocal(np.maximum(data_3d[:,3], data_3d[:,4])), data_3d[:,2],\
'rs-', label='$n=3$')
-if data_3d[0,2] > data_3d[-1,2]:
+if data_3d[0,2] > data_3d[-1,2] and problemType not in [1,7]:
x = 1.0 / np.maximum(data_3d[3,3], data_3d[3,4])
y = 10 ** (0.2 * np.log10(data_3d[3,2]) + 0.8 * np.log10(data_3d[4,2]))
slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
@@ -38,7 +57,7 @@ plt.figure()
plt.loglog(np.reciprocal(data_2d[:,3]), data_2d[:,0], 'bo-', label='$n=2$')
plt.loglog(np.reciprocal(data_3d[:,3]), data_3d[:,0], 'rs-', label='$n=3$')
-if data_3d[0,0] > data_3d[-1,0]:
+if data_3d[0,0] > data_3d[-1,0] and problemType not in [1,7]:
x = 1.0 / data_3d[3,3]
y = 10 ** (0.2 * np.log10(data_3d[3,0]) + 0.8 * np.log10(data_3d[4,0]))
slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
@@ -55,9 +74,9 @@ plt.figure()
plt.loglog(np.reciprocal(data_2d[:,4]), data_2d[:,1], 'bo-', label='$n=2$')
plt.loglog(np.reciprocal(data_3d[:,4]), data_3d[:,1], 'rs-', label='$n=3$')
-if data_3d[0,1] > data_3d[-1,1]:
+if data_3d[0,1] > data_3d[-1,1] and problemType not in [1,7]:
x = 1.0 / data_3d[3,4]
- y = 10 ** (0.1 * np.log10(data_3d[3,1]) + 0.9 * np.log10(data_3d[4,1]))
+ y = 10 ** (0.05 * np.log10(data_3d[3,1]) + 0.95 * np.log10(data_3d[4,1]))
slope_marker((x,y), -2, size_frac = 0.15, invert = True, pad_frac = 0.12, \
poly_kwargs={'ec': 'black', 'fill': False}, \
text_kwargs={'size': 'x-large', 'text': '2', 'usetex': True})