diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh index e6f8d6e2df9af1bd123c53532f437f345d094f57..e6fe7a41834d59102a4dc51fb63325cc2cf281c8 100644 --- a/dune/mmdg/problems/coupleddgproblem.hh +++ b/dune/mmdg/problems/coupleddgproblem.hh @@ -10,7 +10,7 @@ template<class Vector, class Scalar = double> class CoupledDGProblem : public DGProblem<Vector, Scalar> { public: - static constexpr int dimension = Vector::dimension; + static constexpr int dim = Vector::dimension; using Base = DGProblem<Vector, Scalar>; using Matrix = typename Base::Matrix; @@ -37,7 +37,7 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar> { Matrix permeability(0.0); - for (int i = 0; i < dimension; i++) + for (int i = 0; i < dim; i++) { permeability[i][i] = 1.0; } diff --git a/dune/mmdg/problems/dgproblem.hh b/dune/mmdg/problems/dgproblem.hh index 982015eaf089143a0fd010dba1fb30f72efa75d5..13d0bc4be19c808f9cbcbff62d3c0a2e41c272ec 100644 --- a/dune/mmdg/problems/dgproblem.hh +++ b/dune/mmdg/problems/dgproblem.hh @@ -3,15 +3,15 @@ //abstract class providing an interface for a problem that is to be solved with //a discontinuous Galerkin scheme -template<class Vector, class Scalar = double> +template<class Coordinate, class Scalar = double> class DGProblem { public: - static constexpr int dimension = Vector::dimension; - using Matrix = Dune::FieldMatrix<Scalar, dimension, dimension>; + static constexpr int dim = Coordinate::dimension; + using Matrix = Dune::FieldMatrix<Scalar, dim, dim>; //the exact solution at position pos - virtual Scalar exactSolution (const Vector& pos) const + virtual Scalar exactSolution (const Coordinate& pos) const { return Scalar(0.0); } @@ -23,23 +23,23 @@ class DGProblem }; //source term at position pos - virtual Scalar q (const Vector& pos) const + virtual Scalar q (const Coordinate& pos) const { return Scalar(0.0); } //boundary condition at position pos - virtual Scalar boundary (const Vector& pos) const + virtual Scalar boundary (const Coordinate& pos) const { return Scalar(0.0); } //permeability tensor at position pos - virtual Matrix K (const Vector& pos) const + virtual Matrix K (const Coordinate& pos) const { Matrix permeability(0.0); - for (int i = 0; i < dimension; i++) + for (int i = 0; i < dim; i++) { permeability[i][i] = 1.0; } diff --git a/dune/mmdg/problems/inhomogeneousbulkproblem.hh b/dune/mmdg/problems/inhomogeneousbulkproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..8cf6424e4280d900b1bb4028a7947b84a2c763ac --- /dev/null +++ b/dune/mmdg/problems/inhomogeneousbulkproblem.hh @@ -0,0 +1,76 @@ +#ifndef __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__ +#define __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__ + +#include <dune/mmdg/problems/dgproblem.hh> + +template<class Coordinate, class Scalar = double> +class InhomogeneousBulkProblem : public DGProblem<Coordinate, Scalar> +{ + public: + static constexpr int dim = Coordinate::dimension; + using Base = DGProblem<Coordinate, Scalar>; + using Matrix = typename Base::Matrix; + + //the exact solution at position pos + //1d: f(x) = x, + //2d: f(x,y) = x*y, + //3d: f(x,y,z) = x*y*z + virtual Scalar exactSolution (const Coordinate& pos) const + { + return pos * Coordinate(1.0); + } + + //exact solution is implemented + bool hasExactSolution() const + { + return true; + }; + + //source term at position pos + virtual Scalar q (const Coordinate& pos) const + { + return pos * Coordinate(2.0); + } + + //boundary condition + virtual Scalar boundary (const Coordinate& pos) const + { + return exactSolution(pos); + } + + //permeability tensor at position pos + virtual Matrix K (const Coordinate& pos) const + { + Matrix permeability(0.0); + + for (int i = 0; i < dim; i++) + { + permeability[i][i] = 1.0 + pos[i] * pos[i]; + } + + return permeability; + } + + //returns the recommended quadrature order to compute an integral + //over x * q(x) + virtual int quadratureOrder_q () const + { + return 10; + } + + //returns the recommended quadrature order to compute an integral + //over x * boundary(x) + virtual 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 + { + return 10; + } +}; + +#endif diff --git a/dune/mmdg/problems/poisson.hh b/dune/mmdg/problems/poisson.hh index 8781cb38180aae3584f1f48f517da1420d7a9c5c..e6f375743be6428ae56f51844c14bc427cde7453 100644 --- a/dune/mmdg/problems/poisson.hh +++ b/dune/mmdg/problems/poisson.hh @@ -1,65 +1,39 @@ #ifndef __DUNE_MMDG_POISSON_HH__ #define __DUNE_MMDG_POISSON_HH__ -#include <cmath> -#include <string> - #include <dune/mmdg/problems/dgproblem.hh> - -//Poisson problem with right side q = -1 -template<int dim, class Scalar = double> -class Poisson : public DGProblem<Dune::FieldVector<Scalar, dim>, Scalar> +//a Poisson problem +template<class Coordinate, class Scalar = double> +class Poisson : public DGProblem<Coordinate, Scalar> { public: - using Vector = Dune::FieldVector<Scalar, dim>; - //the exact solution at position pos - virtual Scalar exactSolution (const Vector& pos) const + //1d: f(x) = x, + //2d: f(x,y) = x*y, + //3d: f(x,y,z) = x*y*z + virtual Scalar exactSolution (const Coordinate& pos) const { - if (dim == 1) - { - return 0.5*pos[0]*(pos[0] - 1); - } + Scalar solution(1.0); - if (dim == 2) + for (int i = 0; i < pos.dim(); i++) { - const int cutoff = 21; //NOTE: error? - double solution = 0.0; - for (int m = 1; m <= cutoff; m = m + 2) - { - for (int n = 1; n <= cutoff; n = n + 2) - { - solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1]) - / (m * n * (m*m + n*n)); - } - } - solution *= -16 / pow(M_PI, 4); - return solution; + solution *= pos[i]; } - DUNE_THROW(Dune::Exception, - "Exact Solution not implemented for dimension = " + - std::to_string(dim) + "!"); + return solution; } - //exact solution is implemented for dim = 1, 2 + //exact solution is implemented bool hasExactSolution() const { - return (dim == 1 || dim == 2); + return true; }; - //source term - Scalar q (const Vector& pos) const - { - return Scalar(-1.0); - } - - //NOTE: necessary or already default? - //homogeneous Dirichlet boundary conditions - virtual Scalar boundary (const Vector& pos) const + //boundary condition + virtual Scalar boundary (const Coordinate& pos) const { - return Scalar(0.0); + return exactSolution(pos); } }; diff --git a/dune/mmdg/problems/zeropoisson.hh b/dune/mmdg/problems/zeropoisson.hh new file mode 100644 index 0000000000000000000000000000000000000000..dce45f35a7ee00f8aef960c9ae49cb798477c82d --- /dev/null +++ b/dune/mmdg/problems/zeropoisson.hh @@ -0,0 +1,59 @@ +#ifndef __DUNE_MMDG_ZEROPOISSON_HH__ +#define __DUNE_MMDG_ZEROPOISSON_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> +{ + public: + static constexpr int dim = Coordinate::dimension; + + //the exact solution at position pos + virtual Scalar exactSolution (const Coordinate& pos) const + { + if constexpr (dim == 1) + { + return 0.5*pos[0]*(pos[0] - 1); + } + + if constexpr (dim == 2) //approximate series + { + const int cutoff = 21; + double solution = 0.0; + for (int m = 1; m <= cutoff; m = m + 2) + { + for (int n = 1; n <= cutoff; n = n + 2) + { + solution += sin(M_PI * m * pos[0]) * sin(M_PI * n * pos[1]) + / (m * n * (m*m + n*n)); + } + } + solution *= -16 / pow(M_PI, 4); + return solution; + } + + DUNE_THROW(Dune::Exception, + "Exact Solution not implemented for dimension = " + + std::to_string(dim) + "!"); + } + + //exact solution is implemented for dim = 1, 2 + bool hasExactSolution() const + { + return (dim == 1 || dim == 2); + }; + + //source term + Scalar q (const Coordinate& pos) const + { + return Scalar(-1.0); + } +}; + +#endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5088d1777b4549558220e82eaca0ea1d0c420a14..7f1cc3f048a06cc669ac8b46b789986ed1318198 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,10 @@ add_executable("dg-2d" dg.cc) target_compile_definitions("dg-2d" PRIVATE GRIDDIM=2) target_link_dune_default_libraries("dg-2d") +add_executable("dg-3d" dg.cc) +target_compile_definitions("dg-3d" PRIVATE GRIDDIM=3) +target_link_dune_default_libraries("dg-3d") + add_executable("mmdg-2d" mmdg.cc) target_compile_definitions("mmdg-2d" PRIVATE GRIDDIM=2) target_link_dune_default_libraries("mmdg-2d") diff --git a/src/dg.cc b/src/dg.cc index 9f7f07c5cecb2f2a0ccacd6cf148c512350d824c..4b55a715a6c225af670d92c8deba0fa348978b58 100644 --- a/src/dg.cc +++ b/src/dg.cc @@ -20,8 +20,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> int main(int argc, char** argv) { @@ -62,9 +63,17 @@ int main(int argc, char** argv) Problem* problem; //problem to be solved //determine problem type - if (pt["problem"] == "poisson") + if (pt["problem"] == "zeropoisson") + { + problem = new ZeroPoisson<Coordinate>(); + } + else if (pt["problem"] == "poisson") + { + problem = new Poisson<Coordinate>(); + } + else if (pt["problem"] == "inhomogeneousbulk") { - problem = new Poisson<dim>(); + problem = new InhomogeneousBulkProblem<Coordinate>(); } else { @@ -81,7 +90,7 @@ int main(int argc, char** argv) //print total run time const double timeTotal = ClockHelper::elapsedRuntime(timeInitial); std::cout << "Finished with a total run time of " << timeTotal << - " Seconds.\n"; + " Seconds.\nL2-error: " << dg.computeL2error() << ".\n"; return 0; } diff --git a/src/grids/cube1d.dgf b/src/grids/cube1d.dgf index 7f8f178bd159e56cfd5727de8c96e9ea0b3bc002..0121fdffd456169544cb644f3b2ddd1c53800150 100644 --- a/src/grids/cube1d.dgf +++ b/src/grids/cube1d.dgf @@ -3,7 +3,7 @@ DGF Interval 0 1 -100 +1000 # Simplex diff --git a/src/grids/cube2d.dgf b/src/grids/cube2d.dgf index 4ef9b6a3b250bd9dc9feb9ca1a4f0bce49e533a7..d733edc27b154a3c085cb3085a4988315e05a7cc 100644 --- a/src/grids/cube2d.dgf +++ b/src/grids/cube2d.dgf @@ -3,7 +3,7 @@ DGF Interval 0 0 1 1 -20 20 +100 100 # Simplex diff --git a/src/grids/cube3d.dgf b/src/grids/cube3d.dgf index d201d0aa39076bf723bd42eeadf5aebf2432d3f7..0d7af180f32372f9a4f700cef1f5368606a29934 100644 --- a/src/grids/cube3d.dgf +++ b/src/grids/cube3d.dgf @@ -3,7 +3,7 @@ DGF Interval 0 0 0 1 1 1 -20 20 20 +10 10 10 # Simplex