Skip to content
Snippets Groups Projects
Commit 63a01256 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

add problem mmdg7

parent 45af3c7b
No related branches found
No related tags found
No related merge requests found
#ifndef __DUNE_MMDG_MMDGPROBLEM7_HH__
#define __DUNE_MMDG_MMDGPROBLEM7_HH__
#include <cmath>
#include <dune/mmdg/problems/coupleddgproblem.hh>
//a coupled dg problem
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;
static constexpr Scalar dmin = 0.001; //minimum aperture of the fracture
static constexpr Scalar dmax = 0.002; //maximum aperture of the fracture
//constructor
MMDGProblem7 (const Scalar d0) : d0_(d0), xi_(1.0) {}
//constructor
MMDGProblem7 (const Scalar d0, const Scalar xi) : d0_(d0), xi_(xi) {}
//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
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
const Scalar d = aperture(pos);
return ( 1.0 + d ) / d * std::sinh(d);
}
//indicates whether an exact solution is implemented for the problem
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
Scalar qInterface (const Coordinate& pos) const
{
return -( 4.0 + aperturePrimePrime(pos) ) * std::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]) );
}
//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]);
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;
}
//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
Matrix Kparallel (const Coordinate& pos) const
{
Matrix permeability(0.0);
for (int i = 0; i < dim; i++)
{
permeability[i][i] = 1.0;// / aperture(pos);
}
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;
}
private:
Scalar d0_; //prefactor of the aperture
Scalar xi_; //coupling constant
//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]);
}
};
#endif
......@@ -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)
{
......@@ -101,6 +102,12 @@ int main(int argc, char** argv)
problem = new MMDGProblem6<Coordinate>(aperture, xi);
gridType = "mmdg2_C";
}
else if (pt["problem"] == "mmdg7")
{
xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
problem = new MMDGProblem7<Coordinate>(aperture, xi);
gridType = "mmdg2_C";
}
else
{
DUNE_THROW(Dune::Exception,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment