mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-01-21 19:12:58 -06:00
Only add AMG preconditioners to factory if sensible.
Also add test using a new operator class, that would not compile without the change.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
6e87ec6266
commit
d910d42f6e
@ -132,6 +132,8 @@ list (APPEND TEST_DATA_FILES
|
||||
tests/wells_no_perforation.data
|
||||
tests/matr33.txt
|
||||
tests/rhs3.txt
|
||||
tests/matr33rep.txt
|
||||
tests/rhs3rep.txt
|
||||
tests/options_flexiblesolver.json
|
||||
tests/options_flexiblesolver_simple.json
|
||||
)
|
||||
|
||||
@ -288,17 +288,25 @@ private:
|
||||
const double w = prm.get<double>("relaxation", 1.0);
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
|
||||
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 = amgCriterion(prm);
|
||||
auto sargs = amgSmootherArgs<Smoother>(prm);
|
||||
return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
|
||||
}
|
||||
});
|
||||
|
||||
// Only add AMG preconditioners to the factory if the operator
|
||||
// is the overlapping schwarz operator. This could be extended
|
||||
// later, but at this point no other operators are compatible
|
||||
// with the AMG hierarchy construction.
|
||||
if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
|
||||
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 = amgCriterion(prm);
|
||||
auto sargs = amgSmootherArgs<Smoother>(prm);
|
||||
return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
|
||||
}
|
||||
});
|
||||
}
|
||||
|
||||
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
|
||||
assert(weightsCalculator);
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, comm);
|
||||
@ -355,74 +363,79 @@ private:
|
||||
const double w = prm.get<double>("relaxation", 1.0);
|
||||
return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
|
||||
if (smoother == "ILU0" || smoother == "ParOverILU0") {
|
||||
|
||||
// Only add AMG preconditioners to the factory if the operator
|
||||
// is an actual matrix operator.
|
||||
if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
|
||||
if (smoother == "ILU0" || smoother == "ParOverILU0") {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
#else
|
||||
using Smoother = SeqILU0<M, V, V>;
|
||||
using Smoother = SeqILU0<M, V, V>;
|
||||
#endif
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "Jac") {
|
||||
using Smoother = SeqJac<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SOR") {
|
||||
using Smoother = SeqSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SSOR") {
|
||||
using Smoother = SeqSSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "ILUn") {
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "Jac") {
|
||||
using Smoother = SeqJac<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SOR") {
|
||||
using Smoother = SeqSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SSOR") {
|
||||
using Smoother = SeqSSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "ILUn") {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
#else
|
||||
using Smoother = SeqILUn<M, V, V>;
|
||||
using Smoother = SeqILUn<M, V, V>;
|
||||
#endif
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
|
||||
}
|
||||
});
|
||||
doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
|
||||
if (smoother == "ILU0" || smoother == "ParOverILU0") {
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
|
||||
}
|
||||
});
|
||||
doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
|
||||
if (smoother == "ILU0" || smoother == "ParOverILU0") {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
#else
|
||||
using Smoother = SeqILU0<M, V, V>;
|
||||
using Smoother = SeqILU0<M, V, V>;
|
||||
#endif
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "Jac") {
|
||||
using Smoother = SeqJac<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "SOR") {
|
||||
using Smoother = SeqSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
// } else if (smoother == "GS") {
|
||||
// using Smoother = SeqGS<M, V, V>;
|
||||
// return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "SSOR") {
|
||||
using Smoother = SeqSSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "ILUn") {
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "Jac") {
|
||||
using Smoother = SeqJac<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "SOR") {
|
||||
using Smoother = SeqSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
// } else if (smoother == "GS") {
|
||||
// using Smoother = SeqGS<M, V, V>;
|
||||
// return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "SSOR") {
|
||||
using Smoother = SeqSSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else if (smoother == "ILUn") {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
using Smoother = SeqILU<M, V, V>;
|
||||
#else
|
||||
using Smoother = SeqILUn<M, V, V>;
|
||||
using Smoother = SeqILUn<M, V, V>;
|
||||
#endif
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
|
||||
}
|
||||
});
|
||||
doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
auto crit = amgCriterion(prm);
|
||||
Dune::Amg::Parameters parms;
|
||||
parms.setNoPreSmoothSteps(1);
|
||||
parms.setNoPostSmoothSteps(1);
|
||||
return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
|
||||
});
|
||||
return makeAmgPreconditioner<Smoother>(op, prm, true);
|
||||
} else {
|
||||
OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
|
||||
}
|
||||
});
|
||||
doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
auto crit = amgCriterion(prm);
|
||||
Dune::Amg::Parameters parms;
|
||||
parms.setNoPreSmoothSteps(1);
|
||||
parms.setNoPostSmoothSteps(1);
|
||||
return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
|
||||
});
|
||||
}
|
||||
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator);
|
||||
});
|
||||
|
||||
84
tests/matr33rep.txt
Normal file
84
tests/matr33rep.txt
Normal file
@ -0,0 +1,84 @@
|
||||
%%MatrixMarket matrix coordinate real general
|
||||
% ISTL_STRUCT blocked 3 3
|
||||
9 9 81
|
||||
1 1 1
|
||||
2 1 -1
|
||||
3 1 -1
|
||||
4 1 -1
|
||||
5 1 -1
|
||||
6 1 -1
|
||||
7 1 -1
|
||||
8 1 -1
|
||||
9 1 -1
|
||||
1 2 -1
|
||||
2 2 1
|
||||
3 2 -1
|
||||
4 2 -1
|
||||
5 2 -1
|
||||
6 2 -1
|
||||
7 2 -1
|
||||
8 2 -1
|
||||
9 2 -1
|
||||
1 3 -1
|
||||
2 3 -1
|
||||
3 3 1
|
||||
4 3 -1
|
||||
5 3 -1
|
||||
6 3 -1
|
||||
7 3 -1
|
||||
8 3 -1
|
||||
9 3 -1
|
||||
1 4 -1
|
||||
2 4 -1
|
||||
3 4 -1
|
||||
4 4 1
|
||||
5 4 -1
|
||||
6 4 -1
|
||||
7 4 -1
|
||||
8 4 -1
|
||||
9 4 -1
|
||||
1 5 -1
|
||||
2 5 -1
|
||||
3 5 -1
|
||||
4 5 -1
|
||||
5 5 1
|
||||
6 5 -1
|
||||
7 5 -1
|
||||
8 5 -1
|
||||
9 5 -1
|
||||
1 6 -1
|
||||
2 6 -1
|
||||
3 6 -1
|
||||
4 6 -1
|
||||
5 6 -1
|
||||
6 6 1
|
||||
7 6 -1
|
||||
8 6 -1
|
||||
9 6 -1
|
||||
1 7 -1
|
||||
2 7 -1
|
||||
3 7 -1
|
||||
4 7 -1
|
||||
5 7 -1
|
||||
6 7 -1
|
||||
7 7 1
|
||||
8 7 -1
|
||||
9 7 -1
|
||||
1 8 -1
|
||||
2 8 -1
|
||||
3 8 -1
|
||||
4 8 -1
|
||||
5 8 -1
|
||||
6 8 -1
|
||||
7 8 -1
|
||||
8 8 1
|
||||
9 8 -1
|
||||
1 9 -1
|
||||
2 9 -1
|
||||
3 9 -1
|
||||
4 9 -1
|
||||
5 9 -1
|
||||
6 9 -1
|
||||
7 9 -1
|
||||
8 9 -1
|
||||
9 9 1
|
||||
12
tests/rhs3rep.txt
Normal file
12
tests/rhs3rep.txt
Normal file
@ -0,0 +1,12 @@
|
||||
%%MatrixMarket matrix array real general
|
||||
% ISTL_STRUCT blocked 3 1
|
||||
9 1
|
||||
-1
|
||||
-1
|
||||
-1
|
||||
-3
|
||||
-3
|
||||
-3
|
||||
-3
|
||||
-3
|
||||
-3
|
||||
@ -232,6 +232,157 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
|
||||
test3(prm);
|
||||
}
|
||||
|
||||
|
||||
template<class Mat, class Vec>
|
||||
class RepeatingOperator : public Dune::AssembledLinearOperator<Mat, Vec, Vec>
|
||||
{
|
||||
public:
|
||||
using matrix_type = Mat;
|
||||
using domain_type = Vec;
|
||||
using range_type = Vec;
|
||||
using field_type = typename Vec::field_type;
|
||||
|
||||
Dune::SolverCategory::Category category() const override
|
||||
{
|
||||
return Dune::SolverCategory::sequential;
|
||||
}
|
||||
|
||||
RepeatingOperator(const Mat& matrix, const int repeats)
|
||||
: matrix_(matrix)
|
||||
, repeats_(repeats)
|
||||
{
|
||||
}
|
||||
|
||||
// y = A*x;
|
||||
virtual void apply(const Vec& x, Vec& y) const override
|
||||
{
|
||||
y = 0;
|
||||
applyscaleadd(1.0, x, y);
|
||||
}
|
||||
|
||||
// y += \alpha * A * x
|
||||
virtual void applyscaleadd(field_type alpha, const Vec& x, Vec& y) const override
|
||||
{
|
||||
Vec temp1 = x;
|
||||
Vec temp2 = x; // For size.
|
||||
temp2 = 0.0;
|
||||
for (int rr = 0; rr < repeats_; ++rr) {
|
||||
// mv below means: temp2 = matrix_ * temp1;
|
||||
matrix_.mv(temp1, temp2);
|
||||
temp1 = temp2;
|
||||
}
|
||||
temp2 *= alpha;
|
||||
y += temp2;
|
||||
}
|
||||
|
||||
virtual const matrix_type& getmat() const override
|
||||
{
|
||||
return matrix_;
|
||||
}
|
||||
|
||||
protected:
|
||||
const Mat& matrix_;
|
||||
const int repeats_;
|
||||
};
|
||||
|
||||
|
||||
template <int bz>
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>>
|
||||
testPrecRepeating(const boost::property_tree::ptree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
|
||||
{
|
||||
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
|
||||
using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
|
||||
Matrix matrix;
|
||||
{
|
||||
std::ifstream mfile(matrix_filename);
|
||||
if (!mfile) {
|
||||
throw std::runtime_error("Could not read matrix file");
|
||||
}
|
||||
readMatrixMarket(matrix, mfile);
|
||||
}
|
||||
Vector rhs;
|
||||
{
|
||||
std::ifstream rhsfile(rhs_filename);
|
||||
if (!rhsfile) {
|
||||
throw std::runtime_error("Could not read rhs file");
|
||||
}
|
||||
readMatrixMarket(rhs, rhsfile);
|
||||
}
|
||||
using Operator = RepeatingOperator<Matrix, Vector>;
|
||||
Operator op(matrix, 2);
|
||||
using PrecFactory = Opm::PreconditionerFactory<Operator>;
|
||||
|
||||
// Add no-oppreconditioner to factory for block size 1.
|
||||
PrecFactory::addCreator("nothing", [](const Operator&, const pt::ptree&, const std::function<Vector()>&) {
|
||||
return Dune::wrapPreconditioner<NothingPreconditioner<Vector>>();
|
||||
});
|
||||
|
||||
auto prec = PrecFactory::create(op, prm.get_child("preconditioner"));
|
||||
Dune::BiCGSTABSolver<Vector> solver(op, *prec, prm.get<double>("tol"), prm.get<int>("maxiter"), prm.get<int>("verbosity"));
|
||||
Vector x(rhs.size());
|
||||
Dune::InverseOperatorResult res;
|
||||
solver.apply(x, rhs, res);
|
||||
return x;
|
||||
}
|
||||
|
||||
void test1rep(const pt::ptree& prm)
|
||||
{
|
||||
const int bz = 1;
|
||||
auto sol = testPrecRepeating<bz>(prm, "matr33rep.txt", "rhs3rep.txt");
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>> expected {0.285714285714286,
|
||||
0.285714285714286,
|
||||
0.285714285714286,
|
||||
-0.214285714285714,
|
||||
-0.214285714285714,
|
||||
-0.214285714285714,
|
||||
-0.214285714285714,
|
||||
-0.214285714285714,
|
||||
-0.214285714285714};
|
||||
BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
|
||||
for (size_t i = 0; i < sol.size(); ++i) {
|
||||
for (int row = 0; row < bz; ++row) {
|
||||
BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void test3rep(const pt::ptree& prm)
|
||||
{
|
||||
const int bz = 3;
|
||||
auto sol = testPrecRepeating<bz>(prm, "matr33rep.txt", "rhs3rep.txt");
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>> expected {
|
||||
{0.285714285714286, 0.285714285714286, 0.285714285714286},
|
||||
{-0.214285714285714, -0.214285714285714, -0.214285714285714},
|
||||
{-0.214285714285714, -0.214285714285714, -0.214285714285714}
|
||||
};
|
||||
BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
|
||||
for (size_t i = 0; i < sol.size(); ++i) {
|
||||
for (int row = 0; row < bz; ++row) {
|
||||
BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE(TestWithRepeatingOperator)
|
||||
{
|
||||
pt::ptree prm;
|
||||
|
||||
// Read parameters.
|
||||
{
|
||||
std::ifstream file("options_flexiblesolver_simple.json");
|
||||
pt::read_json(file, prm);
|
||||
}
|
||||
|
||||
// Test with 1x1 block solvers.
|
||||
test1rep(prm);
|
||||
|
||||
// Test with 3x3 block solvers.
|
||||
test3rep(prm);
|
||||
}
|
||||
|
||||
|
||||
|
||||
#else
|
||||
|
||||
// Do nothing if we do not have at least Dune 2.6.
|
||||
|
||||
Loading…
Reference in New Issue
Block a user