Add palace::GridFunction to unify mfem::ParGridFunction and mfem::ParComplexGridFunction

[sebastiangrimberg]

We have observed some issues with mfem::Memory aliasing when running on GPUs, which look like error messages of the form:

CUDA error: (cudaMemcpy(dst, src, bytes, cudaMemcpyHostToDevice)) failed with error:
 --> misaligned address
 ... in function: void* mfem::CuMemcpyHtoD(void*, const void*, size_t)
 ... in file: /data/home/hughcars/palace/build_g5_48xlarge/extern/mfem/general/cuda.cpp:116
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

One place this aliasing is used is mfem::ComplexGridFunction and mfem::ParComplexGridFunction. This PR replaces the aliases with two separate ParGridFunction objects for the real and imaginary parts, and resolves the errors in our testing. I'm not certain this is an error with MFEM's MemoryManager or just our use of it without adequate synchronization. But for now, keeping this objects stored separately seems to resolve the problems.

NOTE: This is on top of #193 and resolves the observed instances in testing for GPU support. #194 to follow this PR.

[hughcars]

MFEM_ASSERT(s.imag() == 0.0 || !HasImag(), "Cannot assign complex scalar to a non-complex grid function!");

Assigning an actually complex into a real grid function is almost certainly an error.

[sebastiangrimberg]

Added in 5f0543e

[hughcars]

The alternative here is something like

template <bool Complex = false>
class GridFunction
{
private:
  std::array<mfem::ParGridFunction, Complex ? 2 : 1> gfs;

public:
  GridFunction(mfem::ParFiniteElementSpace &fespace);
  GridFunction(FiniteElementSpace &fespace);

  // Get access to the real and imaginary grid function parts.
  const mfem::ParGridFunction &Real() const { return gfs[0]; }
  mfem::ParGridFunction &Real() { return gf[0]; }

  template <bool Enable = Complex, typename =std::enable_if_t<Complex>>
  const mfem::ParGridFunction &Imag() const
  {
    return gfs[1];
  }

  template <bool Enable = Complex, typename =std::enable_if_t<Complex>>
  mfem::ParGridFunction &Imag()
  {
    return gfs[1];
  }

  // Check if the grid function is suited for storing complex-valued fields.
  bool HasImag() const { return Complex; }

  // Get access to the underlying finite element space (match MFEM interface).
  mfem::FiniteElementSpace *FESpace() { return gfr.FESpace(); }
  const mfem::FiniteElementSpace *FESpace() const { return gfr.FESpace(); }
  mfem::ParFiniteElementSpace *ParFESpace() { return gfr.ParFESpace(); }
  const mfem::ParFiniteElementSpace *ParFESpace() const { return gfr.ParFESpace(); }

  // Set all entries equal to s.

  template <bool Enable = Complex, typename =std::enable_if_t<Enable>>
  GridFunction &operator=(std::complex<double> s);

  GridFunction &operator=(double s);

  // Scale all entries by s.
  GridFunction &operator*=(double s);

  // Transform for space update (for example on mesh change).
  void Update();

  // Get the associated MPI communicator.
  MPI_Comm GetComm() const { return ParFESpace()->GetComm(); }
};

where the enable_if bits remove the imaginary methods (and operators) in the case of non-complex. This would mirror the existing the split between the two grid function types, with all the corresponding dispatch.

I think the runtime flag here is a better approach though, as the template is only going to eliminate if branches at a very high level, so the template isn't really worth it imo.

[sebastiangrimberg]

I considered the template version but came to the same conclusion you have, that the runtime flag is just fine in this case from a performance and from an interface point of view.