OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-07 04:53:03 -05:00
Code Issues Packages Projects Releases Wiki Activity

Add DILU preconditioner

Browse Source
This commit is contained in:
jakobtorben 2023-09-27 15:47:05 +02:00
parent d98de155de
commit ab0ca76194
8 changed files with 1068 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -251,6 +251,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_convergenceoutputconfiguration.cpp
tests/test_convergencereport.cpp
tests/test_deferredlogger.cpp
tests/test_dilu.cpp
tests/test_eclinterregflows.cpp
tests/test_equil.cc
tests/test_extractMatrix.cpp
@ -484,8 +485,10 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/bda/rocsparseWellContributions.hpp
opm/simulators/linalg/bda/WellContributions.hpp
opm/simulators/linalg/amgcpr.hh
opm/simulators/linalg/DILU.hpp
opm/simulators/linalg/twolevelmethodcpr.hh
opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp
opm/simulators/linalg/ExtraSmoothers.hpp
opm/simulators/linalg/FlexibleSolver.hpp
opm/simulators/linalg/FlexibleSolver_impl.hpp
opm/simulators/linalg/FlowLinearSolverParameters.hpp

192
opm/simulators/linalg/DILU.hpp Normal file
View File

@ -0,0 +1,192 @@
/*
Copyright 2022-2023 SINTEF AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_DILU_HEADER_INCLUDED
#define OPM_DILU_HEADER_INCLUDED
#include <config.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
#include <dune/common/unused.hh>
#include <dune/istl/bcrsmatrix.hh>
namespace Dune
{
/*! \brief The sequential DILU preconditioner.
\tparam M The matrix type to operate on
\tparam X Type of the update
\tparam Y Type of the defect
*/
template <class M, class X, class Y>
class SeqDilu : public PreconditionerWithUpdate<X, Y>
{
public:
//! \brief The matrix type the preconditioner is for.
using matrix_type = M;
//! \brief The domain type of the preconditioner.
using domain_type = X;
//! \brief The range type of the preconditioner.
using range_type = Y;
//! \brief The field type of the preconditioner.
using field_type = typename X::field_type;
//! \brief scalar type underlying the field_type
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
*/
SeqDilu(const M& A)
: A_(A)
{
Dinv_.resize(A_.N());
// the Dinv matrix must be initialised
update();
}
/*!
\brief Update the preconditioner.
\copydoc Preconditioner::update()
*/
virtual void update() override
{
OPM_TIMEBLOCK(update);
auto endi = A_.end();
for ( auto row = A_.begin(); row != endi; ++row) {
const auto row_i = row.index();
Dinv_[row_i] = A_[row_i][row_i];
}
for ( auto row = A_.begin(); row != endi; ++row)
{
const auto row_i = row.index();
auto Dinv_temp = Dinv_[row_i];
for (auto a_ij = row->begin(); a_ij.index() < row_i; ++a_ij)
{
const auto col_j = a_ij.index();
const auto a_ji = A_[col_j].find(row_i);
// if A[i, j] != 0 and A[j, i] != 0
if (a_ji != A_[col_j].end()) {
// Dinv_temp -= A[i, j] * d[j] * A[j, i]
Dinv_temp -= (*a_ij) * Dune::FieldMatrix(Dinv_[col_j]) * (*a_ji);
}
}
Dinv_temp.invert();
Dinv_[row_i] = Dinv_temp;
}
}
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre(X& v, Y& d) override
{
DUNE_UNUSED_PARAMETER(v);
DUNE_UNUSED_PARAMETER(d);
}
/*!
\brief Apply the preconditioner.
\copydoc Preconditioner::apply(X&,const Y&)
*/
virtual void apply(X& v, const Y& d) override
{
// M = (D + L_A) D^-1 (D + U_A) (a LU decomposition of M)
// where L_A and U_A are the strictly lower and upper parts of A and M has the properties:
// diag(A) = diag(M)
// Working with defect d = b - Ax and update v = x_{n+1} - x_n
// solving the product M^-1(d) using upper and lower triangular solve
// v = M^{-1}*d = (D + U_A)^{-1} D (D + L_A)^{-1} * d
OPM_TIMEBLOCK(apply);
using Xblock = typename X::block_type;
using Yblock = typename Y::block_type;
// lower triangular solve: (D + L_A) y = d
auto endi = A_.end();
for (auto row = A_.begin(); row != endi; ++row)
{
const auto row_i = row.index();
Yblock rhs = d[row_i];
for (auto a_ij = (*row).begin(); a_ij.index() < row_i; ++a_ij) {
// if A[i][j] != 0
// rhs -= A[i][j]* y[j], where v_j stores y_j
const auto col_j = a_ij.index();
a_ij->mmv(v[col_j], rhs);
}
// y_i = Dinv_i * rhs
// storing y_i in v_i
Dinv_[row_i].mv(rhs, v[row_i]); // (D + L_A)_ii = D_i
}
// upper triangular solve: (D + U_A) v = Dy
auto rendi = A_.beforeBegin();
for (auto row = A_.beforeEnd(); row != rendi; --row)
{
const auto row_i = row.index();
// rhs = 0
Xblock rhs(0.0);
for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij) {
// if A[i][j] != 0
// rhs += A[i][j]*v[j]
const auto col_j = a_ij.index();
a_ij->umv(v[col_j], rhs);
}
// calculate update v = M^-1*d
// v_i = y_i - Dinv_i*rhs
// before update v_i is y_i
Dinv_[row_i].mmv(rhs, v[row_i]);
}
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post(X& x) override
{
DUNE_UNUSED_PARAMETER(x);
}
std::vector<typename M::block_type> getDiagonal()
{
return Dinv_;
}
//! Category of the preconditioner (see SolverCategory::Category)
virtual SolverCategory::Category category() const override
{
return SolverCategory::sequential;
}
private:
//! \brief The matrix we operate on.
const M& A_;
//! \brief The inverse of the diagnal matrix
std::vector<typename M::block_type> Dinv_;
};
} // namespace Dune
#endif

38
opm/simulators/linalg/ExtraSmoothers.hpp Normal file
View File

@ -0,0 +1,38 @@
#ifndef OPM_EXTRASMOOTHERS_HPP
#define OPM_EXTRASMOOTHERS_HPP
#include "DILU.hpp"
namespace Dune
{
template <class M, class X, class Y>
class SeqDilu;
namespace Amg
{
/**
* @brief Policy for the construction of the SeqDilu smoother
*/
template <class M, class X, class Y>
struct ConstructionTraits<SeqDilu<M, X, Y>> {
using Arguments = DefaultConstructionArgs<SeqDilu<M, X, Y>>;
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
static inline std::shared_ptr<SeqDilu<M, X, Y>> construct(Arguments& args) {
return std::make_shared<SeqDilu<M, X, Y>>(args.getMatrix());
}
#else
static inline SeqDilu<M, X, Y>* construct(Arguments& args) {
return new SeqDilu<M, X, Y>(args.getMatrix());
}
static void deconstruct(SeqDilu<M, X, Y>* dilu) {
delete dilu;
}
#endif
};
} // namespace Amg
} // namespace Dune
#endif // OPM_EXTRASMOOTHERS_HPP

2
opm/simulators/linalg/FlowLinearSolverParameters.hpp
View File

@ -307,7 +307,7 @@ namespace Opm
EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver");
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types");
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes, amg or hybrid (experimental). Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolver, "Configuration of solver. Valid options are: ilu0 (default), dilu, cprw, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes, amg or hybrid (experimental). Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverPrintJsonDefinition, "Write the JSON definition of the linear solver setup to the DBG file.");
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate, 4: recreated every CprReuseInterval");
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseInterval, "Reuse preconditioner interval. Used when CprReuseSetup is set to 4, then the preconditioner will be fully recreated instead of reused every N linear solve, where N is this parameter.");

28
opm/simulators/linalg/PreconditionerFactory_impl.hpp
View File

@ -34,7 +34,10 @@
#include <opm/simulators/linalg/PressureTransferPolicy.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/WellOperators.hpp>
#include <opm/simulators/linalg/DILU.hpp>
#include <opm/simulators/linalg/ExtraSmoothers.hpp>
#include <dune/common/unused.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/paamg/amg.hh>
@ -91,6 +94,7 @@ struct AMGSmootherArgsHelper<Opm::ParallelOverlappingILU0<M,V,V,C>>
}
};
template <class Operator, class Comm, class Matrix, class Vector>
typename AMGHelper<Operator, Comm, Matrix, Vector>::Criterion
AMGHelper<Operator,Comm,Matrix,Vector>::criterion(const PropertyTree& prm)
@ -158,6 +162,10 @@ struct StandardPreconditioners
F::addCreator("ILUn", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
});
F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
DUNE_UNUSED_PARAMETER(prm);
return wrapBlockPreconditioner<SeqDilu<M, V, V>>(comm, op.getmat());
});
F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&,
std::size_t, const C& comm) {
const int n = prm.get<int>("repeats", 1);
@ -186,12 +194,23 @@ struct StandardPreconditioners
// with the AMG hierarchy construction.
if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
F::addCreator("amg", [](const O& op, const P& prm, const std::function<V()>&, std::size_t, const C& comm) {
using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<V, V>>;
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
auto crit = AMGHelper<O,C,M,V>::criterion(prm);
auto sargs = AMGSmootherArgsHelper<Smoother>::args(prm);
return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
return prec;
}
else if (smoother == "DILU") {
using SeqSmoother = Dune::SeqDilu<M, V, V>;
using Smoother = Dune::BlockPreconditioner<V, V, C, SeqSmoother>;
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
SmootherArgs sargs;
auto crit = AMGHelper<O,C,M,V>::criterion(prm);
PrecPtr prec = std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
return prec;
} else {
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " + smoother + ".");
}
@ -329,6 +348,10 @@ struct StandardPreconditioners<Operator,Dune::Amg::SequentialInformation>
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
F::addCreator("DILU", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
DUNE_UNUSED_PARAMETER(prm);
return std::make_shared<SeqDilu<M, V, V>>(op.getmat());
});
F::addCreator("Jac", [](const O& op, const P& prm, const std::function<V()>&, std::size_t) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
@ -361,6 +384,9 @@ struct StandardPreconditioners<Operator,Dune::Amg::SequentialInformation>
} else if (smoother == "Jac") {
using Smoother = SeqJac<M, V, V>;
return AMGHelper<O,C,M,V>::template makeAmgPreconditioner<Smoother>(op, prm);
} else if (smoother == "DILU") {
using Smoother = SeqDilu<M, V, V>;
return AMGHelper<O,C,M,V>::template makeAmgPreconditioner<Smoother>(op, prm);
} else if (smoother == "SOR") {
using Smoother = SeqSOR<M, V, V>;
return AMGHelper<O,C,M,V>::template makeAmgPreconditioner<Smoother>(op, prm);

19
opm/simulators/linalg/setupPropertyTree.cpp
View File

@ -96,6 +96,10 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
return setupILU(conf, p);
}
if (conf == "dilu") {
return setupDILU(conf, p);
}
if (conf == "umfpack") {
return setupUMFPack(conf, p);
}
@ -111,7 +115,7 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters
// No valid configuration option found.
OPM_THROW(std::invalid_argument,
conf + " is not a valid setting for --linear-solver-configuration."
" Please use ilu0, cpr, cpr_trueimpes, cpr_quasiimpes or isai");
" Please use ilu0, dilu, cpr, cpr_trueimpes, cpr_quasiimpes or isai");
}
std::string getSolverString(const FlowLinearSolverParameters& p)
@ -266,6 +270,19 @@ setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParamet
return prm;
}
PropertyTree
setupDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
{
using namespace std::string_literals;
PropertyTree prm;
prm.put("tol", p.linear_solver_reduction_);
prm.put("maxiter", p.linear_solver_maxiter_);
prm.put("verbosity", p.linear_solver_verbosity_);
prm.put("solver", getSolverString(p));
prm.put("preconditioner.type", "DILU"s);
return prm;
}
PropertyTree
setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)

1
opm/simulators/linalg/setupPropertyTree.hpp
View File

@ -37,6 +37,7 @@ PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters
PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p);
PropertyTree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
PropertyTree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);
PropertyTree setupDILU(const std::string& conf, const FlowLinearSolverParameters& p);
PropertyTree setupUMFPack(const std::string& conf, const FlowLinearSolverParameters& p);
} // namespace Opm

788
tests/test_dilu.cpp Normal file
View File

@ -0,0 +1,788 @@
/*
Copyright 2022-2023 SINTEF AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#define BOOST_TEST_MODULE TestSeqDILU
#include <config.h>
#include <opm/simulators/linalg/DILU.hpp>
#include <boost/mpl/list.hpp>
#include <boost/test/unit_test.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/preconditioners.hh>
using NumericTypes = boost::mpl::list<double, float>;
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2NoZeros, T, NumericTypes)
{
/*
Tests that the dilu decomposition mathces the expected result
for a 2x2 matrix with no zero blocks, with block size 2x2.
A
| | 3 1| | 1 0| |
| | 2 1| | 0 1| |
| |
| | 2 0| |-1 0| |
| | 0 2| | 0 -1| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 4;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(0);
row.insert(1);
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[0][1][0][0]=1.0;
A[0][1][1][1]=1.0;
A[1][0][0][0]=2.0;
A[1][0][1][1]=2.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_00_inv U_01
auto D_11 = A[1][1] - A[1][0]*D_00_inv*A[0][1];
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
auto Dinv = seqdilu.getDiagonal();
// diagonal stores inverse
auto D_00_dilu = Dinv[0];
D_00_dilu.invert();
auto D_11_dilu = Dinv[1];
D_11_dilu.invert();
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(D_00_dilu[i][j], D_00[i][j], 1e-7);
BOOST_CHECK_CLOSE(D_11_dilu[i][j], D_11[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect2x2, T, NumericTypes)
{
/*
Tests that the dilu decomposition mathces the expected result
for a 2x2 matrix, with block size 2x2.
A
| | 3 1| | 1 0| |
| | 2 1| | 0 1| |
| |
| | 0 0| |-1 0| |
| | 0 0| | 0 -1| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 3;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(row.index());
if (row.index() == 0) {
row.insert(row.index() + 1);
}
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[0][1][0][0]=1.0;
A[0][1][1][1]=1.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_00_inv U_01 = A_11
auto D_11 = A[1][1];
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
auto Dinv = seqdilu.getDiagonal();
// diagonal stores inverse
auto D_00_dilu = Dinv[0];
D_00_dilu.invert();
auto D_11_dilu = Dinv[1];
D_11_dilu.invert();
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(D_00_dilu[i][j], D_00[i][j], 1e-7);
BOOST_CHECK_CLOSE(D_11_dilu[i][j], D_11[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrectNoZeros, T, NumericTypes)
{
/*
Tests that applying the dilu preconditioner mathces the expected result
for a 2x2 matrix with no zero blocks, with block size 2x2.
A x = b
| | 3 1| | 1 0| | | |1| | | |2| |
| | 2 1| | 0 1| | | |2| | | |1| |
| | x | | = | |
| | 2 0| |-1 0| | | |1| | | |3| |
| | 0 2| | 0 -1| | | |1| | | |4| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 4;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(0);
row.insert(1);
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[0][1][0][0]=1.0;
A[0][1][1][1]=1.0;
A[1][0][0][0]=2.0;
A[1][0][1][1]=2.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
Vector x(2);
x[0][0] = 1.0;
x[0][1] = 2.0;
x[1][0] = 1.0;
x[1][1] = 1.0;
Vector b(2);
b[0][0] = 2.0;
b[0][1] = 1.0;
b[1][0] = 3.0;
b[1][1] = 4.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11= A_11 - L_10 D_00_inv U_01
auto D_11 = A[1][1] - A[1][0]*D_00_inv*A[0][1];
auto D_11_inv = D_11;
D_11_inv.invert();
// define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A)
// lower triangular solve: (E + L) y = b - Ax
// y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
Vector y(2);
auto rhs = b[0];
A[0][0].mmv(x[0], rhs);
A[0][1].mmv(x[1], rhs);
D_00_inv.mv(rhs, y[0]);
// y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1) - A_10*y_0)
rhs = b[1];
A[1][0].mmv(x[0], rhs);
A[1][1].mmv(x[1], rhs);
A[1][0].mmv(y[0], rhs);
D_11_inv.mv(rhs, y[1]);
// upper triangular solve: (E + U) z = Ey
// z_1 = y_1
Vector z(2);
z[1] = y[1];
// z_0 = y_0 - D_00_inv*A_01*z_1
z[0] = y[0];
auto temp = D_00_inv*A[0][1];
temp.mmv(z[1], z[0]);
// x_k+1 = x_k + z
Vector new_x = x;
new_x += z;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
seqdilu.apply(x, b);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(x[i][j], new_x[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect1, T, NumericTypes)
{
/*
Tests that applying the dilu preconditioner mathces the expected result
for a 2x2 matrix, with block size 2x2.
A x = b
| | 3 1| | 1 0| | | |1| | | |2| |
| | 2 1| | 0 1| | | |2| | | |1| |
| | x | | = | |
| | 0 0| |-1 0| | | |1| | | |3| |
| | 0 0| | 0 -1| | | |1| | | |4| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 3;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(row.index());
if (row.index() == 0) {
row.insert(row.index() + 1);
}
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[0][1][0][0]=1.0;
A[0][1][1][1]=1.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
Vector x(2);
x[0][0] = 1.0;
x[0][1] = 2.0;
x[1][0] = 1.0;
x[1][1] = 1.0;
Vector b(2);
b[0][0] = 2.0;
b[0][1] = 1.0;
b[1][0] = 3.0;
b[1][1] = 4.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_0_inv U_01 = A_11
auto D_11 = A[1][1];
auto D_11_inv = D_11;
D_11_inv.invert();
// define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A)
// lower triangular solve: (E + L) y = b - Ax
// y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
Vector y(2);
auto rhs = b[0];
A[0][0].mmv(x[0], rhs);
A[0][1].mmv(x[1], rhs);
D_00_inv.mv(rhs, y[0]);
// y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1) - A_10*y_0) = D_11_inv*(b_1 - A_11*x_1)
rhs = b[1];
A[1][1].mmv(x[1], rhs);
D_11_inv.mv(rhs, y[1]);
// upper triangular solve: (E + U) z = Ey
// z_1 = y_1
Vector z(2);
z[1] = y[1];
// z_0 = y_0 - D_00_inv*A_01*z_1
z[0] = y[0];
auto temp = D_00_inv*A[0][1];
temp.mmv(z[1], z[0]);
// x_k+1 = x_k + z
Vector new_x = x;
new_x += z;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
seqdilu.apply(x, b);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(x[i][j], new_x[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect2, T, NumericTypes)
{
/*
Tests that applying the dilu preconditioner mathces the expected result
for a 2x2 matrix, with block size 2x2.
A x = b
| | 3 1| | 0 0| | | |1| | | |2| |
| | 2 1| | 0 0| | | |2| | | |1| |
| | x | | = | |
| | 2 0| |-1 0| | | |1| | | |3| |
| | 0 2| | 0 -1| | | |1| | | |4| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 3;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(row.index());
if (row.index() == 1) {
row.insert(row.index() - 1);
}
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[1][1][0][0]=2.0;
A[1][1][1][1]=2.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
Vector x(2);
x[0][0] = 1.0;
x[0][1] = 2.0;
x[1][0] = 1.0;
x[1][1] = 1.0;
Vector b(2);
b[0][0] = 2.0;
b[0][1] = 1.0;
b[1][0] = 3.0;
b[1][1] = 4.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_0_inv U_01 = A_11
auto D_11 = A[1][1];
auto D_11_inv = D_11;
D_11_inv.invert();
// define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A)
// lower triangular solve: (E + L) y = b - Ax
// y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)] = D_00_inv*[b_0 - A_00*x_0]
Vector y(2);
auto rhs = b[0];
A[0][0].mmv(x[0], rhs);
D_00_inv.mv(rhs, y[0]);
// y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1) - A_10*y_0)
rhs = b[1];
A[1][1].mmv(x[1], rhs);
D_11_inv.mv(rhs, y[1]);
// upper triangular solve: (E + U) z = Ey
// z_1 = y_1
Vector z(2);
z[1] = y[1];
// z_0 = y_0 - D_00_inv*A_01*z_1 = y_0
z[0] = y[0];
// x_k+1 = x_k + z
Vector new_x = x;
new_x += z;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
seqdilu.apply(x, b);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(x[i][j], new_x[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUDiagIsCorrect3x3, T, NumericTypes)
{
/*
Tests that the dilu decomposition mathces the expected result
for a 3x3 matrix, with block size 3x3.
A
| | 3 1 2| | 0 0 0| | 0 0 0| |
| | 2 3 1| | 0 0 0| | 0 0 0| |
| | 2 1 0| | 0 0 0| | 0 0 0| |
| |
| | 0 0 0| | 1 0 1| | 1 0 2| |
| | 0 0 0| | 4 1 0| | 0 1 1| |
| | 0 0 0| | 3 1 3| | 0 1 3| |
| |
| | 0 0 0| | 1 0 2| | 1 3 2| |
| | 0 0 0| | 0 1 4| | 2 1 3| |
| | 0 0 0| | 5 1 1| | 3 1 2| |
*/
const int N = 3;
const int bz = 3;
const int nonZeroes = 5;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
if (row.index() == 0) {
row.insert(row.index());
}
else if (row.index() == 1) {
row.insert(row.index());
row.insert(row.index() + 1);
}
else if (row.index() == 2) {
row.insert(row.index() - 1);
row.insert(row.index());
}
}
A[0][0][0][0]=3.0; A[1][1][0][0]=1.0; A[1][2][0][0]=1.0;
A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0;
A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0;
A[0][0][1][0]=2.0; A[1][1][1][0]=4.0; A[1][2][1][0]=0.0;
A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0;
A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0;
A[0][0][2][0]=2.0; A[1][1][2][0]=3.0; A[1][2][2][0]=0.0;
A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0;
A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.0;
A[2][1][0][0]=1.0; A[2][2][0][0]=1.0;
A[2][1][0][1]=0.0; A[2][2][0][1]=3.0;
A[2][1][0][2]=2.0; A[2][2][0][2]=2.0;
A[2][1][1][0]=0.0; A[2][2][1][0]=2.0;
A[2][1][1][1]=1.0; A[2][2][1][1]=1.0;
A[2][1][1][2]=4.0; A[2][2][1][2]=3.0;
A[2][1][2][0]=5.0; A[2][2][2][0]=3.0;
A[2][1][2][1]=1.0; A[2][2][2][1]=1.0;
A[2][1][2][2]=1.0; A[2][2][2][2]=2.0;
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_00_inv U_01 = A_11
auto D_11 = A[1][1];
auto D_11_inv = D_11;
D_11_inv.invert();
// D_22 = A_22 - A_20 D_00_inv A_02 - A_21 D_11_inv A_12 = A_22 - A_21 D_11_inv A_12
auto D_22 = A[2][2] - A[2][1]*D_11_inv*A[1][2];
auto D_22_inv = D_22;
D_22_inv.invert();
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
auto Dinv = seqdilu.getDiagonal();
// diagonal stores inverse
auto D_00_dilu = Dinv[0];
D_00_dilu.invert();
auto D_11_dilu = Dinv[1];
D_11_dilu.invert();
auto D_22_dilu = Dinv[2];
D_22_dilu.invert();
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
BOOST_CHECK_CLOSE(D_00_dilu[i][j], D_00[i][j], 1e-7);
BOOST_CHECK_CLOSE(D_11_dilu[i][j], D_11[i][j], 1e-7);
BOOST_CHECK_CLOSE(D_22_dilu[i][j], D_22[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsCorrect3, T, NumericTypes)
{
/*
Tests that applying the dilu preconditioner mathces the expected result
for a 3x3 matrix, with block size 3x3.
A x = b
| | 3 1 2| | 0 0 0| | 0 0 0| | | |1| | | |2| |
| | 2 3 1| | 0 0 0| | 0 0 0| | | |2| | | |1| |
| | 2 1 0| | 0 0 0| | 0 0 0| | | |3| | | |2| |
| | | | | |
| | 0 0 0| | 1 0 1| | 1 0 2| | | |1| | | |2| |
| | 0 0 0| | 4 1 0| | 0 1 1| | x | |3| | = | |3| |
| | 0 0 0| | 3 1 3| | 0 1 3| | | |2| | | |2| |
| | | | | |
| | 0 0 0| | 1 0 2| | 1 3 2| | | |1| | | |0| |
| | 0 0 0| | 0 1 4| | 2 1 3| | | |0| | | |2| |
| | 0 0 0| | 5 1 1| | 3 1 2| | | |2| | | |1| |
*/
const int N = 3;
const int bz = 3;
const int nonZeroes = 5;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
if (row.index() == 0) {
row.insert(row.index());
}
else if (row.index() == 1) {
row.insert(row.index());
row.insert(row.index() + 1);
}
else if (row.index() == 2) {
row.insert(row.index() - 1);
row.insert(row.index());
}
}
A[0][0][0][0]=3.0; A[1][1][0][0]=1.0; A[1][2][0][0]=1.0;
A[0][0][0][1]=1.0; A[1][1][0][1]=0.0; A[1][2][0][1]=0.0;
A[0][0][0][2]=2.0; A[1][1][0][2]=1.0; A[1][2][0][2]=2.0;
A[0][0][1][0]=2.0; A[1][1][1][0]=4.0; A[1][2][1][0]=0.0;
A[0][0][1][1]=3.0; A[1][1][1][1]=1.0; A[1][2][1][1]=1.0;
A[0][0][1][2]=1.0; A[1][1][1][2]=0.0; A[1][2][1][2]=1.0;
A[0][0][2][0]=2.0; A[1][1][2][0]=3.0; A[1][2][2][0]=0.0;
A[0][0][2][1]=1.0; A[1][1][2][1]=1.0; A[1][2][2][1]=1.0;
A[0][0][2][2]=0.0; A[1][1][2][2]=3.0; A[1][2][2][2]=3.0;
A[2][1][0][0]=1.0; A[2][2][0][0]=1.0;
A[2][1][0][1]=0.0; A[2][2][0][1]=3.0;
A[2][1][0][2]=2.0; A[2][2][0][2]=2.0;
A[2][1][1][0]=0.0; A[2][2][1][0]=2.0;
A[2][1][1][1]=1.0; A[2][2][1][1]=1.0;
A[2][1][1][2]=4.0; A[2][2][1][2]=3.0;
A[2][1][2][0]=5.0; A[2][2][2][0]=3.0;
A[2][1][2][1]=1.0; A[2][2][2][1]=1.0;
A[2][1][2][2]=1.0; A[2][2][2][2]=2.0;
Vector x(3);
x[0][0] = 1.0; x[1][0] = 1.0; x[2][0] = 1.0;
x[0][1] = 2.0; x[1][1] = 3.0; x[2][1] = 0.0;
x[0][2] = 3.0; x[1][2] = 2.0; x[2][2] = 2.0;
Vector b(3);
b[0][0] = 2.0; b[1][0] = 2.0; b[2][0] = 0.0;
b[0][1] = 1.0; b[1][1] = 3.0; b[2][1] = 2.0;
b[0][2] = 2.0; b[1][2] = 2.0; b[2][2] = 1.0;
// D_00 = A_00
auto D_00 = A[0][0];
auto D_00_inv = D_00;
D_00_inv.invert();
// D_11 = A_11 - L_10 D_00_inv U_01
// = A_11
auto D_11 = A[1][1];
auto D_11_inv = D_11;
D_11_inv.invert();
// D_22 = A_22 - A_20 D_00_inv A_02 - A_21 D_11_inv A_12
// = A_22 - A_21 D_11_inv A_12
auto D_22 = A[2][2] - A[2][1]*D_11_inv*A[1][2];
auto D_22_inv = D_22;
D_22_inv.invert();
// define: z = M^-1(b - Ax)
// where: M = (D + L_A) A^-1 (D + U_A)
// lower triangular solve: (E + L) y = b - Ax
Vector y(3);
// y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
// = D_00_inv*[b_0 - A_00*x_0]
auto rhs = b[0];
A[0][0].mmv(x[0], rhs);
D_00_inv.mv(rhs, y[0]);
// y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1 + A_12*x_2) - A_10*y_0)
// = D_11_inv*(b_1 - A_11*x_1)
rhs = b[1];
A[1][1].mmv(x[1], rhs);
A[1][2].mmv(x[2], rhs);
D_11_inv.mv(rhs, y[1]);
// y_2 = D_22_inv*(b_2 - (A_20*x_0 + A_21*x_1 + A_22*x_2) - (A_20*y_0 + A_21*y_1))
// = D_22_inv*(b_2 - (A_21*x_1 + A_22*x_2) - (A_21*y_1))
rhs = b[2];
A[2][1].mmv(x[1], rhs);
A[2][2].mmv(x[2], rhs);
A[2][1].mmv(y[1], rhs);
D_22_inv.mv(rhs, y[2]);
// upper triangular solve: (E + U) z = Ey
Vector z(3);
// z_2 = y_2
z[2] = y[2];
// z_1 = y_1 - D_11_inv*A_12*z_2
z[1] = y[1];
auto temp = D_11_inv*A[1][2];
temp.mmv(z[2], z[1]);
// z_0 = y_0 - D_00_inv(A_01*z_1 + A_02*z_2)
// z_0 = y_0
z[0] = y[0];
// x_k+1 = x_k + z
Vector new_x = x;
new_x += z;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
seqdilu.apply(x, b);
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
BOOST_CHECK_CLOSE(x[i][j], new_x[i][j], 1e-7);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(SeqDILUApplyIsEqualToDuneSeqILUApply, T, NumericTypes)
{
/*
Tests that applying the DILU preconditioner is equivalent to applying a ILU preconditioner
for a block diagonal matrix.
A x = b
| | 3 1| | 0 0| | | |1| | | |2| |
| | 2 1| | 0 0| | | |2| | | |1| |
| | x | | = | |
| | 0 0| |-1 0| | | |1| | | |3| |
| | 0 0| | 0 -1| | | |1| | | |4| |
*/
const int N = 2;
const int bz = 2;
const int nonZeroes = 2;
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
Matrix A(N, N, nonZeroes, Matrix::row_wise);
for (auto row = A.createbegin(); row != A.createend(); ++row) {
row.insert(row.index());
}
A[0][0][0][0]=3.0;
A[0][0][0][1]=1.0;
A[0][0][1][0]=2.0;
A[0][0][1][1]=1.0;
A[1][1][0][0]=-1.0;
A[1][1][1][1]=-1.0;
Dune::SeqDilu<Matrix, Vector, Vector> seqdilu(A);
Dune::SeqILU<Matrix, Vector, Vector> seqilu(A, 1.0);
Vector dilu_x(2);
dilu_x[0][0] = 1.0;
dilu_x[0][1] = 2.0;
dilu_x[1][0] = 1.0;
dilu_x[1][1] = 1.0;
Vector dilu_b(2);
dilu_b[0][0] = 2.0;
dilu_b[0][1] = 1.0;
dilu_b[1][0] = 3.0;
dilu_b[1][1] = 4.0;
Vector ilu_x(2);
ilu_x[0][0] = 1.0;
ilu_x[0][1] = 2.0;
ilu_x[1][0] = 1.0;
ilu_x[1][1] = 1.0;
Vector ilu_b(2);
ilu_b[0][0] = 2.0;
ilu_b[0][1] = 1.0;
ilu_b[1][0] = 3.0;
ilu_b[1][1] = 4.0;
seqdilu.apply(dilu_x, dilu_b);
seqilu.apply(ilu_x, ilu_b);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
BOOST_CHECK_CLOSE(dilu_x[i][j], ilu_x[i][j], 1e-7);
}
}
}