Geometric multigrid with Galerkin triple product

[nilsfriess]

Hi, I'm currently trying to implement an algorithm that is in some sense equivalent to geometric multigrid in PETSc and now I'd like to use mfem to build the grid hierarchy and prolongation/restriction operators etc.

I followed ex26 to first implement a grid/operator hierarchy where the operators on each level are constructed by building new bilinear forms and rediscretising. That worked but now I'd like to instead construct the coarse operator by Galerkin projection, so A_coarse = P^T * A_fine * P, where P is the prolongation operator.
What would be the best approach to do this?

[tzanio]

This functionality is called RAP in MFEM (with $R = P^T$) and there are various functionalities for it depending on the type of the input objects.

For example, we have RAPOperator that takes as inputs general operators, or a version for dense matrices.

I suspect your inputs are HypreParMatrix objects, in which case you can use this.

[v-dobrev]

RAP is also implemented in the PETSc interface, see

mfem/linalg/petsc.hpp

Lines 576 to 580 in 9171062

	/// Returns the matrix P^t * A * P
	PetscParMatrix * RAP(PetscParMatrix *A, PetscParMatrix *P);
	/// Returns the matrix P^t * A * P
	PetscParMatrix * RAP(HypreParMatrix *A, PetscParMatrix *P);

Converting HypreParMatrix to PetscParMatrix can be done with:

mfem/linalg/petsc.hpp

Lines 377 to 385 in 9171062

	/** @brief Convert a HypreParMatrix @a ha to a PetscParMatrix in the given
	PETSc format @a tid. */
	/** The supported type ids are: Operator::PETSC_MATAIJ,
	Operator::PETSC_MATIS, Operator::PETSC_MATSHELL and
	Operator::PETSC_MATHYPRE
	@a ha can be destroyed unless tid == PETSC_MATSHELL or
	tid == PETSC_MATHYPRE */
	explicit PetscParMatrix(const HypreParMatrix *ha,
	Operator::Type tid = Operator::PETSC_MATAIJ);

[nilsfriess]

Thanks. I'm trying to use this as follows:

int main(int argc, char *argv[]) {
  PetscCall(PetscInitialize(&argc, &argv, nullptr, nullptr));

  if (argc < 2) {
    std::cout << "Usage: " << argv[0] << " meshfile\n";
    return 1;
  }

  std::string meshFile(argv[1]);
  mfem::Mesh mesh(meshFile, 1, 1);

  mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
  mesh.Clear();
  pmesh.UniformRefinement();
  pmesh.UniformRefinement();

  // Build FE space hierarchy
  mfem::H1_FECollection fec(1, 2);
  mfem::ParFiniteElementSpace coarseFespace(&pmesh, &fec);

  const int nLevels = 5;
  mfem::ParFiniteElementSpaceHierarchy fespaces(&pmesh, &coarseFespace, false, false);
  for (int l = 0; l < nLevels - 1; ++l)
    fespaces.AddUniformlyRefinedLevel();

  mfem::Array<int> essBdr(pmesh.bdr_attributes.Max());
  if (pmesh.bdr_attributes.Size())
    essBdr = 1;

  // Assemble hierarchy of operators
  std::vector<mfem::PetscParMatrix *> ops(nLevels, nullptr);

  mfem::ConstantCoefficient one(1.0);

  mfem::ParBilinearForm fineForm(&fespaces.GetFinestFESpace());
  fineForm.AddDomainIntegrator(new mfem::DiffusionIntegrator(one));
  fineForm.SetOperatorType(mfem::Operator::PETSC_MATAIJ);
  fineForm.Assemble();
  fineForm.Finalize();

  mfem::Array<int> essTdofs;
  fespaces.GetFinestFESpace().GetEssentialTrueDofs(essBdr, essTdofs);
  auto *a = new mfem::PetscParMatrix;
  fineForm.FormSystemMatrix(essTdofs, *a);
  ops[fespaces.GetFinestLevelIndex()] = a;

  for (int l = fespaces.GetFinestLevelIndex(); l > 0; --l) {
    mfem::PetscParMatrix pp(MPI_COMM_WORLD, fespaces.GetFESpaceAtLevel(l).GetProlongationMatrix(),
                            mfem::Operator::PETSC_MATAIJ);
    ops[l - 1] = mfem::RAP(ops[l], &pp);

    std::cout << "Size of P: " << pp.M() << "x" << pp.N() << "\n";
    std::cout << "Fine: " << ops[l]->M() << ", coarse: " << ops[l - 1]->M() << "\n";
  }

  PetscCall(PetscFinalize());
}

But after the first loop I get an error in RAP saying that the size of the prolongation matrix and the size of A are not compatible. Then I added the debug prints, and it seems the prolongation matrix is square; I had assumed it to be rectangular. Is this the wrong approach?

[nilsfriess]

@tzanio @v-dobrev Do you have some suggestions for my problem?
Is fespaces.GetFESpaceAtLevel(l).GetProlongationMatrix() not the correct way to get the prolongation matrix at level l?

[v-dobrev]

ParFiniteElementSpace::GetProlongationMatrix() is the wrong prolongation matrix. The prolongation you want is FiniteElementSpaceHierarchy::GetProlongationAtLevel(int level). However, this is a matrix-free version and it is not easy to convert to a PetscParMatrix. To get the prolongation matrix as a HypreParMatrix, you can use the method ParFiniteElementSpace::GetTrueTransferOperator where you need to set the type of the second argument to Operator::Hypre_ParCSR, e.g. something like this:

   for (int l = fespaces.GetFinestLevelIndex(); l > 0; --l)
   {
      OperatorHandle P(Operator::Hypre_ParCSR);
      fespaces.GetFESpaceAtLevel(l).GetTrueTransferOperator(fespaces.GetFESpaceAtLevel(l-1), P);
      HypreParMatrix *hP = P.As<HypreParMatrix>();
      // convert hP to PetscParMatrix ...
      // ...
   }

Let us know if you run into other issues with this approach.

[nilsfriess]

That worked, thank you!