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

add dmin, dmax as member variables; adapt analysis scripts to new problems

parent 37e6bf9a
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/mmdg/problems/dgproblem4.hh
+ 13
− 17
View file @ 99ed2a23
......@@ -12,20 +12,16 @@ class DGProblem4 : public DGProblem<Coordinate, Scalar>
using Base = DGProblem<Coordinate, Scalar>;
using Matrix = typename Base::Matrix;
// static constexpr Scalar dmin_ = 0.05; //minimum aperture of the fracture
// static constexpr Scalar dmax_ = 0.2; //maximum aperture of the fracture
//require dim == 2
DGProblem4 (const Scalar dmin, const Scalar dmax)
: dmin_(dmin), dmax_(dmax) {}
//{
// if (dim != 2)
// {
// DUNE_THROW(Dune::Exception,
// "Problem dg4 is not implemented for dimension = " +
// std::to_string(dim) + "!");
// }
//}
//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
......@@ -118,21 +114,21 @@ class DGProblem4 : public DGProblem<Coordinate, Scalar>
//returns the recommended quadrature order to compute an integral
//over x * q(x)
virtual int quadratureOrder_q () const
int quadratureOrder_q () const
{
return 10;
}
//returns the recommended quadrature order to compute an integral
//over x * boundary(x) and K(x) * boundary(x)
virtual int quadratureOrderBoundary () const
int quadratureOrderBoundary () const
{
return 10;
}
//returns the recommended quadrature order to compute an integral
//over an entry K[i][j] of the permeability tensor
virtual int quadratureOrder_K () const
int quadratureOrder_K () const
{
return 10;
}
......
dune/mmdg/problems/dgproblem5.hh
+ 44
− 43
View file @ 99ed2a23
......@@ -12,22 +12,21 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
using Base = DGProblem<Coordinate, Scalar>;
using Matrix = typename Base::Matrix;
static constexpr Scalar dmin = 0.05; //minimum aperture of the fracture
static constexpr Scalar dmax = 0.1; //maximum aperture of the fracture
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
//require dim == 2
DGProblem5 ()
{
if (dim != 2)
{
DUNE_THROW(Dune::Exception,
"Problem dg4 is not implemented for dimension = " +
std::to_string(dim) + "!");
}
}
//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) + "!");
// }
// }
//the exact solution at position pos
Scalar exactSolution (const Coordinate& pos) const
......@@ -57,24 +56,20 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
{
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 ) *
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 ) *
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];
// return -2.0 * K_perp * pos[1] - ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) *
// ( -1.0 + 4.0 * M_PI * ( dmax - dmin ) * ( 8.0 * M_PI *
// std::cos(8.0 * M_PI * pos[1]) - std::sin(8.0 * M_PI * pos[1]) )
// / aperture(pos) );
}
//permeability tensor at position pos
......@@ -88,12 +83,6 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
for (int i = 0; i < dim; i++)
{
permeability[i][i] = Kbulk;
// aperture(pos) / ( aperture(pos) - 4.0 * pos[1] * M_PI * (dmax - dmin) * std::sin(8.0 * M_PI * pos[1]) );
if (permeability[i][i] <= 0.0)
{
std::cout << "WARNING!" << std::endl;
}
}
}
else
......@@ -103,37 +92,49 @@ class DGProblem5 : public DGProblem<Coordinate, Scalar>
for (int i = 1; i < dim; i++)
{
permeability[i][i] = K_par;
// 1.0 - pos[1] * 4.0 * M_PI * ( dmax - dmin ) *
// std::sin(8.0 * M_PI * pos[1]) / aperture(pos);
if (permeability[i][i] <= 0.0)
{
std::cout << "WARNING" << std::endl;
}
}
}
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 ) *
return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
//return dmin + ( dmax - dmin ) * std::cos(8.0 * M_PI * pos[1]);
// const Scalar y = pos[1] - 0.25 * floor(4.0 * pos[1]);
//
// if ( y <= 0.125 )
// { //0 <= y <= 0.125
// return dmin + ( 1.0 - 8.0 * y ) * ( dmax - dmin );
// }
//
// //0.125 <= y <= 0.25
// return dmin + ( 8.0 * y - 1.0 ) * ( dmax - dmin );
}
//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
dune/mmdg/problems/mmdgproblem7.hh
+ 7
− 13
View file @ 99ed2a23
......@@ -13,14 +13,9 @@ public:
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) {}
MMDGProblem7 (const Scalar dmin, const Scalar dmax)
: dmin_(dmin), dmax_(dmax) {}
//the exact bulk solution at position pos
Scalar exactSolution (const Coordinate& pos) const
......@@ -78,7 +73,7 @@ public:
//aperture d of the fracture at position pos
Scalar aperture (const Coordinate& pos) const
{
return dmin + 0.5 * ( dmax - dmin ) *
return dmin_ + 0.5 * ( dmax_ - dmin_ ) *
( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
}
......@@ -86,7 +81,7 @@ public:
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_ ) * std::sin(8.0 * M_PI * pos[1]);
return gradD;
}
......@@ -182,16 +177,15 @@ public:
}
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 )
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
};
......
src/dgAnalysis.py
+ 20
− 8
View file @ 99ed2a23
......@@ -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
......
src/mmdgAnalysis.py
+ 1
− 1
View file @ 99ed2a23
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment