Hi everyone,

I am using MatrixFree framework with geometric multigrid in dealii. I am 
trying to implement a matrix-based coarse grid solver that uses PETSc AMG 
as a preconditioner. However, I have encountered a problem: under certain 
conditions (which I describe below), 
the MatrixFreeTools::compute_matrix() function fails to compute the coarse 
system matrix. I reproduced the same error using a modified version of the 
tutorial code of step-75. The error is as follows:

```
[3]PETSC ERROR: --------------------- Error Message 
--------------------------------------------------------------
[3]PETSC ERROR: Argument out of range
[3]PETSC ERROR: Inserting a new nonzero at global row/column (30, 35) into 
matrix
[3]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[3]PETSC ERROR: Petsc Release Version 3.21.4, Jul 30, 2024 
[3]PETSC ERROR: step-75 on a arch-linux-c-opt named paragon by dg615 Fri 
Nov 28 22:32:54 2025
[3]PETSC ERROR: Configure options --with-cc=mpicc --with-cxx=mpicxx 
--with-fc=mpif90 --with-debugging=0 COPTFLAGS="-O3 -march=native 
-mtune=native" CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3 
-march=native -mtune=native" --with-blaslapack-dir=/opt/intel/mkl 
--download-superlu_dist --download-superlu --download-metis 
--download-parmetis --download-mumps --download-scalapack --download-hypre 
--download-hdf5 --download-fftw --download-suitesparse --download-cmake
[3]PETSC ERROR: #1 MatSetValues_MPIAIJ() at 
/home/dg615/proj/petsc/petsc-3.21.4/src/mat/impls/aij/mpi/mpiaij.c:592
[3]PETSC ERROR: #2 MatAssemblyEnd_MPIAIJ() at 
/home/dg615/proj/petsc/petsc-3.21.4/src/mat/impls/aij/mpi/mpiaij.c:777
[3]PETSC ERROR: #3 MatAssemblyEnd() at 
/home/dg615/proj/petsc/petsc-3.21.4/src/mat/interface/matrix.c:5820
```

The block of code that produces the error is the following:

```
  template <int dim, typename number>
  const PETScWrappers::MPI::SparseMatrix &
  LaplaceOperator<dim, number>::get_system_matrix() const
  {
    if (system_matrix.m() == 0 && system_matrix.n() == 0)
      {
        const auto &dof_handler = this->matrix_free.get_dof_handler();

        // MODIFIED starts
        // 
=======================================================================
        const IndexSet &owned_dofs = dof_handler.locally_owned_dofs();
        const IndexSet &relevant_dofs = DoFTools::
extract_locally_relevant_dofs(dof_handler);
        const IndexSet &used_dofs = owned_dofs; // using relevant_dofs here 
gives an error
        DynamicSparsityPattern dsp(used_dofs);

        DoFTools::make_sparsity_pattern(dof_handler, dsp, this->constraints)
;

        SparsityPattern sp;
        sp.copy_from(dsp);
        system_matrix.reinit(
          used_dofs,
          used_dofs,
          sp,
          dof_handler.get_triangulation().get_mpi_communicator()
        );
        // 
=======================================================================
        // MODIFIED ends

        // The error is thrown from the function below
        MatrixFreeTools::compute_matrix(
          matrix_free,
          constraints,
          system_matrix,
          &LaplaceOperator::do_cell_integral_local,
          this);
      }

    return this->system_matrix;
  }
```
I've also attached the modified step-75.cc file. You can find the modified 
area of the code by searching MODIFIED keyword. I also changed the code to 
run with PETSc (originally Trilinos; NOTE: you might have to modify the 
CMakeLists.txt file not to require Trilinos if you don't have Trilinos 
installed (I don't) and are pasting the step-75.cc file into the tutorial 
folder).

Error condition: in my code the coarse grid is a cuboid with some number of 
global refinements, the error occurs when the number of those refinements 
is > 2. In step-75 modified code the error occurs when the variable 
`minLevel` is set to anything above 0.

I have tried using locally_relevant_dofs instead of locally_owned_dofs but 
in that case i get a different error:

```
An error occurred in line <203> of file 
</home/dg615/proj/petsc/dealii-9.7.0/source/lac/petsc_parallel_sparse_matrix.cc>
 
in function
    void dealii::PETScWrappers::MPI::SparseMatrix::do_reinit(MPI_Comm, 
const dealii::IndexSet&, const dealii::IndexSet&, const 
SparsityPatternType&) [with SparsityPatternType = dealii::SparsityPattern; 
MPI_Comm = ompi_communicator_t*]
The violated condition was: 
    local_rows.is_contiguous() && local_columns.is_contiguous()
Additional information: 
    PETSc only supports contiguous row/column ranges
```

Can someone please let me know how I can fix this issue? Is it internal to 
dealii/PETSc or is it something I've missed? I know I could assemble the 
coarse matrices using the standard matrix-based assembly logic, but I was 
wondering if it could be done entirely in the MF infrastructure.

Thank you.

Best wishes,
Davit Gyulamiryan

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to [email protected].
To view this discussion visit 
https://groups.google.com/d/msgid/dealii/ceed97dd-b786-4874-9756-2a39297f6a2fn%40googlegroups.com.

Reply via email to