Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #ifdef MOOSE_MFEM_ENABLED
      11             : 
      12             : #pragma once
      13             : #include "libmesh/ignore_warnings.h"
      14             : #include "mfem/miniapps/common/pfem_extras.hpp"
      15             : #include "libmesh/restore_warnings.h"
      16             : #include "MFEMIntegratedBC.h"
      17             : #include "MFEMEssentialBC.h"
      18             : #include "MFEMContainers.h"
      19             : #include "MFEMKernel.h"
      20             : #include "MFEMMixedBilinearFormKernel.h"
      21             : #include "ScaleIntegrator.h"
      22             : 
      23             : namespace Moose::MFEM
      24             : {
      25             : 
      26             : /**
      27             :  * Class to store weak form components (bilinear and linear forms, and optionally
      28             :  * mixed and nonlinear forms) and build methods
      29             :  */
      30             : class EquationSystem : public mfem::Operator
      31             : {
      32             : 
      33             : public:
      34             :   friend class EquationSystemProblemOperator;
      35             :   friend class TimeDomainEquationSystemProblemOperator;
      36             : 
      37         311 :   EquationSystem() = default;
      38             :   ~EquationSystem() override;
      39             : 
      40             :   /// Add test variable to EquationSystem.
      41             :   virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
      42             :   /// Add coupled variable to EquationSystem.
      43             :   virtual void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name);
      44             : 
      45             :   /// Add kernels.
      46             :   virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel);
      47             :   virtual void AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> kernel);
      48             :   virtual void AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc);
      49             : 
      50             :   /// Initialise
      51             :   virtual void Init(Moose::MFEM::GridFunctions & gridfunctions,
      52             :                     const Moose::MFEM::FESpaces & fespaces,
      53             :                     mfem::AssemblyLevel assembly_level);
      54             : 
      55             :   /// Build linear forms and eliminate constrained DoFs
      56             :   virtual void BuildLinearForms();
      57             :   virtual void ApplyEssentialBCs();
      58             :   virtual void EliminateCoupledVariables();
      59             : 
      60             :   /// Build bilinear forms
      61             :   virtual void BuildBilinearForms();
      62             :   virtual void BuildMixedBilinearForms();
      63             :   virtual void BuildEquationSystem();
      64             : 
      65             :   /// Form linear system, with essential boundary conditions accounted for
      66             :   virtual void FormLinearSystem(mfem::OperatorHandle & op,
      67             :                                 mfem::BlockVector & trueX,
      68             :                                 mfem::BlockVector & trueRHS);
      69             : 
      70             :   virtual void
      71             :   FormSystem(mfem::OperatorHandle & op, mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
      72             :   virtual void FormLegacySystem(mfem::OperatorHandle & op,
      73             :                                 mfem::BlockVector & trueX,
      74             :                                 mfem::BlockVector & trueRHS);
      75             : 
      76             :   /// Build linear system, with essential boundary conditions accounted for
      77             :   virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
      78             : 
      79             :   /// Compute residual y = Mu
      80             :   void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
      81             : 
      82             :   /// Compute J = M + grad_H(u)
      83             :   mfem::Operator & GetGradient(const mfem::Vector & u) const override;
      84             : 
      85             :   /// Update variable from solution vector after solve
      86             :   virtual void RecoverFEMSolution(mfem::BlockVector & trueX,
      87             :                                   Moose::MFEM::GridFunctions & gridfunctions);
      88             : 
      89             :   std::vector<mfem::Array<int>> _ess_tdof_lists;
      90             : 
      91         311 :   const std::vector<std::string> & TrialVarNames() const { return _trial_var_names; }
      92         311 :   const std::vector<std::string> & TestVarNames() const { return _test_var_names; }
      93             : 
      94             : private:
      95             :   /// Disallowed inherited method
      96             :   using mfem::Operator::RecoverFEMSolution;
      97             : 
      98             :   /// Set trial variable names from subset of coupled variables that have an associated test variable.
      99             :   virtual void SetTrialVariableNames();
     100             : 
     101             : protected:
     102             :   /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
     103             :   /// and then proceeds to delete all dynamically allocated memory for _h_blocks
     104             :   /// itself, resetting all dimensions to zero.
     105             :   void DeleteAllBlocks();
     106             : 
     107             :   bool VectorContainsName(const std::vector<std::string> & the_vector,
     108             :                           const std::string & name) const;
     109             : 
     110             :   // Test variables are associated with LinearForms,
     111             :   // whereas trial variables are associated with gridfunctions.
     112             : 
     113             :   /// Names of all trial variables of kernels and boundary conditions
     114             :   /// added to this EquationSystem.
     115             :   std::vector<std::string> _coupled_var_names;
     116             :   /// Subset of _coupled_var_names of all variables corresponding to gridfunctions with degrees of
     117             :   /// freedom that comprise the state vector of this EquationSystem. This will differ from
     118             :   /// _coupled_var_names when time derivatives or other eliminated variables are present.
     119             :   std::vector<std::string> _trial_var_names;
     120             :   /// Names of all coupled variables without a corresponding test variable.
     121             :   std::vector<std::string> _eliminated_var_names;
     122             :   /// Pointers to coupled variables not part of the reduced EquationSystem.
     123             :   Moose::MFEM::GridFunctions _eliminated_variables;
     124             :   /// Names of all test variables corresponding to linear forms in this equation system
     125             :   std::vector<std::string> _test_var_names;
     126             :   /// Pointers to finite element spaces associated with test variables.
     127             :   std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
     128             :   /// Pointers to finite element spaces associated with coupled variables.
     129             :   std::vector<mfem::ParFiniteElementSpace *> _coupled_pfespaces;
     130             : 
     131             :   // Components of weak form. // Named according to test variable
     132             :   Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
     133             :   Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _lfs;
     134             :   Moose::MFEM::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
     135             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
     136             :       _mblfs; // named according to trial variable
     137             : 
     138             :   /**
     139             :    * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
     140             :    * or MixedBilinearForm
     141             :    */
     142             :   template <class FormType>
     143             :   void ApplyDomainBLFIntegrators(
     144             :       const std::string & trial_var_name,
     145             :       const std::string & test_var_name,
     146             :       std::shared_ptr<FormType> form,
     147             :       Moose::MFEM::NamedFieldsMap<
     148             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
     149             : 
     150             :   void ApplyDomainLFIntegrators(
     151             :       const std::string & test_var_name,
     152             :       std::shared_ptr<mfem::ParLinearForm> form,
     153             :       Moose::MFEM::NamedFieldsMap<
     154             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
     155             : 
     156             :   template <class FormType>
     157             :   void ApplyBoundaryBLFIntegrators(
     158             :       const std::string & trial_var_name,
     159             :       const std::string & test_var_name,
     160             :       std::shared_ptr<FormType> form,
     161             :       Moose::MFEM::NamedFieldsMap<
     162             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     163             :           integrated_bc_map);
     164             : 
     165             :   void ApplyBoundaryLFIntegrators(
     166             :       const std::string & test_var_name,
     167             :       std::shared_ptr<mfem::ParLinearForm> form,
     168             :       Moose::MFEM::NamedFieldsMap<
     169             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     170             :           integrated_bc_map);
     171             : 
     172             :   /// Gridfunctions holding essential constraints from Dirichlet BCs
     173             :   std::vector<std::unique_ptr<mfem::ParGridFunction>> _var_ess_constraints;
     174             : 
     175             :   mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks;
     176             : 
     177             :   /// Arrays to store kernels to act on each component of weak form.
     178             :   /// Named according to test and trial variables.
     179             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
     180             :       _kernels_map;
     181             :   /// Arrays to store integrated BCs to act on each component of weak form.
     182             :   /// Named according to test and trial variables.
     183             :   Moose::MFEM::NamedFieldsMap<
     184             :       Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>
     185             :       _integrated_bc_map;
     186             :   /// Arrays to store essential BCs to act on each component of weak form.
     187             :   /// Named according to test variable.
     188             :   Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
     189             : 
     190             :   mutable mfem::OperatorHandle _jacobian;
     191             : 
     192             :   mfem::AssemblyLevel _assembly_level;
     193             : };
     194             : 
     195             : template <class FormType>
     196             : void
     197        1624 : EquationSystem::ApplyDomainBLFIntegrators(
     198             :     const std::string & trial_var_name,
     199             :     const std::string & test_var_name,
     200             :     std::shared_ptr<FormType> form,
     201             :     Moose::MFEM::NamedFieldsMap<
     202             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
     203             : {
     204        1624 :   if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
     205             :   {
     206        1602 :     auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
     207        3285 :     for (auto & kernel : kernels)
     208             :     {
     209        1683 :       mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
     210        1683 :       if (integ != nullptr)
     211             :       {
     212        1632 :         kernel->isSubdomainRestricted()
     213        1632 :             ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
     214        1627 :             : form->AddDomainIntegrator(std::move(integ));
     215             :       }
     216             :     }
     217        1602 :   }
     218        1624 : }
     219             : 
     220             : inline void
     221         884 : EquationSystem::ApplyDomainLFIntegrators(
     222             :     const std::string & test_var_name,
     223             :     std::shared_ptr<mfem::ParLinearForm> form,
     224             :     Moose::MFEM::NamedFieldsMap<
     225             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
     226             : {
     227         884 :   if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
     228             :   {
     229         876 :     auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
     230        1833 :     for (auto & kernel : kernels)
     231             :     {
     232         957 :       mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
     233         957 :       if (integ != nullptr)
     234             :       {
     235          51 :         kernel->isSubdomainRestricted()
     236          51 :             ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
     237          51 :             : form->AddDomainIntegrator(std::move(integ));
     238             :       }
     239             :     }
     240         876 :   }
     241         884 : }
     242             : 
     243             : template <class FormType>
     244             : void
     245        1603 : EquationSystem::ApplyBoundaryBLFIntegrators(
     246             :     const std::string & trial_var_name,
     247             :     const std::string & test_var_name,
     248             :     std::shared_ptr<FormType> form,
     249             :     Moose::MFEM::NamedFieldsMap<
     250             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     251             :         integrated_bc_map)
     252             : {
     253        1731 :   if (integrated_bc_map.Has(test_var_name) &&
     254         128 :       integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
     255             :   {
     256         128 :     auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
     257         256 :     for (auto & bc : bcs)
     258             :     {
     259         128 :       mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
     260         128 :       if (integ != nullptr)
     261             :       {
     262         112 :         bc->isBoundaryRestricted()
     263         112 :             ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
     264           0 :             : form->AddBoundaryIntegrator(std::move(integ));
     265             :       }
     266             :     }
     267         128 :   }
     268        1603 : }
     269             : 
     270             : inline void
     271         884 : EquationSystem::ApplyBoundaryLFIntegrators(
     272             :     const std::string & test_var_name,
     273             :     std::shared_ptr<mfem::ParLinearForm> form,
     274             :     Moose::MFEM::NamedFieldsMap<
     275             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     276             :         integrated_bc_map)
     277             : {
     278         956 :   if (integrated_bc_map.Has(test_var_name) &&
     279          72 :       integrated_bc_map.Get(test_var_name)->Has(test_var_name))
     280             :   {
     281          72 :     auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
     282         144 :     for (auto & bc : bcs)
     283             :     {
     284          72 :       mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
     285          72 :       if (integ != nullptr)
     286             :       {
     287          72 :         bc->isBoundaryRestricted()
     288          72 :             ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
     289           8 :             : form->AddBoundaryIntegrator(std::move(integ));
     290             :       }
     291             :     }
     292          72 :   }
     293         884 : }
     294             : 
     295             : /**
     296             :  * Class to store weak form components for time dependent PDEs
     297             :  */
     298             : class TimeDependentEquationSystem : public EquationSystem
     299             : {
     300             : public:
     301             :   TimeDependentEquationSystem();
     302             : 
     303             :   void AddCoupledVariableNameIfMissing(const std::string & coupled_var_name) override;
     304             : 
     305             :   virtual void SetTimeStep(double dt);
     306             :   virtual void UpdateEquationSystem();
     307             : 
     308             :   virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel) override;
     309             :   virtual void BuildBilinearForms() override;
     310             :   virtual void FormLegacySystem(mfem::OperatorHandle & op,
     311             :                                 mfem::BlockVector & truedXdt,
     312             :                                 mfem::BlockVector & trueRHS) override;
     313             :   virtual void FormSystem(mfem::OperatorHandle & op,
     314             :                           mfem::BlockVector & truedXdt,
     315             :                           mfem::BlockVector & trueRHS) override;
     316             : 
     317             : protected:
     318             :   /// Coefficient for timestep scaling
     319             :   mfem::ConstantCoefficient _dt_coef;
     320             : 
     321             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
     322             :       _td_kernels_map;
     323             :   /// Container to store contributions to weak form of the form (F du/dt, v)
     324             :   Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _td_blfs;
     325             : 
     326             : private:
     327             :   /// Set trial variable names from subset of coupled variables that have an associated test variable.
     328             :   virtual void SetTrialVariableNames() override;
     329             : };
     330             : 
     331             : } // namespace Moose::MFEM
     332             : 
     333             : #endif