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         236 :   EquationSystem() = default;
      38             :   ~EquationSystem() override;
      39             : 
      40             :   /// Add test variable to EquationSystem.
      41             :   virtual void AddTestVariableNameIfMissing(const std::string & test_var_name);
      42             :   /// Add trial variable to EquationSystem.
      43             :   virtual void AddTrialVariableNameIfMissing(const std::string & trial_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             :   virtual void ApplyEssentialBCs();
      50             : 
      51             :   /// Build forms
      52             :   virtual void Init(Moose::MFEM::GridFunctions & gridfunctions,
      53             :                     const Moose::MFEM::FESpaces & fespaces,
      54             :                     mfem::AssemblyLevel assembly_level);
      55             :   virtual void BuildLinearForms();
      56             :   virtual void BuildBilinearForms();
      57             :   virtual void BuildMixedBilinearForms();
      58             :   virtual void BuildEquationSystem();
      59             : 
      60             :   /// Form linear system, with essential boundary conditions accounted for
      61             :   virtual void FormLinearSystem(mfem::OperatorHandle & op,
      62             :                                 mfem::BlockVector & trueX,
      63             :                                 mfem::BlockVector & trueRHS);
      64             : 
      65             :   virtual void
      66             :   FormSystem(mfem::OperatorHandle & op, mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
      67             :   virtual void FormLegacySystem(mfem::OperatorHandle & op,
      68             :                                 mfem::BlockVector & trueX,
      69             :                                 mfem::BlockVector & trueRHS);
      70             : 
      71             :   /// Build linear system, with essential boundary conditions accounted for
      72             :   virtual void BuildJacobian(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS);
      73             : 
      74             :   /// Compute residual y = Mu
      75             :   void Mult(const mfem::Vector & u, mfem::Vector & residual) const override;
      76             : 
      77             :   /// Compute J = M + grad_H(u)
      78             :   mfem::Operator & GetGradient(const mfem::Vector & u) const override;
      79             : 
      80             :   /// Update variable from solution vector after solve
      81             :   virtual void RecoverFEMSolution(mfem::BlockVector & trueX,
      82             :                                   Moose::MFEM::GridFunctions & gridfunctions);
      83             : 
      84             :   std::vector<mfem::Array<int>> _ess_tdof_lists;
      85             : 
      86         236 :   const std::vector<std::string> & TrialVarNames() const { return _trial_var_names; }
      87         236 :   const std::vector<std::string> & TestVarNames() const { return _test_var_names; }
      88             : 
      89             : private:
      90             :   /// Disallowed inherited method
      91             :   using mfem::Operator::RecoverFEMSolution;
      92             : 
      93             : protected:
      94             :   /// Deletes the HypreParMatrix associated with any pointer stored in _h_blocks,
      95             :   /// and then proceeds to delete all dynamically allocated memory for _h_blocks
      96             :   /// itself, resetting all dimensions to zero.
      97             :   void DeleteAllBlocks();
      98             : 
      99             :   bool VectorContainsName(const std::vector<std::string> & the_vector,
     100             :                           const std::string & name) const;
     101             : 
     102             :   // Test variables are associated with LinearForms,
     103             :   // whereas trial variables are associated with gridfunctions.
     104             : 
     105             :   /// Names of all variables corresponding to gridfunctions. This may differ
     106             :   /// from test_var_names when time derivatives are present.
     107             :   std::vector<std::string> _trial_var_names;
     108             :   /// Pointers to trial variables.
     109             :   Moose::MFEM::GridFunctions _trial_variables;
     110             :   /// Names of all test variables corresponding to linear forms in this equation system
     111             :   std::vector<std::string> _test_var_names;
     112             :   /// Pointers to finite element spaces associated with test variables.
     113             :   std::vector<mfem::ParFiniteElementSpace *> _test_pfespaces;
     114             : 
     115             :   // Components of weak form. // Named according to test variable
     116             :   Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _blfs;
     117             :   Moose::MFEM::NamedFieldsMap<mfem::ParLinearForm> _lfs;
     118             :   Moose::MFEM::NamedFieldsMap<mfem::ParNonlinearForm> _nlfs;
     119             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>
     120             :       _mblfs; // named according to trial variable
     121             : 
     122             :   /**
     123             :    * Template method for applying BilinearFormIntegrators on domains from kernels to a BilinearForm,
     124             :    * or MixedBilinearForm
     125             :    */
     126             :   template <class FormType>
     127             :   void ApplyDomainBLFIntegrators(
     128             :       const std::string & trial_var_name,
     129             :       const std::string & test_var_name,
     130             :       std::shared_ptr<FormType> form,
     131             :       Moose::MFEM::NamedFieldsMap<
     132             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
     133             : 
     134             :   void ApplyDomainLFIntegrators(
     135             :       const std::string & test_var_name,
     136             :       std::shared_ptr<mfem::ParLinearForm> form,
     137             :       Moose::MFEM::NamedFieldsMap<
     138             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map);
     139             : 
     140             :   template <class FormType>
     141             :   void ApplyBoundaryBLFIntegrators(
     142             :       const std::string & trial_var_name,
     143             :       const std::string & test_var_name,
     144             :       std::shared_ptr<FormType> form,
     145             :       Moose::MFEM::NamedFieldsMap<
     146             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     147             :           integrated_bc_map);
     148             : 
     149             :   void ApplyBoundaryLFIntegrators(
     150             :       const std::string & test_var_name,
     151             :       std::shared_ptr<mfem::ParLinearForm> form,
     152             :       Moose::MFEM::NamedFieldsMap<
     153             :           Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     154             :           integrated_bc_map);
     155             : 
     156             :   /// Gridfunctions for setting Dirichlet BCs
     157             :   std::vector<std::unique_ptr<mfem::ParGridFunction>> _xs;
     158             :   std::vector<std::unique_ptr<mfem::ParGridFunction>> _dxdts;
     159             : 
     160             :   mfem::Array2D<const mfem::HypreParMatrix *> _h_blocks;
     161             : 
     162             :   /// Arrays to store kernels to act on each component of weak form.
     163             :   /// Named according to test and trial variables.
     164             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
     165             :       _kernels_map;
     166             :   /// Arrays to store integrated BCs to act on each component of weak form.
     167             :   /// Named according to test and trial variables.
     168             :   Moose::MFEM::NamedFieldsMap<
     169             :       Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>
     170             :       _integrated_bc_map;
     171             :   /// Arrays to store essential BCs to act on each component of weak form.
     172             :   /// Named according to test variable.
     173             :   Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMEssentialBC>>> _essential_bc_map;
     174             : 
     175             :   mutable mfem::OperatorHandle _jacobian;
     176             : 
     177             :   mfem::AssemblyLevel _assembly_level;
     178             : };
     179             : 
     180             : template <class FormType>
     181             : void
     182         541 : EquationSystem::ApplyDomainBLFIntegrators(
     183             :     const std::string & trial_var_name,
     184             :     const std::string & test_var_name,
     185             :     std::shared_ptr<FormType> form,
     186             :     Moose::MFEM::NamedFieldsMap<
     187             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
     188             : {
     189         541 :   if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(trial_var_name))
     190             :   {
     191         517 :     auto kernels = kernels_map.GetRef(test_var_name).GetRef(trial_var_name);
     192        1107 :     for (auto & kernel : kernels)
     193             :     {
     194         590 :       mfem::BilinearFormIntegrator * integ = kernel->createBFIntegrator();
     195         590 :       if (integ != nullptr)
     196             :       {
     197         547 :         kernel->isSubdomainRestricted()
     198         547 :             ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
     199         547 :             : form->AddDomainIntegrator(std::move(integ));
     200             :       }
     201             :     }
     202         517 :   }
     203         541 : }
     204             : 
     205             : inline void
     206         336 : EquationSystem::ApplyDomainLFIntegrators(
     207             :     const std::string & test_var_name,
     208             :     std::shared_ptr<mfem::ParLinearForm> form,
     209             :     Moose::MFEM::NamedFieldsMap<
     210             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
     211             : {
     212         336 :   if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
     213             :   {
     214         328 :     auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
     215         729 :     for (auto & kernel : kernels)
     216             :     {
     217         401 :       mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();
     218         401 :       if (integ != nullptr)
     219             :       {
     220          43 :         kernel->isSubdomainRestricted()
     221          43 :             ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
     222          43 :             : form->AddDomainIntegrator(std::move(integ));
     223             :       }
     224             :     }
     225         328 :   }
     226         336 : }
     227             : 
     228             : template <class FormType>
     229             : void
     230         525 : EquationSystem::ApplyBoundaryBLFIntegrators(
     231             :     const std::string & trial_var_name,
     232             :     const std::string & test_var_name,
     233             :     std::shared_ptr<FormType> form,
     234             :     Moose::MFEM::NamedFieldsMap<
     235             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     236             :         integrated_bc_map)
     237             : {
     238         605 :   if (integrated_bc_map.Has(test_var_name) &&
     239          80 :       integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
     240             :   {
     241          80 :     auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(trial_var_name);
     242         160 :     for (auto & bc : bcs)
     243             :     {
     244          80 :       mfem::BilinearFormIntegrator * integ = bc->createBFIntegrator();
     245          80 :       if (integ != nullptr)
     246             :       {
     247          64 :         bc->isBoundaryRestricted()
     248          64 :             ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
     249           0 :             : form->AddBoundaryIntegrator(std::move(integ));
     250             :       }
     251             :     }
     252          80 :   }
     253         525 : }
     254             : 
     255             : inline void
     256         336 : EquationSystem::ApplyBoundaryLFIntegrators(
     257             :     const std::string & test_var_name,
     258             :     std::shared_ptr<mfem::ParLinearForm> form,
     259             :     Moose::MFEM::NamedFieldsMap<
     260             :         Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
     261             :         integrated_bc_map)
     262             : {
     263         384 :   if (integrated_bc_map.Has(test_var_name) &&
     264          48 :       integrated_bc_map.Get(test_var_name)->Has(test_var_name))
     265             :   {
     266          48 :     auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
     267          96 :     for (auto & bc : bcs)
     268             :     {
     269          48 :       mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();
     270          48 :       if (integ != nullptr)
     271             :       {
     272          48 :         bc->isBoundaryRestricted()
     273          48 :             ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
     274           8 :             : form->AddBoundaryIntegrator(std::move(integ));
     275             :       }
     276             :     }
     277          48 :   }
     278         336 : }
     279             : 
     280             : /**
     281             :  * Class to store weak form components for time dependent PDEs
     282             :  */
     283             : class TimeDependentEquationSystem : public EquationSystem
     284             : {
     285             : public:
     286             :   TimeDependentEquationSystem();
     287             : 
     288             :   void AddTrialVariableNameIfMissing(const std::string & trial_var_name) override;
     289             : 
     290             :   virtual void SetTimeStep(double dt);
     291             :   virtual void UpdateEquationSystem();
     292             : 
     293             :   virtual void AddKernel(std::shared_ptr<MFEMKernel> kernel) override;
     294             :   virtual void BuildBilinearForms() override;
     295             :   virtual void FormLegacySystem(mfem::OperatorHandle & op,
     296             :                                 mfem::BlockVector & truedXdt,
     297             :                                 mfem::BlockVector & trueRHS) override;
     298             :   virtual void FormSystem(mfem::OperatorHandle & op,
     299             :                           mfem::BlockVector & truedXdt,
     300             :                           mfem::BlockVector & trueRHS) override;
     301             : 
     302             :   const std::vector<std::string> & TrialVarTimeDerivativeNames() const;
     303             : 
     304             : protected:
     305             :   /// Coefficient for timestep scaling
     306             :   mfem::ConstantCoefficient _dt_coef;
     307             :   std::vector<std::string> _trial_var_time_derivative_names;
     308             : 
     309             :   Moose::MFEM::NamedFieldsMap<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>
     310             :       _td_kernels_map;
     311             :   /// Container to store contributions to weak form of the form (F du/dt, v)
     312             :   Moose::MFEM::NamedFieldsMap<mfem::ParBilinearForm> _td_blfs;
     313             : };
     314             : 
     315             : inline const std::vector<std::string> &
     316             : TimeDependentEquationSystem::TrialVarTimeDerivativeNames() const
     317             : {
     318             :   return _trial_var_time_derivative_names;
     319             : }
     320             : 
     321             : } // namespace Moose::MFEM
     322             : 
     323             : #endif