diff --git a/src/amdis/BiLinearForm.hpp b/src/amdis/BiLinearForm.hpp
index f7d3b8e2a99d13432d8d331d153cf1a9e1fec0f4..8971f493b4c3fe4dde3d215be3e4dbfae0c58541 100644
--- a/src/amdis/BiLinearForm.hpp
+++ b/src/amdis/BiLinearForm.hpp
@@ -43,20 +43,41 @@ namespace AMDiS
using ElementMatrix = FlatMatrix<CoefficientType>;
public:
- /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into
- /// a new shared_ptr.
+ /// Constructor. Stores the row and column basis in a local `shared_ptr` to const
BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
: Super(*rowBasis, *colBasis)
, rowBasis_(rowBasis)
, colBasis_(colBasis)
{
- operators_.init(*rowBasis, *colBasis);
+ operators_.init(*rowBasis_, *colBasis_);
- auto const rowSize = rowBasis->localView().maxSize();
- auto const colSize = colBasis->localView().maxSize();
+ auto const rowSize = rowBasis_->localView().maxSize();
+ auto const colSize = colBasis_->localView().maxSize();
elementMatrix_.resize(rowSize, colSize);
}
+ /// Wraps the passed global bases into (non-destroying) shared_ptr
+ template <class RB_, class CB_,
+ REQUIRES(Concepts::Similar<RB_,RB>),
+ REQUIRES(Concepts::Similar<CB_,CB>)>
+ BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
+ : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
+ {}
+
+ /// Constructor for rowBasis == colBasis
+ template <class RB_ = RB, class CB_ = CB,
+ REQUIRES(std::is_same<RB_,CB_>::value)>
+ explicit BiLinearForm(std::shared_ptr<RB> const& rowBasis)
+ : BiLinearForm(rowBasis, rowBasis)
+ {}
+
+ /// Wraps the passed row-basis into a (non-destroying) shared_ptr
+ template <class RB_,
+ REQUIRES(Concepts::Similar<RB_,RB>)>
+ explicit BiLinearForm(RB_&& rowBasis)
+ : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)))
+ {}
+
std::shared_ptr<RowBasis const> const& rowBasis() const { return rowBasis_; }
std::shared_ptr<ColBasis const> const& colBasis() const { return colBasis_; }
@@ -108,7 +129,6 @@ namespace AMDiS
void assemble(SymmetryStructure symmetry = SymmetryStructure::unknown);
protected:
-
/// Dense matrix to store coefficients during \ref assemble()
ElementMatrix elementMatrix_;
@@ -124,15 +144,26 @@ namespace AMDiS
template <class RB, class CB>
BiLinearForm(RB&&, CB&&)
-> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
+
+ template <class RB>
+ BiLinearForm(RB&&)
+ -> BiLinearForm<Underlying_t<RB>, Underlying_t<RB>>;
#endif
- template <class RB, class CB, class... Args>
+ template <class T = double, class RB, class CB>
auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
{
- using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
+ using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>, T>;
return BLF{FWD(rowBasis), FWD(colBasis)};
}
+ template <class T = double, class RB>
+ auto makeBiLinearForm(RB&& rowBasis)
+ {
+ using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<RB>, T>;
+ return BLF{FWD(rowBasis)};
+ }
+
} // end namespace AMDiS
#include <amdis/BiLinearForm.inc.hpp>
diff --git a/src/amdis/LinearForm.hpp b/src/amdis/LinearForm.hpp
index 78d047184bd40a4d19915b11b441cf5c4db578ae..430c50152ee73c4fcdbd5f69b8f145746bdb65a9 100644
--- a/src/amdis/LinearForm.hpp
+++ b/src/amdis/LinearForm.hpp
@@ -49,6 +49,13 @@ namespace AMDiS
elementVector_.resize(localSize);
}
+ /// Wraps the passed global basis into a (non-destroying) shared_ptr
+ template <class GB_,
+ REQUIRES(Concepts::Similar<GB_,GB>)>
+ explicit LinearForm(GB_&& basis)
+ : LinearForm(Dune::wrap_or_move(FWD(basis)))
+ {}
+
std::shared_ptr<GlobalBasis const> const& basis() const { return basis_; }
/// \brief Associate a local operator with this LinearForm
diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp
index b5a297f9d3741f10a7b1fad091bfe21edad9e911..2094d18c18010f51b7c0678e45c9bda20bd49401 100644
--- a/src/amdis/ProblemStatTraits.hpp
+++ b/src/amdis/ProblemStatTraits.hpp
@@ -6,6 +6,7 @@
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/grid/yaspgrid.hh>
+#include <amdis/AdaptiveGrid.hpp>
#include <amdis/common/Logical.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/ParallelGlobalBasis.hpp>
@@ -70,7 +71,7 @@ namespace AMDiS
template <class HostGrid, class PreBasisCreator, class T = double>
struct DefaultBasisCreator
{
- using Grid = AdaptiveGrid<HostGrid>;
+ using Grid = AdaptiveGrid_t<HostGrid>;
using GridView = typename Grid::LeafGridView;
static auto create(std::string const& name, GridView const& gridView)
diff --git a/src/amdis/common/Algorithm.hpp b/src/amdis/common/Algorithm.hpp
index dabbe0331cc384a9584e1e306d06faae8147871b..7805a66f42f5ba1c8b02c11566069496e59f46d2 100644
--- a/src/amdis/common/Algorithm.hpp
+++ b/src/amdis/common/Algorithm.hpp
@@ -1,36 +1,67 @@
#pragma once
#include <algorithm>
+#include <functional>
+#include <utility>
namespace AMDiS
{
- template <class InputIter, class T, class Func>
- void split(InputIter first, InputIter end, T const& t, Func f)
+ /// \brief Split a sequence `[first,last)` by the separators `sep` and pass the tokens as
+ /// begin-end iterator pair to the provided functor `f = void(InputIterator, InputIterator)`
+ template <class InputIter, class Tp, class BinaryFunc>
+ void split(InputIter first, InputIter last, Tp sep, BinaryFunc f)
{
- if (first == end)
+ if (first == last)
return;
while (true) {
- InputIter found = std::find(first, end, t);
+ InputIter found = std::find(first, last, sep);
f(first, found);
- if (found == end)
+ if (found == last)
break;
first = ++found;
}
}
- template <class InputIter, class SeparaterIter, class Func>
- void split(InputIter first, InputIter end, SeparaterIter s_first, SeparaterIter s_end, Func f)
+ /// \brief Split a sequence `[first,last)` by any of the separators `[s_first, s_last)` and pass the tokens as
+ /// begin-end iterator pair to the provided functor `f = void(InputIterator, InputIterator)`
+ template <class InputIter, class SeparaterIter, class BinaryFunc>
+ void split(InputIter first, InputIter last, SeparaterIter s_first, SeparaterIter s_last, BinaryFunc f)
{
- if (first == end)
+ if (first == last)
return;
while (true) {
- InputIter found = std::find_first_of(first, end, s_first, s_end);
+ InputIter found = std::find_first_of(first, last, s_first, s_last);
f(first, found);
- if (found == end)
+ if (found == last)
break;
first = ++found;
}
}
+
+ /// \brief Output the cumulative sum of one range to a second range
+ // NOTE: backport of std::exclusive_scan from c++17
+ template <class InputIter, class OutputIter, class Tp, class BinaryOperation>
+ OutputIter exclusive_scan(InputIter first, InputIter last,
+ OutputIter result, Tp init, BinaryOperation binary_op)
+ {
+ while (first != last) {
+ auto v = init;
+ init = binary_op(init, *first);
+ ++first;
+ *result++ = std::move(v);
+ }
+ return result;
+ }
+
+ /// \brief Output the cumulative sum of one range to a second range
+ // NOTE: backport of std::exclusive_scan from c++17
+ template <class InputIter, class OutputIter, class Tp>
+ inline OutputIter exclusive_scan(InputIter first, InputIter last,
+ OutputIter result, Tp init)
+ {
+ return AMDiS::exclusive_scan(first, last, result, std::move(init), std::plus<>());
+ }
+
}
\ No newline at end of file
diff --git a/src/amdis/linearalgebra/LinearSolver.hpp b/src/amdis/linearalgebra/LinearSolver.hpp
index 3facaa4a2fa6220fbed0c886c4dc14ea213b33e1..1eb051b8b5dee1e65e42dfa775bb73cfbbe60acb 100644
--- a/src/amdis/linearalgebra/LinearSolver.hpp
+++ b/src/amdis/linearalgebra/LinearSolver.hpp
@@ -44,13 +44,17 @@ namespace AMDiS
public:
/// Constructor
- explicit LinearSolver(std::string prefix)
- : runner_(prefix)
+ template <class... Args>
+ explicit LinearSolver(std::string prefix, Args&&... args)
+ : runner_(prefix, FWD(args)...)
{}
+ Runner& runner() { return runner_; }
+ Runner const& runner() const { return runner_; }
+
private:
/// Implements \ref LinearSolverInterface::solveSystemImpl()
- void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm& comm, SolverInfo& solverInfo) override
+ void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo) override
{
Dune::Timer t;
if (solverInfo.doCreateMatrixData()) {
diff --git a/src/amdis/linearalgebra/LinearSolverInterface.hpp b/src/amdis/linearalgebra/LinearSolverInterface.hpp
index bbaee462ae6ecfefc0f8772c69d6f2bfafda56b6..14e255206beea3ec20b9852909ac5222eb8a1e59 100644
--- a/src/amdis/linearalgebra/LinearSolverInterface.hpp
+++ b/src/amdis/linearalgebra/LinearSolverInterface.hpp
@@ -38,14 +38,14 @@ namespace AMDiS
* \p x A [block-]vector for the unknown components.
* \p b A [block-]vector for the right-hand side of the linear system.
**/
- void solve(Mat const& A, Vec& x, Vec const& b, Comm& comm, SolverInfo& solverInfo)
+ void solve(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo)
{
solveImpl(A, x, b, comm, solverInfo);
}
private:
/// main methods that all solvers must implement
- virtual void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm& comm, SolverInfo& solverInfo) = 0;
+ virtual void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo) = 0;
};
} // end namespace AMDiS
diff --git a/src/amdis/linearalgebra/RunnerInterface.hpp b/src/amdis/linearalgebra/RunnerInterface.hpp
index 4f46d3e456c0d23b32b4dc451ba72cca81b7e5a9..7fcf48d7dea87f105edb018006ea0d31af880201 100644
--- a/src/amdis/linearalgebra/RunnerInterface.hpp
+++ b/src/amdis/linearalgebra/RunnerInterface.hpp
@@ -20,7 +20,7 @@ namespace AMDiS
virtual ~RunnerInterface() = default;
/// Is called at the beginning of a solution procedure
- virtual void init(Mat const& A, Comm& comm) = 0;
+ virtual void init(Mat const& A, Comm const& comm) = 0;
/// Is called at the end of a solution procedure
virtual void exit() = 0;
diff --git a/src/amdis/linearalgebra/eigen/DirectRunner.hpp b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
index 61f74dade1152b2389677f15ab38146b208771f7..47649fc6c2faad44b9207d44cf35fbccc1042de3 100644
--- a/src/amdis/linearalgebra/eigen/DirectRunner.hpp
+++ b/src/amdis/linearalgebra/eigen/DirectRunner.hpp
@@ -35,7 +35,7 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::init()
- void init(Mat const& A, Comm& comm) override
+ void init(Mat const& A, Comm const& comm) override
{
DUNE_UNUSED_PARAMETER(comm);
@@ -49,9 +49,11 @@ namespace AMDiS
"Error in solver.compute(matrix)");
}
-
/// Implements \ref RunnerInterface::exit()
- void exit() override {}
+ void exit() override
+ {
+ initialized_ = false;
+ }
/// Implements \ref RunnerInterface::solve()
int solve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
diff --git a/src/amdis/linearalgebra/eigen/IterativeRunner.hpp b/src/amdis/linearalgebra/eigen/IterativeRunner.hpp
index 45541907b9d7a8a89d0ac13c0ce0b0e6075c5ae6..ad66327324adca85af122dd5360cd533fb68c7a8 100644
--- a/src/amdis/linearalgebra/eigen/IterativeRunner.hpp
+++ b/src/amdis/linearalgebra/eigen/IterativeRunner.hpp
@@ -30,7 +30,7 @@ namespace AMDiS
/// Implements \ref RunnerInterface::init()
- void init(Mat const& A, Comm& comm) override
+ void init(Mat const& A, Comm const& comm) override
{
DUNE_UNUSED_PARAMETER(comm);
@@ -45,7 +45,10 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::exit()
- void exit() override {}
+ void exit() override
+ {
+ initialized_ = false;
+ }
/// Implements \ref RunnerInterface::solve()
int solve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
diff --git a/src/amdis/linearalgebra/istl/ISTLRunner.hpp b/src/amdis/linearalgebra/istl/ISTLRunner.hpp
index 7eaf846a41f4feba63bedf24548ad6e2e37975b5..5c688535a5c43de46ef901272407d99a5da23dcb 100644
--- a/src/amdis/linearalgebra/istl/ISTLRunner.hpp
+++ b/src/amdis/linearalgebra/istl/ISTLRunner.hpp
@@ -32,7 +32,7 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::init()
- void init(Mat const& mat, typename Traits::Comm& comm) override
+ void init(Mat const& mat, typename Traits::Comm const& comm) override
{
solver_ = solverCreator_->create(mat, comm);
}
diff --git a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
index 3087971bc11c93c533a3999d3ede8554a624424d..3ec3ee6aae7572ca4336c2351a9e7e5d185d8aa3 100644
--- a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
@@ -41,13 +41,14 @@ namespace AMDiS
Parameters::get(prefix_ + "->relative tolerance", rTol_);
Parameters::get(prefix_ + "->max iteration", maxIter_);
Parameters::get(prefix_ + "->print cycle", printCycle_);
+
+ initPrecon(prefix_);
}
/// Implementation of \ref RunnerInterface::init()
- void init(Matrix const& A, Comm& comm) override
+ void init(Matrix const& A, Comm const& comm) override
{
DUNE_UNUSED_PARAMETER(comm);
- initPrecon(prefix_);
P_->init(A);
}
diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
index 4703b6ca4670888e87d5998c1705635657955d20..4c9065f34df1191ce98fa78c694b4226a7b4fa95 100644
--- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
@@ -45,7 +45,7 @@ namespace AMDiS
}
/// Implementation of \ref RunnerInterface::init()
- void init(Matrix const& matrix, Comm&) override
+ void init(Matrix const& matrix, Comm const&) override
{
try {
if (bool(solver_) && storeSymbolic_)
diff --git a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
index e72a33fbe9b5d0927cd76bc56f4e58e49fd748eb..e00b5874ee5f98bef06baff5995503b3de9a3f8b 100644
--- a/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
+++ b/src/amdis/linearalgebra/petsc/MatrixBackend.hpp
@@ -13,11 +13,13 @@
namespace AMDiS
{
+ template <class> class Constraints;
+
/// \brief The basic container that stores a base matrix
template <class DofMap>
class PetscMatrix
{
- friend struct Constraints<PetscMatrix<DofMap>>;
+ template <class> friend class Constraints;
public:
/// The matrix type of the underlying base matrix
diff --git a/src/amdis/linearalgebra/petsc/PetscRunner.hpp b/src/amdis/linearalgebra/petsc/PetscRunner.hpp
index 55f6ea893a92810f501bd26409b47dc2139733ce..7a233c5cd698365b9aa9d6a291f25299d6e5c1ad 100644
--- a/src/amdis/linearalgebra/petsc/PetscRunner.hpp
+++ b/src/amdis/linearalgebra/petsc/PetscRunner.hpp
@@ -68,7 +68,7 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::init()
- void init(typename Traits::Mat const& mat, typename Traits::Comm& comm) override
+ void init(typename Traits::Mat const& mat, typename Traits::Comm const& comm) override
{
exit();
#if HAVE_MPI
@@ -110,23 +110,26 @@ namespace AMDiS
protected:
// initialize the KSP solver from the initfile
- void initKSP(KSP ksp, std::string prefix) const
+ virtual void initKSP(KSP ksp, std::string prefix) const
{
// see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPType.html
- auto kspTypeStr = Parameters::get<std::string>(prefix);
- if (kspTypeStr) {
- if (kspTypeStr.value() == "direct")
- initDirectSolver(ksp, prefix);
- else if (kspTypeStr.value() != "default") {
- KSPSetType(ksp, kspTypeStr.value().c_str());
-
- // initialize some KSP specific parameters
- initKSPParameters(ksp, kspTypeStr.value().c_str(), prefix);
-
- // set initial guess to nonzero only for non-preonly ksp type
- if (kspTypeStr.value() != "preonly")
- KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
- }
+ auto kspType = Parameters::get<std::string>(prefix + "->ksp");
+ std::string kspTypeStr = kspType.value_or("default");
+
+ if (!kspType)
+ Parameters::get(prefix, kspTypeStr);
+
+ if (kspTypeStr == "direct")
+ initDirectSolver(ksp, prefix);
+ else if (kspTypeStr != "default") {
+ KSPSetType(ksp, kspTypeStr.c_str());
+
+ // initialize some KSP specific parameters
+ initKSPParameters(ksp, kspTypeStr.c_str(), prefix);
+
+ // set initial guess to nonzero only for non-preonly ksp type
+ if (kspTypeStr != "preonly")
+ KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
}
// set a KSPMonitor if info > 0
@@ -167,7 +170,7 @@ namespace AMDiS
// initialize a direct solver from the initfile
- void initDirectSolver(KSP ksp, std::string prefix) const
+ virtual void initDirectSolver(KSP ksp, std::string prefix) const
{
KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
KSPSetType(ksp, KSPRICHARDSON);
@@ -192,7 +195,7 @@ namespace AMDiS
// initialize the preconditioner pc from the initfile
- void initPC(PC pc, std::string prefix) const
+ virtual void initPC(PC pc, std::string prefix) const
{
// see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCType.html
auto pcType = Parameters::get<std::string>(prefix);
@@ -234,7 +237,7 @@ namespace AMDiS
// provide initfile parameters for some PETSc KSP parameters
- void initKSPParameters(KSP ksp, char const* ksptype, std::string prefix) const
+ virtual void initKSPParameters(KSP ksp, char const* ksptype, std::string prefix) const
{
// parameters for the Richardson solver
if (std::strcmp(ksptype, KSPRICHARDSON) == 0)
@@ -288,7 +291,7 @@ namespace AMDiS
}
}
- private:
+ protected:
std::string prefix_;
int info_ = 0;