Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #include "hcurl.hpp"
       5              : 
       6              : #include <mfem.hpp>
       7              : #include "fem/bilinearform.hpp"
       8              : #include "fem/fespace.hpp"
       9              : #include "fem/integrator.hpp"
      10              : #include "linalg/ams.hpp"
      11              : #include "linalg/gmg.hpp"
      12              : #include "linalg/iterative.hpp"
      13              : #include "linalg/rap.hpp"
      14              : #include "models/materialoperator.hpp"
      15              : 
      16              : namespace palace
      17              : {
      18              : 
      19              : namespace
      20              : {
      21              : 
      22              : template <typename OperType>
      23              : auto BuildLevelParOperator(std::unique_ptr<Operator> &&a,
      24              :                            const FiniteElementSpace &fespace);
      25              : 
      26              : template <>
      27              : auto BuildLevelParOperator<Operator>(std::unique_ptr<Operator> &&a,
      28              :                                      const FiniteElementSpace &fespace)
      29              : {
      30            0 :   return std::make_unique<ParOperator>(std::move(a), fespace);
      31              : }
      32              : 
      33              : template <>
      34              : auto BuildLevelParOperator<ComplexOperator>(std::unique_ptr<Operator> &&a,
      35              :                                             const FiniteElementSpace &fespace)
      36              : {
      37            0 :   return std::make_unique<ComplexParOperator>(std::move(a), nullptr, fespace);
      38              : }
      39              : 
      40              : }  // namespace
      41              : 
      42              : template <typename VecType>
      43            0 : WeightedHCurlNormSolver<VecType>::WeightedHCurlNormSolver(
      44              :     const MaterialOperator &mat_op, FiniteElementSpaceHierarchy &nd_fespaces,
      45              :     FiniteElementSpaceHierarchy &h1_fespaces,
      46              :     const std::vector<mfem::Array<int>> &nd_dbc_tdof_lists,
      47              :     const std::vector<mfem::Array<int>> &h1_dbc_tdof_lists, double tol, int max_it,
      48              :     int print)
      49              : {
      50            0 :   MFEM_VERIFY(h1_fespaces.GetNumLevels() == nd_fespaces.GetNumLevels(),
      51              :               "Multigrid hierarchy mismatch for auxiliary space preconditioning!");
      52            0 :   const auto n_levels = nd_fespaces.GetNumLevels();
      53              :   {
      54              :     constexpr bool skip_zeros = false;
      55            0 :     MaterialPropertyCoefficient muinv_func(mat_op.GetAttributeToMaterial(),
      56              :                                            mat_op.GetInvPermeability());
      57            0 :     MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
      58              :                                              mat_op.GetPermittivityReal());
      59              :     BilinearForm a(nd_fespaces.GetFinestFESpace()), a_aux(h1_fespaces.GetFinestFESpace());
      60            0 :     a.AddDomainIntegrator<CurlCurlMassIntegrator>(muinv_func, epsilon_func);
      61            0 :     a_aux.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
      62              :     // a.AssembleQuadratureData();
      63              :     // a_aux.AssembleQuadratureData();
      64            0 :     auto a_vec = a.Assemble(nd_fespaces, skip_zeros);
      65            0 :     auto a_aux_vec = a_aux.Assemble(h1_fespaces, skip_zeros);
      66            0 :     auto A_mg = std::make_unique<BaseMultigridOperator<OperType>>(n_levels);
      67            0 :     for (bool aux : {false, true})
      68              :     {
      69            0 :       for (std::size_t l = 0; l < n_levels; l++)
      70              :       {
      71            0 :         const auto &fespace_l =
      72              :             aux ? h1_fespaces.GetFESpaceAtLevel(l) : nd_fespaces.GetFESpaceAtLevel(l);
      73            0 :         const auto &dbc_tdof_lists_l = aux ? h1_dbc_tdof_lists[l] : nd_dbc_tdof_lists[l];
      74            0 :         auto A_l = BuildLevelParOperator<OperType>(std::move(aux ? a_aux_vec[l] : a_vec[l]),
      75              :                                                    fespace_l);
      76            0 :         A_l->SetEssentialTrueDofs(dbc_tdof_lists_l, Operator::DiagonalPolicy::DIAG_ONE);
      77            0 :         if (aux)
      78              :         {
      79            0 :           A_mg->AddAuxiliaryOperator(std::move(A_l));
      80              :         }
      81              :         else
      82              :         {
      83            0 :           A_mg->AddOperator(std::move(A_l));
      84              :         }
      85              :       }
      86              :     }
      87              :     A = std::move(A_mg);
      88            0 :   }
      89              : 
      90              :   // The system matrix K + M is real and SPD. We use Hypre's AMS solver as the coarse-level
      91              :   // multigrid solve.
      92            0 :   auto ams = std::make_unique<MfemWrapperSolver<OperType>>(std::make_unique<HypreAmsSolver>(
      93            0 :       nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetFESpaceAtLevel(0), 1, 1, false, true,
      94            0 :       false, 0));
      95            0 :   std::unique_ptr<Solver<OperType>> pc;
      96            0 :   if (n_levels > 1)
      97              :   {
      98            0 :     const auto G = nd_fespaces.GetDiscreteInterpolators(h1_fespaces);
      99            0 :     const int mg_smooth_order =
     100            0 :         std::max(nd_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2);
     101            0 :     pc = std::make_unique<GeometricMultigridSolver<OperType>>(
     102            0 :         nd_fespaces.GetFinestFESpace().GetComm(), std::move(ams),
     103            0 :         nd_fespaces.GetProlongationOperators(), &G, 1, 1, mg_smooth_order, 1.0, 0.0, true);
     104              :   }
     105              :   else
     106              :   {
     107              :     pc = std::move(ams);
     108              :   }
     109              : 
     110            0 :   auto pcg =
     111            0 :       std::make_unique<CgSolver<OperType>>(nd_fespaces.GetFinestFESpace().GetComm(), print);
     112            0 :   pcg->SetInitialGuess(false);
     113              :   pcg->SetRelTol(tol);
     114              :   pcg->SetMaxIter(max_it);
     115              : 
     116            0 :   ksp = std::make_unique<BaseKspSolver<OperType>>(std::move(pcg), std::move(pc));
     117            0 :   ksp->SetOperators(*A, *A);
     118            0 : }
     119              : 
     120              : template class WeightedHCurlNormSolver<Vector>;
     121              : template class WeightedHCurlNormSolver<ComplexVector>;
     122              : 
     123              : }  // namespace palace