From 99ed2a2371b6b33f151c53e78cdc573eb2dab347 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 22:00:18 +0200
Subject: [PATCH] add dmin, dmax as member variables; adapt analysis scripts to
new problems
---
dune/mmdg/problems/dgproblem4.hh | 30 +++++------
dune/mmdg/problems/dgproblem5.hh | 87 +++++++++++++++---------------
dune/mmdg/problems/mmdgproblem7.hh | 20 +++----
src/dgAnalysis.py | 28 +++++++---
src/mmdgAnalysis.py | 2 +-
5 files changed, 85 insertions(+), 82 deletions(-)
diff --git a/dune/mmdg/problems/dgproblem4.hh b/dune/mmdg/problems/dgproblem4.hh
index d00936a..7ea27cb 100644
--- a/dune/mmdg/problems/dgproblem4.hh
+++ b/dune/mmdg/problems/dgproblem4.hh
@@ -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;
}
diff --git a/dune/mmdg/problems/dgproblem5.hh b/dune/mmdg/problems/dgproblem5.hh
index 1a6e4d5..880f6b9 100644
--- a/dune/mmdg/problems/dgproblem5.hh
+++ b/dune/mmdg/problems/dgproblem5.hh
@@ -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
diff --git a/dune/mmdg/problems/mmdgproblem7.hh b/dune/mmdg/problems/mmdgproblem7.hh
index 8477a75..19b8984 100644
--- a/dune/mmdg/problems/mmdgproblem7.hh
+++ b/dune/mmdg/problems/mmdgproblem7.hh
@@ -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
};
diff --git a/src/dgAnalysis.py b/src/dgAnalysis.py
index 4a6edd2..a54c24d 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/mmdgAnalysis.py b/src/mmdgAnalysis.py
index 4aecbc6..d403120 100644
--- a/src/mmdgAnalysis.py
+++ b/src/mmdgAnalysis.py
@@ -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
--
GitLab