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