Line data    Source code

       1             : #ifdef MOOSE_MFEM_ENABLED
       2             : 
       3             : #include "ComplexEquationSystem.h"
       4             : #include "libmesh/int_range.h"
       5             : 
       6             : namespace Moose::MFEM
       7             : {
       8             : 
       9             : void
      10          63 : ComplexEquationSystem::Init(GridFunctions & gridfunctions,
      11             :                             ComplexGridFunctions & cmplx_gridfunctions,
      12             :                             mfem::AssemblyLevel assembly_level)
      13             : {
      14          63 :   _assembly_level = assembly_level;
      15             : 
      16          63 :   if (gridfunctions.size())
      17           0 :     mooseError("Mixing real and complex variables is currently not supported.");
      18             : 
      19         105 :   for (auto & test_var_name : _test_var_names)
      20             :   {
      21          42 :     if (!cmplx_gridfunctions.Has(test_var_name))
      22             :     {
      23           0 :       mooseError("MFEM complex variable ",
      24             :                  test_var_name,
      25             :                  " requested by equation system during initialization was "
      26             :                  "not found in gridfunctions");
      27             :     }
      28             :     // Store pointers to test FESpaces
      29          42 :     _test_pfespaces.push_back(cmplx_gridfunctions.Get(test_var_name)->ParFESpace());
      30             :     // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
      31          42 :     _cmplx_var_ess_constraints.emplace_back(std::make_unique<mfem::ParComplexGridFunction>(
      32          84 :         cmplx_gridfunctions.Get(test_var_name)->ParFESpace()));
      33             :   }
      34             : 
      35             :   // Store pointers to FESpaces of all coupled variables
      36         105 :   for (auto & coupled_var_name : _coupled_var_names)
      37          42 :     _coupled_pfespaces.push_back(cmplx_gridfunctions.Get(coupled_var_name)->ParFESpace());
      38             : 
      39             :   // Extract which coupled variables are to be trivially eliminated and which are trial variables
      40          63 :   SetTrialVariableNames();
      41             : 
      42             :   // Store pointers to coupled variable ComplexGridFunctions that are to be eliminated prior to
      43             :   // forming the jacobian
      44          63 :   for (auto & eliminated_var_name : _eliminated_var_names)
      45           0 :     _cmplx_eliminated_variables.Register(eliminated_var_name,
      46           0 :                                          cmplx_gridfunctions.GetShared(eliminated_var_name));
      47             : 
      48             :   // Get a reference to the complex GridFunctions
      49          63 :   _complex_gfuncs = &cmplx_gridfunctions;
      50          63 : }
      51             : 
      52             : void
      53          42 : ComplexEquationSystem::BuildEquationSystem()
      54             : {
      55          42 :   BuildBilinearForms();
      56          42 :   BuildLinearForms();
      57          42 : }
      58             : 
      59             : void
      60          42 : ComplexEquationSystem::BuildLinearForms()
      61             : {
      62             :   // Register linear forms
      63          84 :   for (const auto i : index_range(_test_var_names))
      64             :   {
      65          42 :     auto test_var_name = _test_var_names.at(i);
      66          42 :     _clfs.Register(test_var_name,
      67          84 :                    std::make_shared<mfem::ParComplexLinearForm>(_test_pfespaces.at(i)));
      68          42 :     _clfs.GetRef(test_var_name) = 0.0;
      69          42 :   }
      70             :   // Apply boundary conditions
      71          42 :   ApplyEssentialBCs();
      72             : 
      73          84 :   for (auto & test_var_name : _test_var_names)
      74             :   {
      75             :     // Apply kernels
      76          42 :     auto clf = _clfs.GetShared(test_var_name);
      77          42 :     ApplyDomainLFIntegrators(test_var_name, clf, _cmplx_kernels_map);
      78          42 :     ApplyBoundaryLFIntegrators(test_var_name, clf, _cmplx_integrated_bc_map);
      79          42 :     clf->Assemble();
      80          42 :   }
      81          42 : }
      82             : 
      83             : void
      84          42 : ComplexEquationSystem::BuildBilinearForms()
      85             : {
      86             :   // Register bilinear forms
      87          84 :   for (const auto i : index_range(_test_var_names))
      88             :   {
      89          42 :     auto test_var_name = _test_var_names.at(i);
      90          42 :     _slfs.Register(test_var_name,
      91          84 :                    std::make_shared<mfem::ParSesquilinearForm>(_test_pfespaces.at(i)));
      92             : 
      93             :     // Apply kernels
      94          42 :     auto slf = _slfs.GetShared(test_var_name);
      95          42 :     slf->SetAssemblyLevel(_assembly_level);
      96          42 :     ApplyBoundaryBLFIntegrators<mfem::ParSesquilinearForm>(
      97          42 :         test_var_name, test_var_name, slf, _cmplx_integrated_bc_map);
      98          42 :     ApplyDomainBLFIntegrators<mfem::ParSesquilinearForm>(
      99          42 :         test_var_name, test_var_name, slf, _cmplx_kernels_map);
     100             :     // Assemble
     101          42 :     slf->Assemble();
     102          42 :   }
     103          42 : }
     104             : 
     105             : void
     106          42 : ComplexEquationSystem::ApplyComplexEssentialBC(const std::string & var_name,
     107             :                                                mfem::ParComplexGridFunction & trial_gf,
     108             :                                                mfem::Array<int> & global_ess_markers)
     109             : {
     110          42 :   if (_cmplx_essential_bc_map.Has(var_name))
     111             :   {
     112          42 :     auto & bcs = _cmplx_essential_bc_map.GetRef(var_name);
     113         105 :     for (auto & bc : bcs)
     114             :     {
     115             :       // Set constrained DoFs values on essential boundaries
     116          63 :       bc->ApplyBC(trial_gf);
     117             :       // Fetch marker array labelling essential boundaries of current BC
     118          63 :       mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
     119             :       // Add these boundary markers to the set of markers labelling all essential boundaries
     120         392 :       for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
     121         329 :         global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
     122          63 :     }
     123             :   }
     124          42 : }
     125             : 
     126             : void
     127          42 : ComplexEquationSystem::ApplyEssentialBCs()
     128             : {
     129          42 :   _ess_tdof_lists.resize(_trial_var_names.size());
     130          84 :   for (const auto i : index_range(_trial_var_names))
     131             :   {
     132          42 :     const auto & trial_var_name = _trial_var_names.at(i);
     133          42 :     mfem::ParComplexGridFunction & trial_gf = *_cmplx_var_ess_constraints.at(i);
     134             : 
     135             :     // Make sure we update the size, if this mesh has changed recently for instance
     136          42 :     trial_gf.Update();
     137             : 
     138             :     // Initial guess for iterative solvers (initial condition or the previous time step solution)
     139          42 :     static_cast<mfem::Vector &>(trial_gf) = _complex_gfuncs->GetRef(trial_var_name);
     140             : 
     141          42 :     mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
     142          42 :     global_ess_markers = 0;
     143             :     // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
     144             :     // essential boundaries to the global_ess_markers array
     145          42 :     ApplyComplexEssentialBC(trial_var_name, trial_gf, global_ess_markers);
     146          42 :     trial_gf.FESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
     147          42 :   }
     148          42 : }
     149             : 
     150             : void
     151          63 : ComplexEquationSystem::AddComplexKernel(std::shared_ptr<MFEMComplexKernel> kernel)
     152             : {
     153          63 :   const auto & trial_var_name = kernel->getTrialVariableName();
     154          63 :   const auto & test_var_name = kernel->getTestVariableName();
     155          63 :   AddCoupledVariableNameIfMissing(trial_var_name);
     156          63 :   AddTestVariableNameIfMissing(test_var_name);
     157             :   // Register new complex kernels map if not present for the test variable
     158          63 :   if (!_cmplx_kernels_map.Has(test_var_name))
     159             :   {
     160             :     auto kernel_field_map =
     161          42 :         std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexKernel>>>>();
     162          42 :     _cmplx_kernels_map.Register(test_var_name, std::move(kernel_field_map));
     163          42 :   }
     164             :   // Register new complex kernels map if not present for the test/trial variable pair
     165          63 :   if (!_cmplx_kernels_map.Get(test_var_name)->Has(trial_var_name))
     166             :   {
     167          42 :     auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMComplexKernel>>>();
     168          42 :     _cmplx_kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
     169          42 :   }
     170          63 :   _cmplx_kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
     171          63 : }
     172             : 
     173             : void
     174          14 : ComplexEquationSystem::AddComplexIntegratedBC(std::shared_ptr<MFEMComplexIntegratedBC> bc)
     175             : {
     176          14 :   const auto & trial_var_name = bc->getTrialVariableName();
     177          14 :   const auto & test_var_name = bc->getTestVariableName();
     178          14 :   AddCoupledVariableNameIfMissing(trial_var_name);
     179          14 :   AddTestVariableNameIfMissing(test_var_name);
     180             :   // Register new complex integrated bc map if not present for the test variable
     181          14 :   if (!_cmplx_integrated_bc_map.Has(test_var_name))
     182             :   {
     183             :     auto integrated_bc_field_map =
     184           7 :         std::make_shared<NamedFieldsMap<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>>();
     185           7 :     _cmplx_integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
     186           7 :   }
     187             :   // Register new complex integrated bc map if not present for the test/trial variable pair
     188          14 :   if (!_cmplx_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
     189             :   {
     190           7 :     auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexIntegratedBC>>>();
     191           7 :     _cmplx_integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
     192           7 :   }
     193          14 :   _cmplx_integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
     194          14 : }
     195             : 
     196             : void
     197          63 : ComplexEquationSystem::AddComplexEssentialBCs(std::shared_ptr<MFEMComplexEssentialBC> bc)
     198             : {
     199          63 :   const auto & test_var_name = bc->getTestVariableName();
     200          63 :   AddTestVariableNameIfMissing(test_var_name);
     201             :   // Register new complex essential bc map if not present for the test variable
     202          63 :   if (!_cmplx_essential_bc_map.Has(test_var_name))
     203             :   {
     204          42 :     auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMComplexEssentialBC>>>();
     205          42 :     _cmplx_essential_bc_map.Register(test_var_name, std::move(bcs));
     206          42 :   }
     207          63 :   _cmplx_essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
     208          63 : }
     209             : 
     210             : void
     211           0 : ComplexEquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
     212             :                                           mfem::BlockVector & trueX,
     213             :                                           mfem::BlockVector & trueRHS)
     214             : {
     215           0 :   auto & test_var_name = _test_var_names.at(0);
     216           0 :   mfem::Vector aux_x, aux_rhs;
     217           0 :   mfem::OperatorPtr aux_a;
     218             : 
     219           0 :   auto slf = _slfs.Get(test_var_name);
     220           0 :   slf->FormLinearSystem(_ess_tdof_lists.at(0),
     221           0 :                         *_cmplx_var_ess_constraints.at(0),
     222           0 :                         *_clfs.Get(test_var_name),
     223             :                         aux_a,
     224             :                         aux_x,
     225             :                         aux_rhs,
     226             :                         /*copy_interior=*/true);
     227             : 
     228           0 :   trueX.GetBlock(0) = aux_x;
     229           0 :   trueRHS.GetBlock(0) = aux_rhs;
     230           0 :   trueX.SyncFromBlocks();
     231           0 :   trueRHS.SyncFromBlocks();
     232             : 
     233           0 :   op.Reset(aux_a.Ptr());
     234           0 :   aux_a.SetOperatorOwner(false);
     235           0 : }
     236             : 
     237             : void
     238          42 : ComplexEquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
     239             :                                         mfem::BlockVector & trueX,
     240             :                                         mfem::BlockVector & trueRHS)
     241             : {
     242             : 
     243             :   // Allocate block operator
     244          42 :   DeleteHBlocks();
     245          42 :   _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
     246          42 :   _h_blocks = nullptr;
     247             :   // Zero out RHS and sync memory
     248          42 :   trueRHS = 0.0;
     249          42 :   trueRHS.SyncToBlocks();
     250             : 
     251             :   // Form diagonal blocks.
     252          84 :   for (const auto i : index_range(_test_var_names))
     253             :   {
     254          42 :     auto & test_var_name = _test_var_names.at(i);
     255             : 
     256          42 :     mfem::Vector aux_x, aux_rhs;
     257          42 :     mfem::OperatorHandle aux_a;
     258             : 
     259          42 :     auto slf = _slfs.Get(test_var_name);
     260          42 :     slf->FormLinearSystem(_ess_tdof_lists.at(i),
     261          42 :                           *_cmplx_var_ess_constraints.at(i),
     262          42 :                           *_clfs.Get(test_var_name),
     263             :                           aux_a,
     264             :                           aux_x,
     265             :                           aux_rhs,
     266             :                           /*copy_interior=*/true);
     267          42 :     trueX.GetBlock(i) = aux_x;
     268          42 :     trueRHS.GetBlock(i) = aux_rhs;
     269          42 :     _h_blocks(i, i) = aux_a.As<mfem::ComplexHypreParMatrix>()->GetSystemMatrix();
     270          42 :   }
     271             :   // Sync memory
     272          42 :   trueX.SyncFromBlocks();
     273          42 :   trueRHS.SyncFromBlocks();
     274             : 
     275             :   // Create monolithic matrix
     276          42 :   op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
     277          42 : }
     278             : 
     279             : // Equation system Mult
     280             : void
     281           0 : ComplexEquationSystem::Mult(const mfem::Vector & x, mfem::Vector & residual) const
     282             : {
     283           0 :   _linear_operator->Mult(x, residual);
     284           0 :   x.HostRead();
     285           0 :   residual.HostRead();
     286           0 : }
     287             : 
     288             : void
     289          42 : ComplexEquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
     290             : {
     291          84 :   for (const auto i : index_range(_trial_var_names))
     292             :   {
     293          42 :     auto & trial_var_name = _trial_var_names.at(i);
     294          42 :     trueX.GetBlock(i).SyncMemory(trueX);
     295          42 :     _complex_gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
     296             :   }
     297          42 : }
     298             : 
     299             : }
     300             : 
     301             : #endif