From 37e6bf9abb791ecb12acb857f791fa34430ab682 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 21:58:14 +0200
Subject: [PATCH] add enum SolverType, add parameters solver, dmin, dmax
---
dune/mmdg/dg.hh | 33 +++++++++++++++++++++------
dune/mmdg/mmdg.hh | 24 +++++++++++++++-----
dune/mmdg/solvertype.hh | 5 +++++
src/dg.cc | 43 +++++++++++++++++++++++++++++++++--
src/mmdg.cc | 50 +++++++++++++++++++++++++++++++----------
src/parameterDG.ini | 5 +++--
src/parameterMMDG.ini | 6 +++--
7 files changed, 136 insertions(+), 30 deletions(-)
create mode 100644 dune/mmdg/solvertype.hh
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 623c825..cdc3287 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -14,9 +14,11 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/cholmod.hh>
+#include <dune/istl/umfpack.hh>
#include <dune/mmdg/nonconformingp1vtkfunction.hh>
#include <dune/mmdg/eigenvaluehelper.hh>
+#include <dune/mmdg/solvertype.hh>
template<class GridView, class Mapper, class Problem>
class DG
@@ -40,9 +42,10 @@ public:
using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
//constructor
- DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) :
+ DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
+ const SolverType& solver) :
gridView_(gridView), mapper_(mapper), problem_(problem),
- dof_((1 + dim) * gridView.size(0))
+ 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
@@ -57,9 +60,23 @@ public:
A->compress();
Dune::InverseOperatorResult result;
- Dune::Cholmod<Vector> solver;
- solver.setMatrix(*A);
- solver.apply(d, b, result);
+
+ if (solver_ == SolverType::Cholmod)
+ {
+ Dune::Cholmod<Vector> solver;
+ solver.setMatrix(*A);
+ solver.apply(d, b, result);
+ }
+ else if (solver_ == SolverType::UMFPack)
+ {
+ Dune::UMFPack<Matrix> solver(*A);
+ solver.apply(d, b, result);
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception, "ERROR: Specified solver type not "
+ "implemented");
+ }
//write solution to a vtk file
writeVTKoutput();
@@ -113,8 +130,9 @@ public:
protected:
DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
- const int dof) :
- gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
+ 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
@@ -561,6 +579,7 @@ protected:
std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
+ SolverType solver_; //solver type (Cholmod or UMFPack)
private:
//writes the solution to a vtk file
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index be04ec7..70d6a31 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -30,8 +30,8 @@ public:
//constructor
MMDG (const GridView& gridView, const Mapper& mapper,
const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
- const Problem& problem) :
- Base(gridView, mapper, problem,
+ const Problem& problem, const SolverType& solver) :
+ Base(gridView, mapper, problem, solver,
(1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
iGridView_(iGridView), iMapper_(iMapper) {}
@@ -41,9 +41,23 @@ public:
Base::A->compress();
Dune::InverseOperatorResult result;
- Dune::Cholmod<Vector> solver;
- solver.setMatrix(*Base::A);
- solver.apply(Base::d, Base::b, result);
+
+ if (Base::solver_ == SolverType::Cholmod)
+ {
+ Dune::Cholmod<Vector> solver;
+ solver.setMatrix(*Base::A);
+ solver.apply(Base::d, Base::b, result);
+ }
+ else if (Base::solver_ == SolverType::UMFPack)
+ {
+ Dune::UMFPack<Matrix> solver(*Base::A);
+ solver.apply(Base::d, Base::b, result);
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception, "ERROR: Specified solver type not "
+ "implemented");
+ }
//write solution to a vtk file
writeVTKoutput();
diff --git a/dune/mmdg/solvertype.hh b/dune/mmdg/solvertype.hh
new file mode 100644
index 0000000..560113e
--- /dev/null
+++ b/dune/mmdg/solvertype.hh
@@ -0,0 +1,5 @@
+enum SolverType
+{
+ Cholmod,
+ UMFPack
+};
diff --git a/src/dg.cc b/src/dg.cc
index 5304bf3..2739fce 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -25,6 +25,7 @@
#include <dune/mmdg/problems/dgproblem2.hh>
#include <dune/mmdg/problems/dgproblem3.hh>
#include <dune/mmdg/problems/dgproblem4.hh>
+#include <dune/mmdg/problems/dgproblem5.hh>
int main(int argc, char** argv)
{
@@ -101,14 +102,52 @@ int main(int argc, char** argv)
"file parameterDG.ini and make sure that they are specified.");
}
}
+ else if (pt["problem"] == "dg5")
+ {
+ if ( pt.hasKey("dmin") && pt.hasKey("dmax") )
+ {
+ const double dmin = std::stod(pt["dmin"]);
+ const double dmax = std::stod(pt["dmax"]);
+ problem = new DGProblem5<Coordinate>(dmin, dmax);
+
+ std::cout << "NOTE: You will not see convergence for problem dg5 as it "
+ "only fulfills continuity of the normal component of the Darcy "
+ "velocity with respect to the normal of the hyperplane Gamma = "
+ "{ x_1 = 0.5 } but not with respect to the normal of the actual "
+ "interfaces between the bulk and fracture domains. \n\n";
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg4 requires the parameters 'dmin' and 'dmax'. Check the "
+ "file parameterDG.ini and make sure that they are specified.");
+ }
+ }
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'dg1' to 'dg4')");
+ "parameterDG.ini (use 'dg1' to 'dg5')");
+ }
+
+ //determine solver type
+ SolverType solver = SolverType::Cholmod;
+
+ if (pt.hasKey("solver"))
+ {
+ if (pt["solver"] == "UMFPack")
+ {
+ solver = SolverType::UMFPack;
+ }
+ else if (pt["solver"] != "Cholmod")
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid solver type (use 'Cholmod' or 'UMFPack' for the parameter "
+ "'solver' in the file parameterDG.ini)");
+ }
}
- DG dg(gridView, mapper, *problem);
+ DG dg(gridView, mapper, *problem, solver);
dg(mu0);
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 244d92f..aa8d811 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -58,7 +58,7 @@ int main(int argc, char** argv)
Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
const double mu0 = std::stod(pt["mu0"]); //penalty parameter (prefactor)
- const double aperture = std::stod(pt["d"]); //aperture (prefactor)
+ const double dmax = std::stod(pt["dmax"]); //maximum aperture
double xi; //coupling parameter
Problem* problem; //problem to be solved
@@ -68,51 +68,60 @@ int main(int argc, char** argv)
//determine problem type
if (pt["problem"] == "mmdg1")
{
- problem = new MMDGProblem1<Coordinate>(aperture);
+ problem = new MMDGProblem1<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg2")
{
- problem = new MMDGProblem2<Coordinate>(aperture);
+ problem = new MMDGProblem2<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg3")
{
- problem = new MMDGProblem3<Coordinate>(aperture);
+ problem = new MMDGProblem3<Coordinate>(dmax);
gridType = "mmdg3_C";
xi = 1.0;
}
else if (pt["problem"] == "mmdg4")
{
- problem = new MMDGProblem4<Coordinate>(aperture);
+ problem = new MMDGProblem4<Coordinate>(dmax);
gridType = "mmdg2_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg5")
{
- problem = new MMDGProblem5<Coordinate>(aperture);
+ problem = new MMDGProblem5<Coordinate>(dmax);
gridType = "mmdg5_C";
xi = 0.75;
}
else if (pt["problem"] == "mmdg6")
{
xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
- problem = new MMDGProblem6<Coordinate>(aperture, xi);
+ problem = new MMDGProblem6<Coordinate>(dmax, xi);
gridType = "mmdg2_C";
}
else if (pt["problem"] == "mmdg7")
{
- xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
- problem = new MMDGProblem7<Coordinate>(aperture, xi);
- gridType = "mmdg2_C";
+ if (pt.hasKey("dmin"))
+ {
+ problem = new MMDGProblem7<Coordinate>(std::stod(pt["dmin"]), dmax);
+ gridType = "mmdg2_C";
+ xi = 0.75;
+ }
+ else
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg4 requires the parameter 'dmin'. Check the "
+ "file parameterMMDG.ini and make sure that it is specified.");
+ }
}
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'mmdg1' to 'mmdg6')");
+ "parameterDG.ini (use 'mmdg1' to 'mmdg7')");
}
if (pt.hasKey("gridfile")) //non-default grid type
@@ -125,6 +134,23 @@ int main(int argc, char** argv)
xi = std::stod(pt["xi"]);
}
+ //determine solver type
+ SolverType solver = SolverType::Cholmod;
+
+ if (pt.hasKey("solver"))
+ {
+ if (pt["solver"] == "UMFPack")
+ {
+ solver = SolverType::UMFPack;
+ }
+ else if (pt["solver"] != "Cholmod")
+ {
+ DUNE_THROW(Dune::Exception,
+ "Invalid solver type (use 'Cholmod' or 'UMFPack' for the parameter "
+ "'solver' in the file parameterDG.ini)");
+ }
+ }
+
//create a grid from .dgf file
GridFactory gridFactory( "grids/" + gridType + "_" + std::to_string(dim)
+ "d.msh" );
@@ -137,7 +163,7 @@ int main(int argc, char** argv)
const Mapper mapper(gridView, Dune::mcmgElementLayout());
const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
- MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
+ MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem, solver);
mmdg(mu0, xi);
//errors and total run time
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index f6448b0..11dd074 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,5 +1,6 @@
mu0 = 1000
problem = dg4
+dmin = 1e-2
+dmax = 1e-1
+solver = Cholmod
fastEOC = true
-dmin = 1e-1
-dmax = 3e-1
diff --git a/src/parameterMMDG.ini b/src/parameterMMDG.ini
index 257871d..f572d82 100644
--- a/src/parameterMMDG.ini
+++ b/src/parameterMMDG.ini
@@ -1,6 +1,8 @@
mu0 = 1000
-problem = mmdg2
-d = 1e-2
+problem = mmdg6
+dmax = 1e-3
+dmin = 1e-2
#xi = 0.75
gridfile = mmdg2_C
+solver = Cholmod
fastEOC = true
--
GitLab