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 "divfree.hpp"
5 :
6 : #include <limits>
7 : #include <mfem.hpp>
8 : #include "fem/bilinearform.hpp"
9 : #include "fem/fespace.hpp"
10 : #include "fem/integrator.hpp"
11 : #include "linalg/amg.hpp"
12 : #include "linalg/gmg.hpp"
13 : #include "linalg/iterative.hpp"
14 : #include "linalg/rap.hpp"
15 : #include "models/materialoperator.hpp"
16 : #include "utils/timer.hpp"
17 :
18 : namespace palace
19 : {
20 :
21 : namespace
22 : {
23 :
24 : template <typename OperType>
25 : auto BuildLevelParOperator(std::unique_ptr<Operator> &&a,
26 : const FiniteElementSpace &fespace);
27 :
28 : template <>
29 : auto BuildLevelParOperator<Operator>(std::unique_ptr<Operator> &&a,
30 : const FiniteElementSpace &fespace)
31 : {
32 0 : return std::make_unique<ParOperator>(std::move(a), fespace);
33 : }
34 :
35 : template <>
36 : auto BuildLevelParOperator<ComplexOperator>(std::unique_ptr<Operator> &&a,
37 : const FiniteElementSpace &fespace)
38 : {
39 0 : return std::make_unique<ComplexParOperator>(std::move(a), nullptr, fespace);
40 : }
41 :
42 : } // namespace
43 :
44 : template <typename VecType>
45 0 : DivFreeSolver<VecType>::DivFreeSolver(
46 : const MaterialOperator &mat_op, FiniteElementSpace &nd_fespace,
47 : FiniteElementSpaceHierarchy &h1_fespaces,
48 : const std::vector<mfem::Array<int>> &h1_bdr_tdof_lists, double tol, int max_it,
49 0 : int print)
50 : {
51 0 : BlockTimer bt(Timer::DIV_FREE);
52 :
53 : // If no boundaries on the mesh have been marked, add a single degree of freedom
54 : // constraint so the system for the projection is not singular. This amounts to enforcing
55 : // a scalar potential of 0 at a point in space if it is otherwise completely
56 : // unconstrained.
57 : const auto *ptr_h1_bdr_tdof_lists = &h1_bdr_tdof_lists;
58 : {
59 0 : MFEM_VERIFY(
60 : !h1_bdr_tdof_lists.empty(),
61 : "Unexpected empty list of boundary true dofs for finite element space hierarchy!");
62 0 : HYPRE_BigInt coarse_bdr_tdofs = h1_bdr_tdof_lists[0].Size();
63 : MPI_Comm comm = h1_fespaces.GetFESpaceAtLevel(0).GetComm();
64 : Mpi::GlobalSum(1, &coarse_bdr_tdofs, comm);
65 0 : if (coarse_bdr_tdofs == 0)
66 : {
67 0 : int root = (h1_fespaces.GetFESpaceAtLevel(0).GetTrueVSize() == 0) ? Mpi::Size(comm)
68 : : Mpi::Rank(comm);
69 : Mpi::GlobalMin(1, &root, comm);
70 0 : MFEM_VERIFY(root < Mpi::Size(comm),
71 : "No root process found for single true dof constraint!");
72 0 : if (root == Mpi::Rank(comm))
73 : {
74 0 : aux_tdof_lists.reserve(h1_fespaces.GetNumLevels());
75 0 : for (std::size_t l = 0; l < h1_fespaces.GetNumLevels(); l++)
76 : {
77 0 : auto &tdof_list = aux_tdof_lists.emplace_back(1);
78 0 : tdof_list[0] = 0;
79 : }
80 : ptr_h1_bdr_tdof_lists = &aux_tdof_lists;
81 : }
82 : }
83 : }
84 :
85 : // Create the mass and weak divergence operators for divergence-free projection.
86 0 : MaterialPropertyCoefficient epsilon_func(mat_op.GetAttributeToMaterial(),
87 : mat_op.GetPermittivityReal());
88 : {
89 : constexpr bool skip_zeros = false;
90 : BilinearForm m(h1_fespaces.GetFinestFESpace());
91 0 : m.AddDomainIntegrator<DiffusionIntegrator>(epsilon_func);
92 : // m.AssembleQuadratureData();
93 0 : auto m_vec = m.Assemble(h1_fespaces, skip_zeros);
94 0 : auto M_mg =
95 0 : std::make_unique<BaseMultigridOperator<OperType>>(h1_fespaces.GetNumLevels());
96 0 : for (std::size_t l = 0; l < h1_fespaces.GetNumLevels(); l++)
97 : {
98 : const auto &h1_fespace_l = h1_fespaces.GetFESpaceAtLevel(l);
99 : auto M_l = BuildLevelParOperator<OperType>(std::move(m_vec[l]), h1_fespace_l);
100 0 : M_l->SetEssentialTrueDofs((*ptr_h1_bdr_tdof_lists)[l],
101 : Operator::DiagonalPolicy::DIAG_ONE);
102 0 : if (l == h1_fespaces.GetNumLevels() - 1)
103 : {
104 0 : bdr_tdof_list_M = M_l->GetEssentialTrueDofs();
105 : }
106 0 : M_mg->AddOperator(std::move(M_l));
107 : }
108 : M = std::move(M_mg);
109 0 : }
110 : {
111 : // Weak divergence operator is always partially assembled.
112 : BilinearForm weakdiv(nd_fespace, h1_fespaces.GetFinestFESpace());
113 0 : weakdiv.AddDomainIntegrator<MixedVectorWeakDivergenceIntegrator>(epsilon_func);
114 0 : WeakDiv = std::make_unique<ParOperator>(weakdiv.PartialAssemble(), nd_fespace,
115 0 : h1_fespaces.GetFinestFESpace(), false);
116 : }
117 0 : Grad = &nd_fespace.GetDiscreteInterpolator(h1_fespaces.GetFinestFESpace());
118 :
119 : // The system matrix for the projection is real and SPD.
120 0 : auto amg = std::make_unique<MfemWrapperSolver<OperType>>(
121 0 : std::make_unique<BoomerAmgSolver>(1, 1, true, 0));
122 0 : std::unique_ptr<Solver<OperType>> pc;
123 0 : if (h1_fespaces.GetNumLevels() > 1)
124 : {
125 0 : const int mg_smooth_order =
126 0 : std::max(h1_fespaces.GetFinestFESpace().GetMaxElementOrder(), 2);
127 0 : pc = std::make_unique<GeometricMultigridSolver<OperType>>(
128 0 : h1_fespaces.GetFinestFESpace().GetComm(), std::move(amg),
129 0 : h1_fespaces.GetProlongationOperators(), nullptr, 1, 1, mg_smooth_order, 1.0, 0.0,
130 0 : true);
131 : }
132 : else
133 : {
134 : pc = std::move(amg);
135 : }
136 :
137 0 : auto pcg =
138 0 : std::make_unique<CgSolver<OperType>>(h1_fespaces.GetFinestFESpace().GetComm(), print);
139 0 : pcg->SetInitialGuess(false);
140 : pcg->SetRelTol(tol);
141 : pcg->SetAbsTol(std::numeric_limits<double>::epsilon());
142 : pcg->SetMaxIter(max_it);
143 :
144 0 : ksp = std::make_unique<BaseKspSolver<OperType>>(std::move(pcg), std::move(pc));
145 0 : ksp->SetOperators(*M, *M);
146 :
147 0 : psi.SetSize(h1_fespaces.GetFinestFESpace().GetTrueVSize());
148 0 : rhs.SetSize(h1_fespaces.GetFinestFESpace().GetTrueVSize());
149 0 : psi.UseDevice(true);
150 0 : rhs.UseDevice(true);
151 0 : }
152 :
153 : template <typename VecType>
154 0 : void DivFreeSolver<VecType>::Mult(VecType &y) const
155 : {
156 0 : BlockTimer bt(Timer::DIV_FREE);
157 :
158 : // Compute the divergence of y.
159 : if constexpr (std::is_same<VecType, ComplexVector>::value)
160 : {
161 0 : WeakDiv->Mult(y.Real(), rhs.Real());
162 0 : WeakDiv->Mult(y.Imag(), rhs.Imag());
163 : }
164 : else
165 : {
166 0 : WeakDiv->Mult(y, rhs);
167 : }
168 :
169 : // Apply essential BC and solve the linear system.
170 0 : if (bdr_tdof_list_M)
171 : {
172 0 : linalg::SetSubVector(rhs, *bdr_tdof_list_M, 0.0);
173 : }
174 0 : ksp->Mult(rhs, psi);
175 :
176 : // Compute the irrotational portion of y and subtract.
177 : if constexpr (std::is_same<VecType, ComplexVector>::value)
178 : {
179 0 : Grad->AddMult(psi.Real(), y.Real(), 1.0);
180 0 : Grad->AddMult(psi.Imag(), y.Imag(), 1.0);
181 : }
182 : else
183 : {
184 0 : Grad->AddMult(psi, y, 1.0);
185 : }
186 0 : }
187 :
188 : template class DivFreeSolver<Vector>;
189 : template class DivFreeSolver<ComplexVector>;
190 :
191 : } // namespace palace
|
