amgcl: Stokes problem gives NaN by AMG but GMRES works fine

Sorry about the confusion. The system is a “bad” system where the boundary condtions are not propertely setup. The V5 works fine now after fixing the boundary condition. I will post my latest benchmark soon.

I see. Well, I am out of ideas then. The approach that works (SchurPC with dummy preconditioner for the pressure part) is not strictly multigrid, so I guess the increase in the number of iterations is not unexpected. You can try to increase the number of iterations for the pressure solver (since the unpreconditioned solver is not a very effective one). It somewhat helps with the matrix you provided:

./schur_pressure_correction \
    -A problems/issue-144/condensed/A.mm \
    -f problems/issue-144/condensed/b.mm \
    -m '>15552' -p \
    solver.type=fgmres \
    precond.usolver.solver.maxiter=1 \
    precond.psolver.precond.class=dummy \
    precond.usolver.precond.{class=relaxation,type=ilu0} \
    precond.psolver.solver.maxiter=5
Schur complement (two-stage preconditioner)
  unknowns: 19392(3840)
  nonzeros: 1102992

Iterations: 39
Error:      7.89528e-09

[Profile:                2.831 s] (100.00%)
[  reading:              2.174 s] ( 76.78%)
[  schur_complement:     0.657 s] ( 23.21%)
[    setup:              0.057 s] (  2.03%)
[    solve:              0.599 s] ( 21.15%)

Another thing to try is to disable the use of approximate schur complement:

./schur_pressure_correction \
    -A problems/issue-144/condensed/A.mm \
    -f problems/issue-144/condensed/b.mm \
    -m '>15552' -p \
    solver.type=fgmres \
    precond.usolver.solver.maxiter=1 \
    precond.psolver.precond.class=dummy \
    precond.usolver.precond.{class=relaxation,type=ilu0} \
    precond.psolver.solver.maxiter=5 \
    precond.approx_schur=0
Schur complement (two-stage preconditioner)
  unknowns: 19392(3840)
  nonzeros: 1102992

Iterations: 20
Error:      9.33415e-09

[Profile:                3.460 s] (100.00%)
[  reading:              2.162 s] ( 62.48%)
[  schur_complement:     1.298 s] ( 37.52%)
[    setup:              0.057 s] (  1.64%)
[    solve:              1.240 s] ( 35.85%)

Finally, since the pressure part in the condensed matrix does have a nonzero diagonal, you can replace dummy preconditioner with some simple single level one like spai0 or damped_jacobi:

./schur_pressure_correction \
    -A problems/issue-144/condensed/A.mm \
    -f problems/issue-144/condensed/b.mm \
    -m '>15552' -p \
    solver.type=fgmres \
    precond.usolver.solver.maxiter=1 \
    precond.psolver.precond.{class=relaxation,type=damped_jacobi} \
    precond.usolver.precond.{class=relaxation,type=ilu0} \
    precond.psolver.solver.maxiter=5 \
    precond.approx_schur=0
Schur complement (two-stage preconditioner)
  unknowns: 19392(3840)
  nonzeros: 1102992

Iterations: 24
Error:      2.14533e-09

[Profile:                3.910 s] (100.00%)
[  reading:              2.377 s] ( 60.79%)
[  schur_complement:     1.533 s] ( 39.21%)
[    setup:              0.055 s] (  1.41%)
[    solve:              1.477 s] ( 37.78%)

Not sure how any of these will scale though.

Also, the USolver and PSolver template parameters for the schur_pressure_correction are in the wrong order. The following works for me (I’ve changed the inner iterative solvers to bicgstab as that is what was used in my example above):

#include <iostream>
#include <string>

#ifdef _OPENMP
#include <omp.h>
#endif

#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/preconditioner/dummy.hpp>
#include <amgcl/preconditioner/schur_pressure_correction.hpp>
#include <amgcl/solver/fgmres.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/coarsening/aggregation.hpp>
#include <amgcl/coarsening/ruge_stuben.hpp>
#include <amgcl/relaxation/as_preconditioner.hpp>
#include <amgcl/relaxation/damped_jacobi.hpp>
#include <amgcl/relaxation/ilu0.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/io/binary.hpp>
#include <amgcl/io/mm.hpp>
#include <amgcl/profiler.hpp>

typedef amgcl::backend::builtin<double> Backend;


int main(int argc, char* argv[]) {
    using std::string;
    using std::vector;

    using namespace amgcl;
    namespace io = amgcl::io;

    size_t rows;
    vector<ptrdiff_t> ptr, col;
    vector<double> val, rhs;
    std::vector<char> pm;


    size_t cols;
    std::tie(rows, cols) = io::mm_reader("A.mm")(ptr, col, val);

    size_t n, m;
    std::tie(n, m) = io::mm_reader("b.mm")(rhs);

    typedef
        make_solver<
            preconditioner::schur_pressure_correction<
                make_solver<
                    amg<
                        Backend,
                        coarsening::smoothed_aggregation,
                        relaxation::ilu0
                        >,
                    solver::bicgstab<Backend>
                    >,
                make_solver<
                    preconditioner::dummy<Backend>,
                    solver::bicgstab<Backend>
                    >
                >,
            solver::fgmres<Backend>
        > Solver;

    //Parameter setup
    Solver::params prm;
    double tol = 1.0e-8;
    prm.solver.tol = tol;

    prm.precond.usolver.solver.tol = tol * 10;
    prm.precond.psolver.solver.tol = tol * 10;
    prm.precond.psolver.solver.maxiter = 1;
    prm.precond.usolver.solver.maxiter = 1;

    prm.precond.pmask.resize(n, 0);
    for (size_t i = 0; i < rows; ++i)
        if (i >= 3600) prm.precond.pmask[i] = 1;


    //Solve the system
    Solver solve(std::tie(rows, ptr, col, val), prm);
    std::cout << solve << std::endl;
    std::vector<double> x(n, 0.0);

    int iters;
    double error;
    std::tie(iters, error) = solve(rhs, x);
    printf("#Iters=%d Error=%lf\n", iters, error);
}

The result is:

Solver
======
Type:             FGMRES(30)
Unknowns:         4080
Memory footprint: 1.94 M

Preconditioner
==============
Schur complement (two-stage preconditioner)
  unknowns: 4080(480)
  nonzeros: 201072

#Iters=41 Error=0.000000