EigenproblemEquationSystem.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#ifdef MOOSE_MFEM_ENABLED
 
#include "EigenproblemEquationSystem.h"
#include "MFEMEigensolverBase.h"
#include "libmesh/int_range.h"
 
namespace Moose::MFEM
{
 
void
EigenproblemEquationSystem::ApplyEssentialBCs()
{
  _ess_tdof_lists.resize(1);
  _ess_markers.resize(1);
  mfem::ParGridFunction & trial_gf = *(_var_ess_constraints.at(0));
  trial_gf.Update();
  trial_gf = _gfuncs->GetRef(_trial_var_names.at(0));
  // A mesh without boundary elements, such as a periodic one, has no boundary attributes
  const auto & bdr_attributes = trial_gf.ParFESpace()->GetParMesh()->bdr_attributes;
  _ess_markers.at(0).SetSize(bdr_attributes.Size() ? bdr_attributes.Max() : 0, 0);
  // Set constrained DoF values on user-declared essential boundaries and collect their markers
  ApplyEssentialBC(_trial_var_names.at(0), trial_gf, _ess_markers.at(0));
  trial_gf.FESpace()->GetEssentialTrueDofs(_ess_markers.at(0), _ess_tdof_lists.at(0));
  CheckProblemIsHomogeneous();
}
 
void
EigenproblemEquationSystem::CheckProblemIsHomogeneous()
{
  // Reject nonzero Dirichlet BCs.
  mfem::ParGridFunction & trial_gf = *(_var_ess_constraints.at(0));
  mfem::Vector ess_values;
  trial_gf.GetTrueDofs(ess_values);
  ess_values.SetSubVectorComplement(_ess_tdof_lists.at(0), 0.0);
  mfem::real_t max_ess_value =
      mfem::GlobalLpNorm(mfem::infinity(), ess_values.Normlinf(), trial_gf.ParFESpace()->GetComm());
  // Roundoff guard. Zero coefficients project to exactly zero.
  if (max_ess_value > 10 * std::numeric_limits<mfem::real_t>::epsilon())
    mooseError("Essential boundary conditions on variable '",
               _trial_var_names.at(0),
               "' prescribe a nonzero value. "
               "An eigenproblem is homogeneous so only zero-valued essential boundary conditions "
               "are meaningful. Set the boundary coefficient to zero.");
}
 
void
EigenproblemEquationSystem::FormEigenproblemMatrix()
{
  auto & test_var_name = _test_var_names.at(0);
  auto blf = _blfs.Get(test_var_name);
 
  blf->EliminateEssentialBCDiag(_ess_markers.at(0), 1.0);
  blf->Finalize();
  _jacobian.Reset(blf->ParallelAssemble());
}
 
void
EigenproblemEquationSystem::FormMassMatrix(EigenRHSCoefficient rhs_coefficient)
{
  mfem::ParFiniteElementSpace * fespace = _test_pfespaces.at(0);
  std::unique_ptr<mfem::ParBilinearForm> m = std::make_unique<mfem::ParBilinearForm>(fespace);
 
  if (fespace->GetTypicalFE()->GetRangeType() == mfem::FiniteElement::SCALAR)
  {
    // VectorMassIntegrator reads vdim x vdim entries of the evaluated matrix without bounds checks.
    if (std::holds_alternative<mfem::MatrixCoefficient *>(rhs_coefficient))
    {
      const auto & mat_coef = *std::get<mfem::MatrixCoefficient *>(rhs_coefficient);
      if (mat_coef.GetHeight() != fespace->GetVDim() || mat_coef.GetWidth() != fespace->GetVDim())
        mooseError("A matrix rhs_coefficient on an H1 or L2 finite element space must be ",
                   fespace->GetVDim(),
                   "x",
                   fespace->GetVDim(),
                   " to match the number of components of the space, but it is ",
                   mat_coef.GetHeight(),
                   "x",
                   mat_coef.GetWidth(),
                   ".");
    }
    auto * const integ = std::visit(
        [](auto * coef) { return new mfem::VectorMassIntegrator(*coef); }, rhs_coefficient);
    // Otherwise, for scalar coefficients, VectorMassIntegrator defaults to the space dimension.
    integ->SetVDim(fespace->GetVDim());
    m->AddDomainIntegrator(integ);
  }
  else
    m->AddDomainIntegrator(std::visit(
        [](auto * coef) { return new mfem::VectorFEMassIntegrator(*coef); }, rhs_coefficient));
 
  m->Assemble();
  // Shift the eigenvalue corresponding to eliminated dofs to a large value. The BC DoFs on the
  // stiffness matrix are set to 1 and the mass matrix BC DoFs are set to a small value eps, such
  // that the eigenvalues associate with these DOFs are ~1/eps.
  m->EliminateEssentialBCDiag(_ess_markers.at(0), std::numeric_limits<mfem::real_t>::epsilon());
  m->Finalize();
  _mass_rhs.Reset(m->ParallelAssemble());
}
 
void
EigenproblemEquationSystem::BuildEigenproblemJacobian(mfem::BlockVector & trueX,
                                                      EigenRHSCoefficient rhs_coefficient)
{
  mooseAssert(_test_var_names.size() == 1 && (_test_var_names.size() == _trial_var_names.size()) &&
                  (_test_var_names.at(0) == _trial_var_names.at(0)),
              "Eigensolve is only supported for single-variable, square systems");
 
  height = trueX.Size();
  width = trueX.Size();
  ApplyEssentialBCs();
  FormEigenproblemMatrix();
  FormMassMatrix(rhs_coefficient);
}
 
void
EigenproblemEquationSystem::PrepareEigensolver(EigensolverBase & solver)
{
  solver.SetMassMatrix(*_mass_rhs);
  solver.SetOperator(*_jacobian);
}
 
} // namespace Moose::MFEM
 
#endif