diff --git a/README.md b/README.md index b4736659764719d1674dcef972a01035e3f04046..af64cb26b9cd61d72db8781703d6ea3bea7c8ba4 100644 --- a/README.md +++ b/README.md @@ -51,8 +51,8 @@ You will find the plots `mmdg_bulk.pdf`, `mmdg_interface.pdf` and `mmdg_total.pd Test Problems ------------- -The `dune-mmdg` module includes three test problems, `dg1` to `dg3`, for the bulk scheme and six test problems, `mmdg1` to `mmdg6`, for the coupled scheme. -The test problems `dg2`, `mmdg2`, `mmdg5` and `mmdg6` correspond to the examples 5.1 to 5.4 in [1]. +The `dune-mmdg` module includes five test problems, `dg1` to `dg5`, for the bulk scheme and seven test problems, `mmdg1` to `mmdg7`, for the coupled scheme. +The test problems `dg2`, `mmdg2`, `mmdg5`, `mmdg6` and `dg4`/`mmdg7` correspond to the examples 5.1, 5.2, 5.3, 5.4/5.5 and 5.6 in [1] in the given order. Parameters @@ -61,16 +61,21 @@ Parameters For the bulk scheme the following parameters can be set in the file `parameterDG.ini` in the `src` directory: - `mu0`: the prefactor of the penalty parameter (must be sufficiently large) -- `problem`: the name of a test problem (`dg1` to `dg3`) +- `problem`: the name of a test problem (`dg1` to `dg5`) +- `dmax`: the maximum value of the aperture function (only relevant for the problems `dg4` and `dg5` that include a fracture) +- `dmin`: the minimum value of the aperture function (only relevant for the problems `dg4` and `dg5` that include a fracture) +- `solver`: (optional) the name of the solver that is applied to the system of linear equations (`Cholmod` or `UMFPack`) - `fastEOC`: (optional) exclude computationally expensive grids from the EOC loop For the coupled scheme the following parameters can be set in the file `parameterMMDG.ini` in the `src` directory: - `mu0`: the prefactor of the penalty parameter (must be sufficiently large) -- `problem`: the name of a test problem (`mmdg1` to `mmdg6`) -- `d`: the aperture of the fracture or the prefactor of the aperture function +- `problem`: the name of a test problem (`mmdg1` to `mmdg7`) +- `dmax`: the maximum value of the aperture of the fracture or respective constant aperture value (mandatory for all problems) +- `dmin`: the minimum value of the aperture of the fracture (only relevant for problem `mmdg7`) - `xi`: (optional) the value of the coupling parameter (must be larger than 0.5) - `gridfile`: (optional) the prefix of a grid file name in the directory `src/grids` +- `solver`: (optional) the name of the solver that is applied to the system of linear equations (`Cholmod` or `UMFPack`) - `fastEOC`: (optional) exclude computationally expensive grids from the EOC loop diff --git a/dune/mmdg/clockhelper.hh b/dune/mmdg/clockhelper.hh index 05c3ddc2afd28bba0deb357fc5f9edbcc0ca4790..850e1e49d0e82ca344ceea4ee4a780befd89f38c 100644 --- a/dune/mmdg/clockhelper.hh +++ b/dune/mmdg/clockhelper.hh @@ -3,6 +3,7 @@ #include <chrono> +//contains a function for run time measurements class ClockHelper { public: @@ -12,8 +13,8 @@ class ClockHelper //precision of a microsecond static double elapsedRuntime(const Clock::time_point begin) { - return 1e-6*std::chrono::duration_cast<std::chrono::microseconds> - (Clock::now() - begin).count();; + return 1e-6 * std::chrono::duration_cast<std::chrono::microseconds> + ( Clock::now() - begin ).count();; } }; diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index cdc328742b4edce60906f84def056a8db414b6df..8a00b9bc87c1205c1abb56a5b8e237ad5f8ef674 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -48,7 +48,7 @@ public: dof_((1 + dim) * gridView.size(0)), solver_(solver) { //initialize stiffness matrix A, load vector b and solution vector d - A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO + A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); b = Vector(dof_); d = Vector(dof_); } @@ -59,6 +59,7 @@ public: assembleSLE(mu0); A->compress(); + //solve system of linear equations Dune::InverseOperatorResult result; if (solver_ == SolverType::Cholmod) @@ -82,7 +83,7 @@ public: writeVTKoutput(); } - //compute L2 error using a quadrature rule + //computes the L2 error using a quadrature rule Scalar computeL2error(const int quadratureOrder = 5) { if (!problem_.hasExactSolution()) @@ -129,13 +130,14 @@ public: const int dof_; //degrees of freedom protected: + //constructor with custom value of dof_ DG (const GridView& gridView, const Mapper& mapper, const Problem& problem, const SolverType solver, const int dof) : gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof), solver_(solver) { //initialize stiffness matrix A, load vector b and solution vector d - A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO + A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); b = Vector(dof_); d = Vector(dof_); } @@ -169,7 +171,7 @@ protected: //in the system of linear equations (SLE) Ad = b, //the index elemIdxSLE refers to the basis function phi_elem,0 //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the - //basis function phi_elem,i + //basis functions phi_elem,i const int elemIdxSLE = (dim + 1)*elemIdx; //lambda to fill load vector b using a quadrature rule @@ -230,7 +232,7 @@ protected: //create storage for multiply used integral values //linearIntegrals[i] = int_intersct x[i] ds - //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds //NOTE: is there a better type for a symmetric matrix? + //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds //KIntegrals[i] = int_intersct K[:,i]*normal ds //KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds Coordinate linearIntegrals(0.0); @@ -285,7 +287,7 @@ protected: } }; - //apply quadrature rules in intersct + //apply quadrature rules on intersct applyQuadratureRule(intersctGeo, rule_2ndOrder, lambda_2ndOrder); applyQuadratureRule(intersctGeo, rule_K_Edge, lambda_K_Edge); applyQuadratureRule(intersctGeo, rule_Kplus1_Edge, lambda_Kplus1_Edge); diff --git a/dune/mmdg/eigenvaluehelper.hh b/dune/mmdg/eigenvaluehelper.hh index a3e9cf0110cf2db8cb68bd30bea57c92d60f4a14..7fe56886e13b5457466458d04b91826281c81655 100644 --- a/dune/mmdg/eigenvaluehelper.hh +++ b/dune/mmdg/eigenvaluehelper.hh @@ -10,7 +10,7 @@ /** \file * \brief Eigenvalue computations for the Dune::FieldMatrix class - * NOTE: taken from DUNE-common, fmatrixev.hh, with adapted tolerance + * NOTE: taken from DUNE-common (fmatrixev.hh) with minor changes */ class EigenvalueHelper { diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index 70d6a3100c3b1800131f38fe353559e954c89992..0889dafd7bd1519bd9a1b3c4de3be1d53bebee0b 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -37,9 +37,11 @@ public: const void operator() (const Scalar mu0, const Scalar xi) { + //assemble stiffness matrix A and load vector b assembleSLE(mu0, xi); Base::A->compress(); + //solve system of linear equations Dune::InverseOperatorResult result; if (Base::solver_ == SolverType::Cholmod) @@ -65,7 +67,7 @@ public: //compute L2 errors using quadrature rules, //returns {error(bulk), error(interface), error(total)}, - //the total error is calculated using the norm + //the total error is calculated with respect to the norm // norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) ) const std::array<Scalar, 3> computeL2error(const int quadratureOrder = 5) { @@ -104,8 +106,8 @@ public: Base::applyQuadratureRule(iGeo, qrule, qlambda); } - return {bulkError, sqrt(interfaceError), - sqrt(bulkError * bulkError + interfaceError)}; + return { bulkError, sqrt(interfaceError), + sqrt(bulkError * bulkError + interfaceError) }; } const InterfaceGridView& iGridView_; @@ -161,6 +163,7 @@ private: //and load vector b with bulk entries Base::assembleSLE(mu0); + //fill in coupling and interface entries for (const auto& iElem : elements(iGridView_)) { const int iElemIdx = iMapper_.index(iElem); @@ -892,6 +895,8 @@ private: return iBasis; } + //write exakt and numerical bulk and interface solution as well as the + //aperture function of the problem to vtk files void writeVTKoutput () const { //degrees of freedom of bulk and interface grids @@ -929,6 +934,7 @@ private: Base::problem_.exactInterfaceSolution(iGeo.corner(k)); } + //aperture evaluated at the kth corner of iElem aperture[iElemIdxOutput + k] = Base::problem_.aperture(iGeo.corner(k)); diff --git a/dune/mmdg/nonconformingp1vtkfunction.hh b/dune/mmdg/nonconformingp1vtkfunction.hh index 9ae9ec80968a149730263320846d4e4715cbfb0e..45a69ad7582d8e8af9245c9ff7c33063ad3ba4f8 100644 --- a/dune/mmdg/nonconformingp1vtkfunction.hh +++ b/dune/mmdg/nonconformingp1vtkfunction.hh @@ -1,8 +1,8 @@ // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- // vi: set et ts=4 sw=2 sts=2: -#ifndef DUNE_MMFD_FUNCTION_HH -#define DUNE_MMFD_FUNCTION_HH +#ifndef DUNE_MMDG_FUNCTION_HH +#define DUNE_MMDG_FUNCTION_HH #include <string> @@ -17,6 +17,7 @@ #include <dune/grid/io/file/vtk/common.hh> #include <dune/grid/io/file/vtk/function.hh> +//NOTE: without changes from DUNE-MMFD by Samuel Burbulla namespace Dune { ////////////////////////////////////////////////////////////////////// @@ -143,4 +144,4 @@ namespace Dune } // namespace Dune -#endif // DUNE_MMFD_FUNCTION_HH +#endif // DUNE_MMDG_FUNCTION_HH diff --git a/dune/mmdg/solvertype.hh b/dune/mmdg/solvertype.hh index 560113e31015fd3a8287a27d4f462c101c20678b..6ca83002e71a6f558eeab4ad7036f14ab47a2eb1 100644 --- a/dune/mmdg/solvertype.hh +++ b/dune/mmdg/solvertype.hh @@ -1,3 +1,4 @@ +//determines the type of solver that is applied to a system of linear equations enum SolverType { Cholmod, diff --git a/src/parameterDG.ini b/src/parameterDG.ini index 11dd074d59ed026ad20d0ac76224c67cae72225c..c7ee27f47aab88515617404f860f61d66f9e2618 100644 --- a/src/parameterDG.ini +++ b/src/parameterDG.ini @@ -1,6 +1,6 @@ mu0 = 1000 problem = dg4 -dmin = 1e-2 dmax = 1e-1 +dmin = 1e-2 solver = Cholmod fastEOC = true