From f1c0cbf5651265cc9d982858e76670cc3e5030a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 3 Apr 2020 13:23:06 +0200
Subject: [PATCH] revise quadrature orders for mmdg problems, rename dg
problems
---
dune/mmdg/mmdg.hh | 2 --
dune/mmdg/problems/coupleddgproblem.hh | 18 ++++++++-----
.../problems/{poisson.hh => dgproblem1.hh} | 6 ++---
...omogeneousbulkproblem.hh => dgproblem2.hh} | 6 ++---
.../{zeropoisson.hh => dgproblem3.hh} | 9 +++----
dune/mmdg/problems/mmdgproblem1.hh | 7 -----
dune/mmdg/problems/mmdgproblem3.hh | 14 ----------
dune/mmdg/problems/mmdgproblem4.hh | 22 ++++------------
dune/mmdg/problems/mmdgproblem5.hh | 21 ++++++++++++++-
src/dg.cc | 26 ++++++++++++-------
src/parameterDG.ini | 2 +-
11 files changed, 64 insertions(+), 69 deletions(-)
rename dune/mmdg/problems/{poisson.hh => dgproblem1.hh} (92%)
rename dune/mmdg/problems/{inhomogeneousbulkproblem.hh => dgproblem2.hh} (90%)
rename dune/mmdg/problems/{zeropoisson.hh => dgproblem3.hh} (88%)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index c2040e3..32c5953 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -186,7 +186,6 @@ private:
// -int_iElem beta * phi_iElem,j * avg(phi_elem,0) ds
//for i = 0,...,dim and j = 0,...,dim-1
//and for elem in {elemIn, elemOut} using a quadrature rule
- //TODO symmetrize more efficiently
const auto lambda_Kperp =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
@@ -292,7 +291,6 @@ private:
// -int_iElem beta * phi_iElem,i * avg(phi_elem1,k) ds
//for i,j = 1,...,dim-1 and k,l = 1,...,dim
//and for elem1,elem2 in {elemIn, elemOut} using a quadrature rule
- //TODO symmetrize more efficiently
const auto lambda_KperpPlus1 =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index ed9ad18..9f7774f 100644
--- a/dune/mmdg/problems/coupleddgproblem.hh
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -57,7 +57,8 @@ public:
return permeability;
}
- //permeability per aperture of the fracture in normal direction at position pos
+ //permeability per aperture of the fracture in normal direction
+ //at position pos
virtual Scalar Kperp (const Coordinate& pos) const
{
return Scalar(1.0);
@@ -71,21 +72,26 @@ public:
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
virtual int quadratureOrderInterfaceBoundary () const
{
return 1;
}
- //returns the recommended quadrature order to compute an integral
- //over Kparallel(x)
+ //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
virtual int quadratureOrder_Kparallel () const
{
- return 1;
+ return 2;
}
//returns the recommended quadrature order to compute an integral
- //over x * Kperp(x) / d(x)
+ //over x * Kperp(x)
virtual int quadratureOrder_Kperp () const
{
return 1;
diff --git a/dune/mmdg/problems/poisson.hh b/dune/mmdg/problems/dgproblem1.hh
similarity index 92%
rename from dune/mmdg/problems/poisson.hh
rename to dune/mmdg/problems/dgproblem1.hh
index 135dc58..39f6748 100644
--- a/dune/mmdg/problems/poisson.hh
+++ b/dune/mmdg/problems/dgproblem1.hh
@@ -1,11 +1,11 @@
-#ifndef __DUNE_MMDG_POISSON_HH__
-#define __DUNE_MMDG_POISSON_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM1_HH__
+#define __DUNE_MMDG_DGPROBLEM1_HH__
#include <dune/mmdg/problems/dgproblem.hh>
//a Poisson problem
template<class Coordinate, class Scalar = double>
-class Poisson : public DGProblem<Coordinate, Scalar>
+class DGProblem1 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
diff --git a/dune/mmdg/problems/inhomogeneousbulkproblem.hh b/dune/mmdg/problems/dgproblem2.hh
similarity index 90%
rename from dune/mmdg/problems/inhomogeneousbulkproblem.hh
rename to dune/mmdg/problems/dgproblem2.hh
index 601d0bb..5d82404 100644
--- a/dune/mmdg/problems/inhomogeneousbulkproblem.hh
+++ b/dune/mmdg/problems/dgproblem2.hh
@@ -1,10 +1,10 @@
-#ifndef __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
-#define __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM2_HH__
+#define __DUNE_MMDG_DGPROBLEM2_HH__
#include <dune/mmdg/problems/dgproblem.hh>
template<class Coordinate, class Scalar = double>
-class InhomogeneousBulkProblem : public DGProblem<Coordinate, Scalar>
+class DGProblem2 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
diff --git a/dune/mmdg/problems/zeropoisson.hh b/dune/mmdg/problems/dgproblem3.hh
similarity index 88%
rename from dune/mmdg/problems/zeropoisson.hh
rename to dune/mmdg/problems/dgproblem3.hh
index 56525ab..2942138 100644
--- a/dune/mmdg/problems/zeropoisson.hh
+++ b/dune/mmdg/problems/dgproblem3.hh
@@ -1,15 +1,14 @@
-#ifndef __DUNE_MMDG_ZEROPOISSON_HH__
-#define __DUNE_MMDG_ZEROPOISSON_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM3_HH__
+#define __DUNE_MMDG_DGPROBLEM3_HH__
#include <cmath>
#include <string>
#include <dune/mmdg/problems/dgproblem.hh>
-
//Poisson problem with right side q = -1
template<class Coordinate, class Scalar = double>
-class ZeroPoisson : public DGProblem<Coordinate, Scalar>
+class DGProblem3 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
@@ -19,7 +18,7 @@ class ZeroPoisson : public DGProblem<Coordinate, Scalar>
{
if constexpr (dim == 1)
{
- return 0.5*pos[0]*(pos[0] - 1);
+ return 0.5 * pos[0] * (pos[0] - 1);
}
if constexpr (dim == 2) //approximate series
diff --git a/dune/mmdg/problems/mmdgproblem1.hh b/dune/mmdg/problems/mmdgproblem1.hh
index 6ea0cd8..f35a61f 100644
--- a/dune/mmdg/problems/mmdgproblem1.hh
+++ b/dune/mmdg/problems/mmdgproblem1.hh
@@ -47,13 +47,6 @@ public:
return 2;
}
- //returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
- int quadratureOrderInterfaceBoundary () const
- {
- return 2;
- }
-
private:
Scalar d_; //aperture
};
diff --git a/dune/mmdg/problems/mmdgproblem3.hh b/dune/mmdg/problems/mmdgproblem3.hh
index a6914da..d3f284c 100644
--- a/dune/mmdg/problems/mmdgproblem3.hh
+++ b/dune/mmdg/problems/mmdgproblem3.hh
@@ -79,20 +79,6 @@ public:
return 10;
}
- //returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
- int quadratureOrderInterfaceBoundary () const
- {
- return 10;
- }
-
- //returns the recommended quadrature order to compute an integral
- //over x * qInterface(x)
- int quadratureOrder_qInterface () const
- {
- return 10;
- }
-
private:
Scalar d_; //aperture
};
diff --git a/dune/mmdg/problems/mmdgproblem4.hh b/dune/mmdg/problems/mmdgproblem4.hh
index 1cae124..cf09ef6 100644
--- a/dune/mmdg/problems/mmdgproblem4.hh
+++ b/dune/mmdg/problems/mmdgproblem4.hh
@@ -76,28 +76,16 @@ public:
//over x * boundary(x)
int quadratureOrderBoundary () const
{
- return 10;
+ return 3;
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //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;
+ return 3;
}
private:
diff --git a/dune/mmdg/problems/mmdgproblem5.hh b/dune/mmdg/problems/mmdgproblem5.hh
index 34ade2f..ee61489 100644
--- a/dune/mmdg/problems/mmdgproblem5.hh
+++ b/dune/mmdg/problems/mmdgproblem5.hh
@@ -100,7 +100,9 @@ public:
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //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;
@@ -120,6 +122,23 @@ public:
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 5;
+ }
+
+ //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:
Scalar d_; //aperture
diff --git a/src/dg.cc b/src/dg.cc
index dcc4659..1b21492 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -21,9 +21,9 @@
#include <dune/mmdg/dg.hh>
#include <dune/mmdg/clockhelper.hh>
#include <dune/mmdg/problems/dgproblem.hh>
-#include <dune/mmdg/problems/zeropoisson.hh>
-#include <dune/mmdg/problems/poisson.hh>
-#include <dune/mmdg/problems/inhomogeneousbulkproblem.hh>
+#include <dune/mmdg/problems/dgproblem1.hh>
+#include <dune/mmdg/problems/dgproblem2.hh>
+#include <dune/mmdg/problems/dgproblem3.hh>
int main(int argc, char** argv)
{
@@ -65,23 +65,29 @@ int main(int argc, char** argv)
Problem* problem; //problem to be solved
//determine problem type
- if (pt["problem"] == "zeropoisson")
+ if (pt["problem"] == "dg1")
{
- problem = new ZeroPoisson<Coordinate>();
+ problem = new DGProblem1<Coordinate>();
}
- else if (pt["problem"] == "poisson")
+ else if (pt["problem"] == "dg2")
{
- problem = new Poisson<Coordinate>();
+ problem = new DGProblem2<Coordinate>();
}
- else if (pt["problem"] == "inhomogeneousbulk")
+ else if (pt["problem"] == "dg3")
{
- problem = new InhomogeneousBulkProblem<Coordinate>();
+ problem = new DGProblem3<Coordinate>();
+
+ if constexpr (dim == 2)
+ {
+ std::cout << "NOTE: The exact solution of your problem is given by "
+ "an infinite series; increase the cut-off if necessary!\n\n";
+ }
}
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'poisson' or TODO)");
+ "parameterDG.ini (use 'dg1' to 'dg3')");
}
DG dg(gridView, mapper, *problem);
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 10896d3..e3e1b98 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,2 +1,2 @@
mu0 = 1000
-problem = poisson
+problem = dg1
--
GitLab