EquationSystem.h

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

#pragma once

#include "libmesh/ignore_warnings.h"
#include "mfem/miniapps/common/pfem_extras.hpp"
#include "libmesh/restore_warnings.h"
#include "MFEMIntegratedBC.h"
#include "MFEMEssentialBC.h"
#include "MFEMContainers.h"
#include "MFEMKernel.h"
#include "MFEMMixedBilinearFormKernel.h"
#include "ScaleIntegrator.h"
#include "NLScaleIntegrator.h"

namespace Moose::MFEM
{
class LinearSolverBase;

class EquationSystem : public mfem::Operator
{

public:
  EquationSystem() = default;
  ~EquationSystem() override;

  virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
  virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
  virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);

  virtual void Init(GridFunctions & gridfunctions,
                    ComplexGridFunctions & cmplx_gridfunctions,
                    mfem::AssemblyLevel assembly_level);
  void FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
  void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
  virtual void ComputeNonlinearResidual(const mfem::Vector & u, mfem::Vector & residual) const;
  mfem::Operator & GetGradient(const mfem::Vector & u) const override;

  virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;

  void SetSolverRequiresGradient(bool requires_gradient)
  {
    _solver_requires_gradient = requires_gradient;
  }

  // Test variables are associated with linear forms,
  // whereas trial variables are associated with gridfunctions.
  const std::vector<std::string> & GetTrialVarNames() const { return _trial_var_names; }
  const std::vector<std::string> & GetTestVarNames() const { return _test_var_names; }

  bool Nonlinear() const { return _non_linear; }

  void PrepareLinearSolver(LinearSolverBase & solver);

protected:
  virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
  virtual void AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name);
  virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
  virtual void SetTrialVariableNames();

  void DeleteHBlocks();

  void DeleteJacobianBlocks();

  bool VectorContainsName(const std::vector<std::string> & the_vector,
                          const std::string & name) const;

  virtual void ApplyEssentialBC(const std::string & var_name,
                                mfem::ParGridFunction & trial_gf,
                                mfem::Array<int> & global_ess_markers);
  virtual void ApplyEssentialBCs();
  virtual void EliminateCoupledVariables();
  virtual void BuildLinearForms();
  virtual void BuildNonlinearForms();
  virtual void BuildBilinearForms();
  virtual void BuildMixedBilinearForms();
  virtual void BuildEquationSystem();

  virtual void FormLinearSystem(mfem::OperatorHandle & op,
                                mfem::BlockVector & trueX,
                                mfem::BlockVector & trueRHS);
  using mfem::Operator::FormSystemOperator;
  virtual void FormSystemOperator(mfem::OperatorHandle & op,
                                  mfem::BlockVector & trueX,
                                  mfem::BlockVector & trueRHS);
  virtual void FormSystemMatrix(mfem::OperatorHandle & op,
                                mfem::BlockVector & trueX,
                                mfem::BlockVector & trueRHS);
  void FormJacobianMatrix(const mfem::Vector & u);

  template <class FormType>
  void ApplyDomainBLFIntegrators(
      const std::string & trial_var_name,
      const std::string & test_var_name,
      std::shared_ptr<FormType> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
      std::optional<mfem::real_t> scale_factor = std::nullopt);

  void ApplyDomainLFIntegrators(
      const std::string & test_var_name,
      std::shared_ptr<mfem::ParLinearForm> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);

  void ApplyDomainNLFIntegrators(
      const std::string & test_var_name,
      std::shared_ptr<mfem::ParNonlinearForm> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
      std::optional<mfem::real_t> scale_factor = std::nullopt);

  template <class FormType>
  void ApplyBoundaryBLFIntegrators(
      const std::string & trial_var_name,
      const std::string & test_var_name,
      std::shared_ptr<FormType> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
          integrated_bc_map,
      std::optional<mfem::real_t> scale_factor = std::nullopt);

  void ApplyBoundaryLFIntegrators(
      const std::string & test_var_name,
      std::shared_ptr<mfem::ParLinearForm> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
          integrated_bc_map);

  void ApplyBoundaryNLFIntegrators(
      const std::string & test_var_name,
      std::shared_ptr<mfem::ParNonlinearForm> form,
      NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
          integrated_bc_map,
      std::optional<mfem::real_t> scale_factor = std::nullopt);

  virtual bool Complex() const { return false; }

  std::vector<std::string> _coupled_var_names;
  std::vector<std::string> _trial_var_names;
  std::vector<std::string> _eliminated_var_names;
  Moose::MFEM::GridFunctions _eliminated_variables;
  std::vector<std::string> _test_var_names;
  std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
  std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;

  // Components of weak form, named according to test variable
  NamedFieldsMap<mfem::ParBilinearForm> _blfs;
  NamedFieldsMap<mfem::ParLinearForm> _lfs;
  NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
  NamedFieldsMap<NamedFieldsMap<mfem::ParMixedBilinearForm>> _mblfs; // named according to trial var

  std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
  std::vector<mfem::Array<int>> _ess_tdof_lists;

  mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks, _jacobian_blocks;
  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> _kernels_map;
  NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> _integrated_bc_map;
  NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;

  // Operator handle for the jacobian
  mutable mfem::OperatorHandle _jacobian;
  // Operator handle for the linear components of the system operator
  mutable mfem::OperatorHandle _linear_operator;
  mfem::AssemblyLevel _assembly_level;

  // Pointer to GridFunctions to enable updates during nonlinear iterations
  Moose::MFEM::GridFunctions * _gfuncs;
  // Array storing block offsets of solution and residual vector
  mfem::Array<int> _block_true_offsets;
  // Boolean indicating if EquationSystem contains nonlinear integrators
  bool _non_linear = false;
  // Whether a nonlinear solver exists and whether it requires Jacobian/gradient information.
  bool _solver_requires_gradient = false;

private:
  friend class EquationSystemProblemOperator;
  friend class ::MFEMProblemSolve;
  using mfem::Operator::RecoverFEMSolution;
};

template <class FormType>
void
EquationSystem::ApplyDomainBLFIntegrators(
    const std::string & trial_var_name,
    const std::string & test_var_name,
    std::shared_ptr<FormType> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
    std::optional<mfem::real_t> scale_factor)
{
  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
  {
    auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
    for (auto & kernel : kernels)
    {
      mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();

      if (integ)
      {
        if (scale_factor.has_value())
          integ = new ScaleIntegrator(integ, scale_factor.value(), true);
        kernel->isSubdomainRestricted()
            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
            : form->AddDomainIntegrator(std::move(integ));
      }
    }
  }
}

template <class FormType>
void
EquationSystem::ApplyBoundaryBLFIntegrators(
    const std::string & trial_var_name,
    const std::string & test_var_name,
    std::shared_ptr<FormType> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
        integrated_bc_map,
    std::optional<mfem::real_t> scale_factor)
{
  if (integrated_bc_map.Has(test_var_name) &&
      integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
  {
    auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
    for (auto & bc : bcs)
    {
      mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();

      if (integ)
      {
        if (scale_factor.has_value())
          integ = new ScaleIntegrator(integ, scale_factor.value(), true);
        bc->isBoundaryRestricted()
            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
            : form->AddBoundaryIntegrator(std::move(integ));
      }
    }
  }
}

} // namespace Moose::MFEM

#endif