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 CoefficientManager;
 
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;
  mfem::OperatorHandle & GetLinearOperator() const { return _linear_operator; };
 
  virtual void SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const;
 
  void SetGradientRequired(bool requires_gradient) { _gradient_required = requires_gradient; }
 
  void SetCoefficientManager(CoefficientManager & coefficients)
  {
    _coefficient_manager = &coefficients;
  }
 
  // 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; }
 
  const mfem::Array<int> & GetBlockOffsets() const { return _block_true_offsets; }
 
  mfem::ParBilinearForm & GetBilinearForm(const std::string & test_var_name)
  {
    return _blfs.GetRef(test_var_name);
  }
 
  mfem::ParGridFunction & GetGridFunction(const std::string & trial_var_name)
  {
    return _gfuncs->GetRef(trial_var_name);
  }
 
  bool Nonlinear() const { return _non_linear; }
 
  const mfem::Vector & GetLinearizationPoint() const;
 
  std::shared_ptr<mfem::ParBilinearForm>
  BuildBilinearFormForFESpace(const std::string & var_name,
                              mfem::ParFiniteElementSpace & fespace,
                              mfem::AssemblyLevel assembly_level);
 
  std::shared_ptr<mfem::ParNonlinearForm>
  BuildNonlinearFormForFESpace(const std::string & var_name,
                               mfem::ParFiniteElementSpace & fespace,
                               mfem::AssemblyLevel assembly_level);
 
  mfem::Array<int> & GetEssentialBoundaryMarkers(const std::string & var_name);
 
  virtual bool IsComplex() const { return false; }
  virtual bool IsEigen() const { return false; }
  virtual bool IsTimeDependent() const { return false; }
  bool IsMultivariate() const { return _test_var_names.size() > 1; }
  bool IsNonlinear() const { return _non_linear; }
 
  virtual void BuildEquationSystem();
 
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 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);
 
  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;
  std::vector<mfem::Array<int>> _ess_markers;
 
  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;
  // Non-owning pointer to the vector passed into the most recent GetGradient() call.
  mutable const mfem::Vector * _linearization_point = nullptr;
  // Boolean indicating if EquationSystem contains nonlinear integrators
  bool _non_linear = false;
  // Whether an external object (e.g. solver) requires Jacobian/gradient information.
  bool _gradient_required = false;
  // Coefficient manager notified when trial variables are updated, so that stored projections
  // of solution-dependent coefficients are invalidated.
  CoefficientManager * _coefficient_manager = nullptr;
 
private:
  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->isDGKernel() ? form->AddInteriorFaceIntegrator(std::move(integ))
        : 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->isDGBC() ? bc->isBoundaryRestricted()
                           ? form->AddBdrFaceIntegrator(std::move(integ), bc->getBoundaryMarkers())
                           : form->AddBdrFaceIntegrator(std::move(integ))
        : bc->isBoundaryRestricted()
            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
            : form->AddBoundaryIntegrator(std::move(integ));
      }
    }
  }
}
 
} // namespace Moose::MFEM
 
#endif