Ifpack2: additive Schwarz initialize() should not require matrix values

[jhux2]

@trilinos/ifpack2 @vbrunini @csiefer2 @srajama1

Reported by Sierra/Aria developers:

"We found this a year or two ago [...] but I just noticed the disabled unit test reproducer we wrote at the time and it still fails so I figured I’d send it to you.

The root of the issue is that Ifpack2 additive schwarz with overlap appears to rely on the matrix values being populated before calling initialize(). We’ve worked around this by just calling initialize() right before our first solve, but it would be nice if it appropriately respected the convention that initialize() only requires the graph and compute() is where anything needing the values happens.

Unit test that reproduces the issue (must be run with > 1 mpi rank) is below.

TEST(Schwarz, DISABLED_OverlapBug)
{
  using Comm = Teuchos::MpiComm<int>;
  const auto comm = Teuchos::rcp(new Comm(MPI_COMM_WORLD));
  const int par_size = comm->getSize();
  const int par_rank = comm->getRank();
  const int num_local = 3;
  const int num_global = num_local * par_size;
 
  using GlobalOrdinal = long;
  using Matrix = Tpetra::CrsMatrix<double, int, GlobalOrdinal>;
  using Vector = Tpetra::Vector<double, int, GlobalOrdinal>;
  using Graph = Matrix::crs_graph_type;
  using Map = Graph::map_type;
  auto row_map = Teuchos::rcp(new Map(num_global, num_local, 0, comm));
  auto A = Teuchos::rcp(new Matrix(row_map, 3));
  auto x = Teuchos::rcp(new Vector(row_map));
  auto b = Teuchos::rcp(new Vector(row_map));
 
  b->sync_host();
  b->modify_host();
 
  // Initialize structure of A to tri-diagonal, all filled with 0's initially
  // Fill b with gid + 1
  for (int i = 0; i < num_local; ++i)
  {
    const GlobalOrdinal gid = num_local * par_rank + i;
    const double one = 0.0;
    const double zero = 0.0;
    A->insertGlobalValues(gid,
        Teuchos::ArrayView<const GlobalOrdinal>(&gid, 1),
        Teuchos::ArrayView<const double>(&one, 1));
    if (gid > 0)
    {
      const GlobalOrdinal gid_minus_one = gid - 1;
      A->insertGlobalValues(gid,
          Teuchos::ArrayView<const GlobalOrdinal>(&gid_minus_one, 1),
          Teuchos::ArrayView<const double>(&zero, 1));
    }
 
    if (gid < num_global - 1)
    {
      const GlobalOrdinal gid_plus_one = gid + 1;
      A->insertGlobalValues(gid,
          Teuchos::ArrayView<const GlobalOrdinal>(&gid_plus_one, 1),
          Teuchos::ArrayView<const double>(&zero, 1));
    }
 
    b->replaceLocalValue(i, gid + 1.);
  }
  A->fillComplete();
  b->sync_device();
 
  // Initialize the solver & preconditioner, this is supposed to be valid to do
  // once the Matrix structure is set, but before values are populated.
  auto paramlist = default_params();
  std::string mySolverInfoName = "solver";
  auto iterative_solver_parameters = paramlist.sublist(solver_sublist_name(), true);
  auto preconditioner_parameters = paramlist.sublist(preconditioner_sublist_name(), true);
  auto iterative_solver_type = paramlist.get<std::string>("solution method");
  std::string preconditioner_type = "SCHWARZ";
 
  // Don't hang on to the current matrix, preconditioner, or solver
  // longer than needed, to avoid storing two instances of each.
  Teuchos::RCP<const Matrix> matrix = A;
 
  using MultiVector = Tpetra::MultiVector<double, int, GlobalOrdinal>;
  using Operator = Tpetra::Operator<double, int, GlobalOrdinal>;
  using Belos_SolverFactory = Belos::SolverFactory<double, MultiVector, Operator>;
  Belos_SolverFactory solverFactory;
  auto iterative_solver = solverFactory.create(iterative_solver_type, Teuchos::null);
  using Teuchos::rcpFromRef;
  iterative_solver->setParameters(rcpFromRef(iterative_solver_parameters));
 
  Ifpack2::Factory factory;
  auto preconditioner = factory.create(preconditioner_type, matrix);
  preconditioner->setParameters(preconditioner_parameters);
  preconditioner->initialize();
 
  // Populate A with a basic diffusion matrix
  A->resumeFill();
  for (int i = 0; i < num_local; ++i)
  {
    const GlobalOrdinal gid = num_local * par_rank + i;
    const double mtwo = -2.0;
    const double one = 1.0;
    A->replaceGlobalValues(gid,
        Teuchos::ArrayView<const GlobalOrdinal>(&gid, 1),
        Teuchos::ArrayView<const double>(&mtwo, 1));
    if (gid > 0)
    {
      const GlobalOrdinal gid_minus_one = gid - 1;
      A->replaceGlobalValues(gid,
          Teuchos::ArrayView<const GlobalOrdinal>(&gid_minus_one, 1),
          Teuchos::ArrayView<const double>(&one, 1));
    }
 
    if (gid < num_global - 1)
    {
      const GlobalOrdinal gid_plus_one = gid + 1;
      A->replaceGlobalValues(gid,
          Teuchos::ArrayView<const GlobalOrdinal>(&gid_plus_one, 1),
          Teuchos::ArrayView<const double>(&one, 1));
    }
  }
  A->fillComplete();
 
  preconditioner->compute();
 
  using Belos_LinearProblem = Belos::LinearProblem<double, MultiVector, Operator>;
  auto linear_problem = Teuchos::rcp(new Belos_LinearProblem(matrix, x, b));
  linear_problem->setRightPrec(preconditioner);
  linear_problem->setProblem();
 
  iterative_solver->setProblem(linear_problem);
 
  // This throws when par_size > 1 iff the overlap level is > 0
  // Passes if the preconditioner initialization is done after the second call to A->fillComplte().
  iterative_solver->solve();
}

"

[github-actions]

This issue has had no activity for 365 days and is marked for closure. It will be closed after an additional 30 days of inactivity.
If you would like to keep this issue open please add a comment and/or remove the MARKED_FOR_CLOSURE label.
If this issue should be kept open even with no activity beyond the time limits you can add the label DO_NOT_AUTOCLOSE.
If it is ok for this issue to be closed, feel free to go ahead and close it. Please do not add any comments or change any labels or otherwise touch this issue unless your intention is to reset the inactivity counter for an additional year.

[vbrunini]

This is still a problem (and maybe related to #9606 ?) and should not get closed.

[github-actions]

This issue has had no activity for 365 days and is marked for closure. It will be closed after an additional 30 days of inactivity.
If you would like to keep this issue open please add a comment and/or remove the MARKED_FOR_CLOSURE label.
If this issue should be kept open even with no activity beyond the time limits you can add the label DO_NOT_AUTOCLOSE.
If it is ok for this issue to be closed, feel free to go ahead and close it. Please do not add any comments or change any labels or otherwise touch this issue unless your intention is to reset the inactivity counter for an additional year.

[vbrunini]

This got fixed by @brian-kelley in #10706