Commit 4c41292a authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[linear] make linear algebra traits work for residual vectors with state

Apply the logic from the Newton solver.
Test in porousmediumflow/2p2c/implicit/injection.
parent e022365a
......@@ -24,9 +24,24 @@
#ifndef DUMUX_LINEAR_ALGEBRA_TRAITS_HH
#define DUMUX_LINEAR_ALGEBRA_TRAITS_HH
#include <utility>
#include <type_traits>
#include <dune/istl/bvector.hh>
#include <dune/common/std/type_traits.hh>
namespace Dumux {
namespace Detail {
template<class T, std::enable_if_t<Dune::IsNumber<std::decay_t<T>>::value, int> = 0>
constexpr std::size_t blockSize() { return 1; }
template<class T, std::enable_if_t<!Dune::IsNumber<std::decay_t<T>>::value, int> = 0>
constexpr std::size_t blockSize() { return std::decay_t<T>::size(); }
} // end namespace Detail
template<class M, class V>
struct LinearAlgebraTraits
{
......@@ -37,8 +52,14 @@ struct LinearAlgebraTraits
template<class Assembler>
struct LinearAlgebraTraitsFromAssembler
{
private:
using VectorPossiblyWithState = typename Assembler::ResidualType;
using Scalar = std::decay_t<decltype(std::declval<VectorPossiblyWithState>()[0][0])>;
static constexpr auto blockSize = Detail::blockSize<decltype(std::declval<VectorPossiblyWithState>()[0])>();
using BlockType = Dune::FieldVector<Scalar, blockSize>;
public:
using Vector = Dune::BlockVector<BlockType>;
using Matrix = typename Assembler::JacobianMatrix;
using Vector = typename Assembler::ResidualType;
};
} // end namespace Dumux
......
......@@ -26,6 +26,7 @@
#include <memory>
#include <dune/common/version.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/mpicommunication.hh>
......
......@@ -53,6 +53,7 @@
#include <dumux/io/format.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
#include <dumux/linear/algebratraits.hh>
#include <dumux/assembly/partialreassembler.hh>
#include "newtonconvergencewriter.hh"
......
......@@ -32,7 +32,8 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/algebratraits.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/nonlinear/newtonsolver.hh>
......@@ -128,7 +129,8 @@ int main(int argc, char** argv)
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
// the non-linear solver
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment