diff --git a/README.md b/README.md
index b4736659764719d1674dcef972a01035e3f04046..af64cb26b9cd61d72db8781703d6ea3bea7c8ba4 100644
--- a/README.md
+++ b/README.md
@@ -51,8 +51,8 @@ You will find the plots `mmdg_bulk.pdf`, `mmdg_interface.pdf` and `mmdg_total.pd
Test Problems
-------------
-The `dune-mmdg` module includes three test problems, `dg1` to `dg3`, for the bulk scheme and six test problems, `mmdg1` to `mmdg6`, for the coupled scheme.
-The test problems `dg2`, `mmdg2`, `mmdg5` and `mmdg6` correspond to the examples 5.1 to 5.4 in [1].
+The `dune-mmdg` module includes five test problems, `dg1` to `dg5`, for the bulk scheme and seven test problems, `mmdg1` to `mmdg7`, for the coupled scheme.
+The test problems `dg2`, `mmdg2`, `mmdg5`, `mmdg6` and `dg4`/`mmdg7` correspond to the examples 5.1, 5.2, 5.3, 5.4/5.5 and 5.6 in [1] in the given order.
Parameters
@@ -61,16 +61,21 @@ Parameters
For the bulk scheme the following parameters can be set in the file `parameterDG.ini` in the `src` directory:
- `mu0`: the prefactor of the penalty parameter (must be sufficiently large)
-- `problem`: the name of a test problem (`dg1` to `dg3`)
+- `problem`: the name of a test problem (`dg1` to `dg5`)
+- `dmax`: the maximum value of the aperture function (only relevant for the problems `dg4` and `dg5` that include a fracture)
+- `dmin`: the minimum value of the aperture function (only relevant for the problems `dg4` and `dg5` that include a fracture)
+- `solver`: (optional) the name of the solver that is applied to the system of linear equations (`Cholmod` or `UMFPack`)
- `fastEOC`: (optional) exclude computationally expensive grids from the EOC loop
For the coupled scheme the following parameters can be set in the file `parameterMMDG.ini` in the `src` directory:
- `mu0`: the prefactor of the penalty parameter (must be sufficiently large)
-- `problem`: the name of a test problem (`mmdg1` to `mmdg6`)
-- `d`: the aperture of the fracture or the prefactor of the aperture function
+- `problem`: the name of a test problem (`mmdg1` to `mmdg7`)
+- `dmax`: the maximum value of the aperture of the fracture or respective constant aperture value (mandatory for all problems)
+- `dmin`: the minimum value of the aperture of the fracture (only relevant for problem `mmdg7`)
- `xi`: (optional) the value of the coupling parameter (must be larger than 0.5)
- `gridfile`: (optional) the prefix of a grid file name in the directory `src/grids`
+- `solver`: (optional) the name of the solver that is applied to the system of linear equations (`Cholmod` or `UMFPack`)
- `fastEOC`: (optional) exclude computationally expensive grids from the EOC loop
diff --git a/dune/mmdg/clockhelper.hh b/dune/mmdg/clockhelper.hh
index 05c3ddc2afd28bba0deb357fc5f9edbcc0ca4790..850e1e49d0e82ca344ceea4ee4a780befd89f38c 100644
--- a/dune/mmdg/clockhelper.hh
+++ b/dune/mmdg/clockhelper.hh
@@ -3,6 +3,7 @@
#include <chrono>
+//contains a function for run time measurements
class ClockHelper
{
public:
@@ -12,8 +13,8 @@ class ClockHelper
//precision of a microsecond
static double elapsedRuntime(const Clock::time_point begin)
{
- return 1e-6*std::chrono::duration_cast<std::chrono::microseconds>
- (Clock::now() - begin).count();;
+ return 1e-6 * std::chrono::duration_cast<std::chrono::microseconds>
+ ( Clock::now() - begin ).count();;
}
};
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 623c82561d0d696e62e2478319a2623e94796fd5..8a00b9bc87c1205c1abb56a5b8e237ad5f8ef674 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -14,9 +14,11 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/cholmod.hh>
+#include <dune/istl/umfpack.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
#include <dune/mmdg/eigenvaluehelper.hh>
+#include <dune/mmdg/solvertype.hh>
template<class GridView, class Mapper, class Problem>
class DG
@@ -40,12 +42,13 @@ public:
using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
//constructor
- DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) :
+ DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
+ const SolverType& solver) :
gridView_(gridView), mapper_(mapper), problem_(problem),
- dof_((1 + dim) * gridView.size(0))
+ dof_((1 + dim) * gridView.size(0)), solver_(solver)
{
//initialize stiffness matrix A, load vector b and solution vector d
- A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit);
b = Vector(dof_);
d = Vector(dof_);
}
@@ -56,16 +59,31 @@ public:
assembleSLE(mu0);
A->compress();
+ //solve system of linear equations
Dune::InverseOperatorResult result;
- Dune::Cholmod<Vector> solver;
- solver.setMatrix(*A);
- solver.apply(d, b, result);
+
+ if (solver_ == SolverType::Cholmod)
+ {
+ Dune::Cholmod<Vector> solver;
+ solver.setMatrix(*A);
+ solver.apply(d, b, result);
+ }
+ else if (solver_ == SolverType::UMFPack)
+ {
+ Dune::UMFPack<Matrix> solver(*A);
+ solver.apply(d, b, result);
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception, "ERROR: Specified solver type not "
+ "implemented");
+ }
//write solution to a vtk file
writeVTKoutput();
}
- //compute L2 error using a quadrature rule
+ //computes the L2 error using a quadrature rule
Scalar computeL2error(const int quadratureOrder = 5)
{
if (!problem_.hasExactSolution())
@@ -112,12 +130,14 @@ public:
const int dof_; //degrees of freedom
protected:
+ //constructor with custom value of dof_
DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
- const int dof) :
- gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
+ const SolverType solver, const int dof) :
+ gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof),
+ solver_(solver)
{
//initialize stiffness matrix A, load vector b and solution vector d
- A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
+ A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit);
b = Vector(dof_);
d = Vector(dof_);
}
@@ -151,7 +171,7 @@ protected:
//in the system of linear equations (SLE) Ad = b,
//the index elemIdxSLE refers to the basis function phi_elem,0
//and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
- //basis function phi_elem,i
+ //basis functions phi_elem,i
const int elemIdxSLE = (dim + 1)*elemIdx;
//lambda to fill load vector b using a quadrature rule
@@ -212,7 +232,7 @@ protected:
//create storage for multiply used integral values
//linearIntegrals[i] = int_intersct x[i] ds
- //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds //NOTE: is there a better type for a symmetric matrix?
+ //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds
//KIntegrals[i] = int_intersct K[:,i]*normal ds
//KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds
Coordinate linearIntegrals(0.0);
@@ -267,7 +287,7 @@ protected:
}
};
- //apply quadrature rules in intersct
+ //apply quadrature rules on intersct
applyQuadratureRule(intersctGeo, rule_2ndOrder, lambda_2ndOrder);
applyQuadratureRule(intersctGeo, rule_K_Edge, lambda_K_Edge);
applyQuadratureRule(intersctGeo, rule_Kplus1_Edge, lambda_Kplus1_Edge);
@@ -561,6 +581,7 @@ protected:
std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
+ SolverType solver_; //solver type (Cholmod or UMFPack)
private:
//writes the solution to a vtk file
diff --git a/dune/mmdg/eigenvaluehelper.hh b/dune/mmdg/eigenvaluehelper.hh
index a3e9cf0110cf2db8cb68bd30bea57c92d60f4a14..7fe56886e13b5457466458d04b91826281c81655 100644
--- a/dune/mmdg/eigenvaluehelper.hh
+++ b/dune/mmdg/eigenvaluehelper.hh
@@ -10,7 +10,7 @@
/** \file
* \brief Eigenvalue computations for the Dune::FieldMatrix class
- * NOTE: taken from DUNE-common, fmatrixev.hh, with adapted tolerance
+ * NOTE: taken from DUNE-common (fmatrixev.hh) with minor changes
*/
class EigenvalueHelper {
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index a83541dbf9d2d190a9b2c948f63fc04c0adfa8df..0889dafd7bd1519bd9a1b3c4de3be1d53bebee0b 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -30,20 +30,36 @@ public:
//constructor
MMDG (const GridView& gridView, const Mapper& mapper,
const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
- const Problem& problem) :
- Base(gridView, mapper, problem,
+ const Problem& problem, const SolverType& solver) :
+ Base(gridView, mapper, problem, solver,
(1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
iGridView_(iGridView), iMapper_(iMapper) {}
const void operator() (const Scalar mu0, const Scalar xi)
{
+ //assemble stiffness matrix A and load vector b
assembleSLE(mu0, xi);
Base::A->compress();
+ //solve system of linear equations
Dune::InverseOperatorResult result;
- Dune::Cholmod<Vector> solver;
- solver.setMatrix(*Base::A);
- solver.apply(Base::d, Base::b, result);
+
+ if (Base::solver_ == SolverType::Cholmod)
+ {
+ Dune::Cholmod<Vector> solver;
+ solver.setMatrix(*Base::A);
+ solver.apply(Base::d, Base::b, result);
+ }
+ else if (Base::solver_ == SolverType::UMFPack)
+ {
+ Dune::UMFPack<Matrix> solver(*Base::A);
+ solver.apply(Base::d, Base::b, result);
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception, "ERROR: Specified solver type not "
+ "implemented");
+ }
//write solution to a vtk file
writeVTKoutput();
@@ -51,7 +67,7 @@ public:
//compute L2 errors using quadrature rules,
//returns {error(bulk), error(interface), error(total)},
- //the total error is calculated using the norm
+ //the total error is calculated with respect to the norm
// norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
const std::array<Scalar, 3> computeL2error(const int quadratureOrder = 5)
{
@@ -90,8 +106,8 @@ public:
Base::applyQuadratureRule(iGeo, qrule, qlambda);
}
- return {bulkError, sqrt(interfaceError),
- sqrt(bulkError * bulkError + interfaceError)};
+ return { bulkError, sqrt(interfaceError),
+ sqrt(bulkError * bulkError + interfaceError) };
}
const InterfaceGridView& iGridView_;
@@ -147,6 +163,7 @@ private:
//and load vector b with bulk entries
Base::assembleSLE(mu0);
+ //fill in coupling and interface entries
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdx = iMapper_.index(iElem);
@@ -878,6 +895,8 @@ private:
return iBasis;
}
+ //write exakt and numerical bulk and interface solution as well as the
+ //aperture function of the problem to vtk files
void writeVTKoutput () const
{
//degrees of freedom of bulk and interface grids
@@ -894,6 +913,7 @@ private:
//create the output using a linear interpolation for each interface element
Dune::DynamicVector<Scalar> iPressure(interfaceDOF, 0.0);
Dune::DynamicVector<Scalar> exactIPressure(interfaceDOF, 0.0);
+ Dune::DynamicVector<Scalar> aperture(interfaceDOF, 0.0);
for (const auto& iElem : elements(iGridView_))
{
@@ -914,6 +934,10 @@ private:
Base::problem_.exactInterfaceSolution(iGeo.corner(k));
}
+ //aperture evaluated at the kth corner of iElem
+ aperture[iElemIdxOutput + k] =
+ Base::problem_.aperture(iGeo.corner(k));
+
//contribution of the basis function
// phi_iElem,0 (x) = indicator(iElem);
//at the kth corner of iElem
@@ -934,6 +958,8 @@ private:
vtkWriter(iGridView_, Dune::VTK::nonconforming);
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(iGridView_, iPressure, "interfacePressure")));
+ vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
+ new P1Function(iGridView_, aperture, "aperture")));
if (Base::problem_.hasExactSolution())
{
diff --git a/dune/mmdg/nonconformingp1vtkfunction.hh b/dune/mmdg/nonconformingp1vtkfunction.hh
index 9ae9ec80968a149730263320846d4e4715cbfb0e..45a69ad7582d8e8af9245c9ff7c33063ad3ba4f8 100644
--- a/dune/mmdg/nonconformingp1vtkfunction.hh
+++ b/dune/mmdg/nonconformingp1vtkfunction.hh
@@ -1,8 +1,8 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_MMFD_FUNCTION_HH
-#define DUNE_MMFD_FUNCTION_HH
+#ifndef DUNE_MMDG_FUNCTION_HH
+#define DUNE_MMDG_FUNCTION_HH
#include <string>
@@ -17,6 +17,7 @@
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/io/file/vtk/function.hh>
+//NOTE: without changes from DUNE-MMFD by Samuel Burbulla
namespace Dune
{
//////////////////////////////////////////////////////////////////////
@@ -143,4 +144,4 @@ namespace Dune
} // namespace Dune
-#endif // DUNE_MMFD_FUNCTION_HH
+#endif // DUNE_MMDG_FUNCTION_HH
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index 9f7774ff6e3611f92dfcf8b16d21727815fda430..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);
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
new file mode 100644
index 0000000000000000000000000000000000000000..07b8dd0aea89b1e81d52ece3c0d719337f8bb365
--- /dev/null
+++ b/dune/mmdg/problems/dgproblem4.hh
@@ -0,0 +1,164 @@
+#ifndef __DUNE_MMDG_DGPROBLEM4_HH__
+#define __DUNE_MMDG_DGPROBLEM4_HH__
+
+#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>
+{
+ public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = DGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ //constructor, require dim != 1
+ DGProblem4 (const Scalar dmin, const Scalar dmax) : dmin_(dmin), dmax_(dmax)
+ {
+ if (dim == 1)
+ {
+ 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)
+ int quadratureOrder_q () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * boundary(x) and K(x) * boundary(x)
+ int quadratureOrderBoundary () const
+ {
+ return 10;
+ }
+
+ //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;
+ }
+
+ 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/dune/mmdg/problems/dgproblem5.hh b/dune/mmdg/problems/dgproblem5.hh
new file mode 100644
index 0000000000000000000000000000000000000000..678c0e5283d87b8e6db6558ffc6362a31c65a397
--- /dev/null
+++ b/dune/mmdg/problems/dgproblem5.hh
@@ -0,0 +1,131 @@
+#ifndef __DUNE_MMDG_DGPROBLEM5_HH__
+#define __DUNE_MMDG_DGPROBLEM5_HH__
+
+#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>
+{
+ public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = DGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ static constexpr Scalar Kbulk = 1.0; //bulk permeability
+ static constexpr Scalar K_par = 3.0; //tangential fracture permeability
+ static constexpr Scalar K_perp = 0.5; //normal fracture permeability
+
+ //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
+ {
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+ { //pos in Omega_1 (left of the fracture)
+ return -0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+ }
+
+ if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+ { //pos in Omega_2 (right of the fracture)
+ return 0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) * pos[1];
+ }
+
+ //exact solution is implemented
+ bool hasExactSolution() const
+ {
+ return true;
+ };
+
+ //source term at position pos
+ Scalar q (const Coordinate& pos) const
+ {
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+ { //pos in Omega_1 (left of the fracture)
+ return -4.0 * M_PI * Kbulk * ( dmax_ - dmin_ ) *
+ ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+ + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+ { //pos in Omega_2 (right of the fracture)
+ return 4.0 * M_PI * Kbulk * ( dmax_ - dmin_ ) *
+ ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+ + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return -2.0 * K_perp * pos[1];
+ }
+
+ //permeability tensor at position pos
+ Matrix K (const Coordinate& pos) const
+ {
+ Matrix permeability(0.0);
+
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) ||
+ pos[0] > 0.5 * (1.0 + aperture(pos)))
+ { //pos in Omega_1 or Omega_2 (outside the fracture)
+ for (int i = 0; i < dim; i++)
+ {
+ permeability[i][i] = Kbulk;
+ }
+ }
+ else
+ { //pos in Omega_R (inside the fracture)
+ permeability[0][0] = K_perp;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = K_par;
+ }
+ }
+
+ return permeability;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * q(x)
+ int quadratureOrder_q () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * boundary(x) and K(x) * boundary(x)
+ int quadratureOrderBoundary () const
+ {
+ return 10;
+ }
+
+ private:
+ //aperture of the fracture
+ const Scalar aperture (const Coordinate& pos) const
+ {
+ return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
+ ( 1.0 + 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/dune/mmdg/problems/mmdgproblem1.hh b/dune/mmdg/problems/mmdgproblem1.hh
index f35a61ff999ac3eb073e87408dc67cf92c07a1b1..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;
@@ -37,7 +42,40 @@ public:
//aperture d of the fracture at position pos
Scalar aperture (const Coordinate& pos) const
{
- return d_;
+ return d_ * (pos[1] + 0.01);
+ }
+
+ //interface source term at position pos
+ //NOTE: contains an extra factor d
+ Scalar qInterface (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+ const Coordinate gradD = gradAperture(pos);
+
+ // 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;
+ }
+
+ //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] = aperture(pos);
+ }
+
+ return permeability;
}
//returns the recommended quadrature order to compute an integral
@@ -47,8 +85,34 @@ public:
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 cf09ef6b727e7b946854d364c29605e582e9b6ba..c8adbfff494997ba68a2cad962f7da6bb05dcb07 100644
--- a/dune/mmdg/problems/mmdgproblem4.hh
+++ b/dune/mmdg/problems/mmdgproblem4.hh
@@ -4,17 +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; //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
@@ -23,7 +30,7 @@ public:
if (pos[0] > 0.5)
{
- return solution += d_;
+ return solution += jump;
}
return solution;
@@ -32,41 +39,56 @@ public:
//the exact solution on the interface at position pos
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
- return 0.25 - pos[1] * pos[1] + 0.5 * d_;
+ 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_;
+ const Scalar d = aperture(pos);
+ const Scalar dPrime = sqrt(gradAperture(pos) * gradAperture(pos));
+ const Scalar dPrimePrime = 2.0 * d_;
+
+ 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_;
+ 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] = 2.0 * d_ * pos[1];
+ return gradD;
}
//permeability of the fracture in normal direction at position pos
Scalar Kperp (const Coordinate& pos) const
{
- return Scalar(1.0) / d_;
+ 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 / d_;
+ permeability[i][i] = aperture(pos);
}
return permeability;
@@ -79,13 +101,30 @@ public:
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
//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;
+ 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 311042544d461531f0bdf64b144b02d35dfcfa04..90d267497bc2170f5499ba26a83921a6380b22d4 100644
--- a/dune/mmdg/problems/mmdgproblem6.hh
+++ b/dune/mmdg/problems/mmdgproblem6.hh
@@ -4,18 +4,26 @@
#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>
{
public:
static constexpr int dim = Coordinate::dimension;
+ 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
@@ -23,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;
@@ -48,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 * 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;
}
@@ -83,7 +93,21 @@ 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] = aperture(pos);
+ }
+
+ return permeability;
}
//returns the recommended quadrature order to compute an integral
@@ -127,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
new file mode 100644
index 0000000000000000000000000000000000000000..b0a606327bb2765ea3858f5a0943d9f8886731f6
--- /dev/null
+++ b/dune/mmdg/problems/mmdgproblem7.hh
@@ -0,0 +1,216 @@
+#ifndef __DUNE_MMDG_MMDGPROBLEM7_HH__
+#define __DUNE_MMDG_MMDGPROBLEM7_HH__
+
+#include <cmath>
+#include <dune/mmdg/problems/coupleddgproblem.hh>
+
+//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>
+{
+public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = CoupledDGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ //constructor
+ MMDGProblem7 (const Scalar dmin, const Scalar dmax)
+ : dmin_(dmin), dmax_(dmax) {}
+
+ //the exact bulk solution at position pos
+ Scalar exactSolution (const Coordinate& pos) const
+ {
+ const Scalar d = aperture(pos);
+
+ if ( pos[0] < 0.5 )
+ { //pos in Omega_1 (left of the fracture)
+ return 2.0 * ( pos[0] + d ) * exp(-d);
+ }
+
+ //pos in Omega_2 (right of the fracture)
+ return 2.0 * pos[0] * exp(d);
+ }
+
+ //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 * sinh(d);
+ }
+
+ //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);
+ const Scalar dPrime2 = gradAperture(pos) * gradAperture(pos);
+ const Scalar dPrimePrime = aperturePrimePrime(pos);
+
+ if ( pos[0] < 0.5 )
+ { //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 );
+ }
+
+ //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) ) * 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 + 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_ ) * sin(8.0 * M_PI * pos[1]);
+ return gradD;
+ }
+
+ //bulk 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 )
+ { //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
+ { //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 );
+ }
+ }
+
+ 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
+ {
+ const Scalar d = aperture(pos);
+
+ return 1.0 / ( d * ( 1.0 + d ) );
+ }
+
+ //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;
+ }
+
+ return permeability;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * boundary(x)
+ int quadratureOrderBoundary () const
+ {
+ return 10;
+ }
+
+ //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 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * qInterface(x)
+ int quadratureOrder_qInterface () const
+ {
+ return 10;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over x * q(x)
+ int quadratureOrder_q () const
+ {
+ return 10;
+ }
+
+ //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 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_ ) * 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/dune/mmdg/solvertype.hh b/dune/mmdg/solvertype.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6ca83002e71a6f558eeab4ad7036f14ab47a2eb1
--- /dev/null
+++ b/dune/mmdg/solvertype.hh
@@ -0,0 +1,6 @@
+//determines the type of solver that is applied to a system of linear equations
+enum SolverType
+{
+ Cholmod,
+ UMFPack
+};
diff --git a/src/dg.cc b/src/dg.cc
index 2dedbca292ccc37ac5a9b0d35be06b423aaf1221..d78f919e8e8cfd963f8f135f5da4efe5079b38fd 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -12,9 +12,8 @@
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
-#include <dune/grid/io/file/dgfparser/dgfug.hh>
-#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/io/file/dgfparser/dgfoned.hh>
#include <dune/mmesh/mmesh.hh>
@@ -24,6 +23,8 @@
#include <dune/mmdg/problems/dgproblem1.hh>
#include <dune/mmdg/problems/dgproblem2.hh>
#include <dune/mmdg/problems/dgproblem3.hh>
+#include <dune/mmdg/problems/dgproblem4.hh>
+#include <dune/mmdg/problems/dgproblem5.hh>
int main(int argc, char** argv)
{
@@ -36,8 +37,8 @@ int main(int argc, char** argv)
static constexpr int dim = GRIDDIM;
- //use YaspGrid if dim = 1 and MMesh otherwise
- using Grid = std::conditional<dim == 1, Dune::YaspGrid<dim>,
+ //use OneDGrid if dim = 1 and MMesh otherwise
+ using Grid = std::conditional<dim == 1, Dune::OneDGrid,
Dune::MovingMesh<dim>>::type;
using GridFactory = Dune::DGFGridFactory<Grid>;
using GridView = typename Grid::LeafGridView;
@@ -52,12 +53,21 @@ int main(int argc, char** argv)
Dune::ParameterTree pt;
Dune::ParameterTreeParser::readINITree("parameterDG.ini", pt);
- //assume constant penalty parameter
- const double mu0 = std::stod(pt["mu0"]);
+ double mu0; //penalty parameter (prefactor)
+
+ try
+ {
+ mu0 = std::stod(pt["mu0"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid penalty parameter or parameter 'mu0' not specified in the "
+ "file parameterDG.ini");
+ }
//create a grid from .dgf file
GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
-
const Grid& grid = *gridFactory.grid();
const GridView& gridView = grid.leafGridView();
@@ -87,13 +97,63 @@ int main(int argc, char** argv)
}
else
{
- DUNE_THROW(Dune::Exception,
- "Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'dg1' to 'dg3')");
+ double dmin; //minimum aperture
+ double dmax; //maximum aperture
+
+ try
+ {
+ dmin = std::stod(pt["dmin"]);
+ dmax = std::stod(pt["dmax"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "The problems dg4 and dg5 require the parameters 'dmin' and 'dmax'. "
+ "Check the file parameterDG.ini and make sure that they are "
+ "specified correctly.");
+ }
+
+ if (pt["problem"] == "dg4")
+ {
+ problem = new DGProblem4<Coordinate>(dmin, dmax);
+ }
+ else if (pt["problem"] == "dg5")
+ {
+ problem = new DGProblem5<Coordinate>(dmin, dmax);
+
+ std::cout << "NOTE: You will not see convergence for problem dg5 as it"
+ " only fulfills continuity of the normal component of the Darcy "
+ "velocity 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. \n\n";
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid problem type or parameter 'problem' not specified in file "
+ "parameterDG.ini (use 'dg1' to 'dg5')");
+ }
}
- DG dg(gridView, mapper, *problem);
+ //determine solver type
+ SolverType solver = SolverType::Cholmod;
+
+ if (pt.hasKey("solver"))
+ {
+ if (pt["solver"] == "UMFPack")
+ {
+ solver = SolverType::UMFPack;
+ }
+ else if (pt["solver"] != "Cholmod")
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid solver type (use 'Cholmod' or 'UMFPack' for the parameter "
+ "'solver' in the file parameterDG.ini)");
+ }
+ }
+ //apply DG scheme
+ DG dg(gridView, mapper, *problem, solver);
dg(mu0);
//error and total run time
@@ -107,7 +167,7 @@ int main(int argc, char** argv)
largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
}
- //write error, runtime and maximum edge length
+ //write error, runtime and maximum edge length to file
std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
if (errorFile.is_open())
diff --git a/src/dgAnalysis.py b/src/dgAnalysis.py
index 4a6edd2bc48ce20c9aca62c095e5137798e42c4b..a54c24d9871000717ffbf7360a502b8b15b829d1 100644
--- a/src/dgAnalysis.py
+++ b/src/dgAnalysis.py
@@ -47,22 +47,26 @@ def readData():
return np.array([error, timeTotal, numberOfElements])
-# determines the parameter "fastEOC" from the ini file
-def determineFastEOC ():
+# determines the problem type and the parameter "fastEOC" from the ini file
+def getProblemType ():
+ # storage for problem type (1 to 5)
+ problemType = None
+
# boolean that determines whether to exclude computationally expensive grids
fastEOC = False
- # read ini file content
+ # 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("fastEOC") and \
+ if (lines[i].startswith("problem")):
+ problemType = int(lines[i].rstrip()[-1])
+ elif (lines[i].startswith("fastEOC") and \
lines[i].rstrip().lower().endswith("true")):
fastEOC = True
- break
- return fastEOC
+ return problemType, fastEOC
# === main ===
@@ -78,13 +82,21 @@ numberOfMeans = 1
numElemsPerDir_12 = np.unique(np.logspace(0, 2.7, 15).astype(int))
numElemsPerDir_3 = np.unique(np.logspace(0, 1.5, 8).astype(int))
-if determineFastEOC(): # exclude computationally expensive grids
+# get problem type (integer from 1 to 5) from ini file
+problemType, fastEOC = getProblemType()
+
+if fastEOC: # exclude computationally expensive grids
numElemsPerDir_12 = numElemsPerDir_12[:-3]
numElemsPerDir_3 = numElemsPerDir_3[:-2]
else:
print("fastEOC = false. This may take some minutes.\n")
-for dim in [1, 2, 3]:
+dimensions = [1, 2, 3]
+
+if problemType in [4, 5]:
+ dimensions = [2, 3]
+
+for dim in dimensions:
numElemsPerDir = numElemsPerDir_3 if (dim == 3) else numElemsPerDir_12
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 7a2bca12179fef4e816757b7f015b7776fc5d5c0..6897c0d2304293067cb869aaa76407f06a54450a 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -26,6 +26,7 @@
#include <dune/mmdg/problems/mmdgproblem4.hh>
#include <dune/mmdg/problems/mmdgproblem5.hh>
#include <dune/mmdg/problems/mmdgproblem6.hh>
+#include <dune/mmdg/problems/mmdgproblem7.hh>
int main(int argc, char** argv)
{
@@ -56,56 +57,97 @@ int main(int argc, char** argv)
Dune::ParameterTree pt;
Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
- const double mu0 = std::stod(pt["mu0"]); //penalty parameter (prefactor)
- const double aperture = std::stod(pt["d"]); //aperture (prefactor)
+ double mu0; //penalty parameter (prefactor)
+ double dmax; //maximum aperture
double xi; //coupling parameter
+ try
+ {
+ mu0 = std::stod(pt["mu0"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid penalty parameter or parameter 'mu0' not specified in the "
+ "file parameterDG.ini");
+ }
+
+ try
+ {
+ dmax = std::stod(pt["dmax"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid maximum aperture or parameter 'dmax' not specified in the "
+ "file parameterDG.ini");
+ }
+
Problem* problem; //problem to be solved
- std::string gridType;
+ std::string gridType; //prefix of the grid file name
//determine problem type
if (pt["problem"] == "mmdg1")
{
- problem = new MMDGProblem1<Coordinate>(aperture);
+ problem = new MMDGProblem1<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg2")
{
- problem = new MMDGProblem2<Coordinate>(aperture);
+ problem = new MMDGProblem2<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg3")
{
- problem = new MMDGProblem3<Coordinate>(aperture);
+ problem = new MMDGProblem3<Coordinate>(dmax);
gridType = "mmdg3_C";
xi = 1.0;
}
else if (pt["problem"] == "mmdg4")
{
- problem = new MMDGProblem4<Coordinate>(aperture);
+ problem = new MMDGProblem4<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg5")
{
- problem = new MMDGProblem5<Coordinate>(aperture);
+ problem = new MMDGProblem5<Coordinate>(dmax);
gridType = "mmdg5_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg6")
{
xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
- problem = new MMDGProblem6<Coordinate>(aperture, xi);
+ problem = new MMDGProblem6<Coordinate>(dmax, xi);
+ gridType = "mmdg2_C";
+ }
+ else if (pt["problem"] == "mmdg7")
+ {
+ double dmin; //minimum aperture
+
+ try
+ {
+ dmin = std::stod(pt["dmin"]);
+ }
+ catch (std::exception& e)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg4 requires the parameter 'dmin'. Check the "
+ "file parameterMMDG.ini and make sure that it is specified.");
+ }
+
+ problem = new MMDGProblem7<Coordinate>(dmin, dmax);
gridType = "mmdg2_C";
+ xi = 0.75;
}
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'mmdg1' to 'mmdg6')");
+ "parameterDG.ini (use 'mmdg1' to 'mmdg7')");
}
if (pt.hasKey("gridfile")) //non-default grid type
@@ -118,7 +160,24 @@ int main(int argc, char** argv)
xi = std::stod(pt["xi"]);
}
- //create a grid from .dgf file
+ //determine solver type
+ SolverType solver = SolverType::Cholmod;
+
+ if (pt.hasKey("solver"))
+ {
+ if (pt["solver"] == "UMFPack")
+ {
+ solver = SolverType::UMFPack;
+ }
+ else if (pt["solver"] != "Cholmod")
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid solver type (use 'Cholmod' or 'UMFPack' for the parameter "
+ "'solver' in the file parameterDG.ini)");
+ }
+ }
+
+ //create the grid from the .msh file
GridFactory gridFactory( "grids/" + gridType + "_" + std::to_string(dim)
+ "d.msh" );
const Grid& grid = *gridFactory.grid();
@@ -130,7 +189,8 @@ int main(int argc, char** argv)
const Mapper mapper(gridView, Dune::mcmgElementLayout());
const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
- MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
+ //apply the coupled DG scheme
+ MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem, solver);
mmdg(mu0, xi);
//errors and total run time
diff --git a/src/mmdgAnalysis.py b/src/mmdgAnalysis.py
index 4aecbc6a4113f44c9c0a66fba986f37334ebcd90..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
@@ -113,7 +113,7 @@ gridfiles_3d = ["A", "B", "C", "D", "E", "F", "G"]
# get problem type (integer from 1 to 6) from ini file
problemType, fastEOC = getProblemType()
-if problemType in [1, 4, 6]: # these problems use the same grids
+if problemType in [1, 4, 6, 7]: # these problems use the same grids
problemType = 2
if fastEOC: # exclude computationally expensive grids
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 275babcaa75b797c5558068c7a394eb716598827..c7ee27f47aab88515617404f860f61d66f9e2618 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,3 +1,6 @@
mu0 = 1000
-problem = dg2
+problem = dg4
+dmax = 1e-1
+dmin = 1e-2
+solver = Cholmod
fastEOC = true
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
index 257871d69143b4d3315df49397f181c785ad5fab..f572d82a712d1c711ad7a6b727344ab9b9723f50 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -1,6 +1,8 @@
mu0 = 1000
-problem = mmdg2
-d = 1e-2
+problem = mmdg6
+dmax = 1e-3
+dmin = 1e-2
#xi = 0.75
gridfile = mmdg2_C
+solver = Cholmod
fastEOC = true
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})