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