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"

namespace Moose::MFEM
{

class EquationSystem : public mfem::Operator
{

public:
  friend class EquationSystemProblemOperator;
  friend class TimeDomainEquationSystemProblemOperator;

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

  virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
  virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);

  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(Moose::MFEM::GridFunctions & gridfunctions, mfem::AssemblyLevel assembly_level);

  virtual void BuildLinearForms();

  virtual void ApplyEssentialBC(const std::string & test_var_name,
                                mfem::ParGridFunction & trial_gf,
                                mfem::Array<int> & global_ess_markers);
  virtual void ApplyEssentialBCs();

  virtual void EliminateCoupledVariables();

  virtual void BuildBilinearForms();
  virtual void BuildMixedBilinearForms();
  virtual void BuildEquationSystem();

  void AssembleJacobian(
      Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> & jac_blfs,
      Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>> &
          jac_mblfs,
      Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> & rhs_lfs,
      std::vector<mfem::Array<int>> & ess_tdof_lists,
      std::vector<std::unique_ptr<mfem::ParGridFunction>> & var_ess_constraints,
      mfem::OperatorHandle & op,
      mfem::BlockVector & trueX,
      mfem::BlockVector & trueRHS);

  virtual void FormLinearSystem(mfem::OperatorHandle & op,
                                mfem::BlockVector & trueX,
                                mfem::BlockVector & trueRHS);

  virtual void
  FormSystem(mfem::OperatorHandle & op, mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
  virtual void FormLegacySystem(mfem::OperatorHandle & op,
                                mfem::BlockVector & trueX,
                                mfem::BlockVector & trueRHS);

  virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);

  void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;

  mfem::Operator & GetGradient(const mfem::Vector & u) const override;

  virtual void RecoverFEMSolution(mfem::BlockVector & trueX,
                                  Moose::MFEM::GridFunctions & gridfunctions);

  std::vector<mfem::Array<int>> _ess_tdof_lists;

  const std::vector<std::string> & TrialVarNames() const { return _trial_var_names; }
  const std::vector<std::string> & TestVarNames() const { return _test_var_names; }

private:
  using mfem::Operator::RecoverFEMSolution;

  virtual void SetTrialVariableNames();

protected:
  void DeleteAllBlocks();

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

  // Test variables are associated with LinearForms,
  // whereas trial variables are associated with gridfunctions.

  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
  Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
  Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _lfs;
  Moose::MFEM::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
  Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
      _mblfs; // named according to trial variable

  template <class FormType>
  void ApplyDomainBLFIntegrators(
      const std::string & trial_var_name,
      const std::string & test_var_name,
      std::shared_ptr<FormType> form,
      Moose::MFEM::NamedFieldsMap<
          Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);

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

  template <class FormType>
  void ApplyBoundaryBLFIntegrators(
      const std::string & trial_var_name,
      const std::string & test_var_name,
      std::shared_ptr<FormType> form,
      Moose::MFEM::NamedFieldsMap<
          Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
          integrated_bc_map);

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

  std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;

  mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks;

  Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
      _kernels_map;
  Moose::MFEM::NamedFieldsMap<
      Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>
      _integrated_bc_map;
  Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;

  mutable mfem::OperatorHandle _jacobian;

  mfem::AssemblyLevel _assembly_level;
};

template <class FormType>
void
EquationSystem::ApplyDomainBLFIntegrators(
    const std::string & trial_var_name,
    const std::string & test_var_name,
    std::shared_ptr<FormType> form,
    Moose::MFEM::NamedFieldsMap<
        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
{
  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 != nullptr)
      {
        kernel->isSubdomainRestricted()
            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
            : form->AddDomainIntegrator(std::move(integ));
      }
    }
  }
}

inline void
EquationSystem::ApplyDomainLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParLinearForm> form,
    Moose::MFEM::NamedFieldsMap<
        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
{
  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
  {
    auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
    for (auto & kernel : kernels)
    {
      mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
      if (integ != nullptr)
      {
        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,
    Moose::MFEM::NamedFieldsMap<
        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
        integrated_bc_map)
{
  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 != nullptr)
      {
        bc->isBoundaryRestricted()
            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
            : form->AddBoundaryIntegrator(std::move(integ));
      }
    }
  }
}

inline void
EquationSystem::ApplyBoundaryLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParLinearForm> form,
    Moose::MFEM::NamedFieldsMap<
        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
        integrated_bc_map)
{
  if (integrated_bc_map.Has(test_var_name) &&
      integrated_bc_map.Get(test_var_name)->Has(test_var_name))
  {
    auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
    for (auto & bc : bcs)
    {
      mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
      if (integ != nullptr)
      {
        bc->isBoundaryRestricted()
            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
            : form->AddBoundaryIntegrator(std::move(integ));
      }
    }
  }
}

class TimeDependentEquationSystem : public EquationSystem
{
public:
  TimeDependentEquationSystem(const Moose::MFEM::TimeDerivativeMap & time_derivative_map);

  virtual void Init(Moose::MFEM::GridFunctions & gridfunctions,
                    mfem::AssemblyLevel assembly_level) override;

  virtual void SetTimeStep(mfem::real_t dt);
  virtual void UpdateEquationSystem();

  virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel) override;
  virtual void BuildBilinearForms() override;
  virtual void BuildMixedBilinearForms() override;
  virtual void ApplyEssentialBCs() override;
  virtual void EliminateCoupledVariables() override;
  virtual void FormLegacySystem(mfem::OperatorHandle & op,
                                mfem::BlockVector & truedXdt,
                                mfem::BlockVector & trueRHS) override;
  virtual void FormSystem(mfem::OperatorHandle & op,
                          mfem::BlockVector & truedXdt,
                          mfem::BlockVector & trueRHS) override;

  template <class FormType>
  void ScaleAndAddBLFIntegrators(std::shared_ptr<FormType> source_blf,
                                 std::shared_ptr<FormType> target_blf,
                                 mfem::real_t scale_factor)
  {
    // Add and scale all domain integrators existing on source_blf to target_blf
    auto domain_integrators = source_blf->GetDBFI();
    auto domain_markers = source_blf->GetDBFI_Marker();
    for (int i = 0; i < domain_integrators->Size(); ++i)
      target_blf->AddDomainIntegrator(
          new ScaleIntegrator(*domain_integrators[i], scale_factor, false), *(*domain_markers[i]));
    // Add and scale all boundary integrators existing on source_blf to target_blf
    auto boundary_integrators = source_blf->GetBBFI();
    auto boundary_markers = source_blf->GetBBFI_Marker();
    for (int i = 0; i < boundary_integrators->Size(); ++i)
      target_blf->AddBoundaryIntegrator(
          new ScaleIntegrator(*boundary_integrators[i], scale_factor, false),
          *(*boundary_markers[i]));
  }

protected:
  mfem::ConstantCoefficient _dt_coef;

  Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
      _td_kernels_map;
  Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _td_blfs;
  Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
      _td_mblfs; // named according to trial variable

  std::vector<std::unique_ptr<mfem::ParGridFunction>> _td_var_ess_constraints;

  const Moose::MFEM::TimeDerivativeMap & _time_derivative_map;

private:
  virtual void SetTrialVariableNames() override;
};

} // namespace Moose::MFEM

#endif