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 "ams.hpp"
5 :
6 : #include "fem/bilinearform.hpp"
7 : #include "fem/fespace.hpp"
8 : #include "fem/integrator.hpp"
9 : #include "linalg/hypre.hpp"
10 : #include "linalg/rap.hpp"
11 : #include "utils/omp.hpp"
12 :
13 : namespace palace
14 : {
15 :
16 0 : HypreAmsSolver::HypreAmsSolver(FiniteElementSpace &nd_fespace,
17 : FiniteElementSpace &h1_fespace, int cycle_it, int smooth_it,
18 : bool vector_interp, bool singular_op, bool agg_coarsen,
19 0 : int print)
20 : : mfem::HypreSolver(),
21 : // From the Hypre docs for AMS: cycles 1, 5, 8, 11, 13 are fastest, 7 yields fewest its
22 : // (MFEM default is 13). 14 is similar to 11/13 but is cheaper in that is uses additive
23 : // scalar Pi-space corrections.
24 0 : cycle_type(vector_interp ? 1 : 14), space_dim(nd_fespace.SpaceDimension()),
25 : // When used as the coarse solver of geometric multigrid, always do only a single
26 : // V-cycle.
27 0 : ams_it(cycle_it), ams_smooth_it(smooth_it),
28 : // If we know the operator is singular (no mass matrix, for magnetostatic problems),
29 : // internally the AMS solver will avoid G-space corrections.
30 0 : ams_singular(singular_op),
31 : // For positive (SPD) operators, we will use aggressive coarsening but not for frequency
32 : // domain problems when the preconditioner matrix is not SPD.
33 0 : agg_coarsen(agg_coarsen), print((print > 1) ? print - 1 : 0)
34 : {
35 : // From MFEM: The AMS preconditioner may sometimes require inverting singular matrices
36 : // with BoomerAMG, which are handled correctly in Hypre's Solve method, but can produce
37 : // Hypre errors in the Setup (specifically in the row l1-norm computation). See the
38 : // documentation of MFEM's SetErrorMode() for more details.
39 0 : error_mode = IGNORE_HYPRE_ERRORS;
40 :
41 : // Set up the AMS solver.
42 0 : ConstructAuxiliaryMatrices(nd_fespace, h1_fespace);
43 0 : InitializeSolver();
44 0 : }
45 :
46 0 : HypreAmsSolver::~HypreAmsSolver()
47 : {
48 0 : HYPRE_AMSDestroy(ams);
49 0 : }
50 :
51 0 : void HypreAmsSolver::ConstructAuxiliaryMatrices(FiniteElementSpace &nd_fespace,
52 : FiniteElementSpace &h1_fespace)
53 : {
54 : // Set up the auxiliary space objects for the preconditioner. Mostly the same as MFEM's
55 : // HypreAMS:Init. Start with the discrete gradient matrix. We don't skip zeros for the
56 : // full assembly to accelerate things on GPU and since they shouldn't affect the sparsity
57 : // pattern of the parallel G^T A G matrix (computed by Hypre).
58 0 : const bool skip_zeros_interp = !mfem::Device::Allows(mfem::Backend::DEVICE_MASK);
59 : {
60 : const auto *PtGP =
61 0 : dynamic_cast<const ParOperator *>(&nd_fespace.GetDiscreteInterpolator(h1_fespace));
62 0 : MFEM_VERIFY(
63 : PtGP,
64 : "HypreAmsSolver requires the discrete gradient matrix as a ParOperator operator!");
65 0 : G = &PtGP->ParallelAssemble(skip_zeros_interp);
66 : }
67 :
68 : // Vertex coordinates for the lowest order case, or Nedelec interpolation matrix or
69 : // matrices for order > 1. Expects that Mesh::SetVerticesFromNodes has been called at some
70 : // point to avoid calling GridFunction::GetNodalValues here.
71 : mfem::ParMesh &mesh = h1_fespace.GetParMesh();
72 0 : if (h1_fespace.GetMaxElementOrder() == 1)
73 : {
74 0 : mfem::ParGridFunction x_coord(&h1_fespace.Get()), y_coord(&h1_fespace.Get()),
75 0 : z_coord(&h1_fespace.Get());
76 0 : MFEM_VERIFY(x_coord.Size() == mesh.GetNV(),
77 : "Unexpected size for vertex coordinates in AMS setup!");
78 0 : PalacePragmaOmp(parallel for schedule(static))
79 : for (int i = 0; i < mesh.GetNV(); i++)
80 : {
81 : x_coord(i) = mesh.GetVertex(i)[0];
82 : if (space_dim > 1)
83 : {
84 : y_coord(i) = mesh.GetVertex(i)[1];
85 : }
86 : if (space_dim > 2)
87 : {
88 : z_coord(i) = mesh.GetVertex(i)[2];
89 : }
90 : }
91 0 : x.reset(x_coord.ParallelProject());
92 0 : x->HypreReadWrite();
93 0 : if (space_dim > 1)
94 : {
95 0 : y.reset(y_coord.ParallelProject());
96 0 : y->HypreReadWrite();
97 : }
98 0 : if (space_dim > 2)
99 : {
100 0 : z.reset(z_coord.ParallelProject());
101 0 : z->HypreReadWrite();
102 : }
103 0 : }
104 : else
105 : {
106 : // Fall back to MFEM legacy assembly for identity interpolator.
107 0 : FiniteElementSpace h1d_fespace(h1_fespace.GetMesh(), &h1_fespace.GetFEColl(), space_dim,
108 0 : mfem::Ordering::byVDIM);
109 : mfem::DiscreteLinearOperator pi(&h1d_fespace.Get(), &nd_fespace.Get());
110 0 : pi.AddDomainInterpolator(new mfem::IdentityInterpolator);
111 0 : pi.SetAssemblyLevel(mfem::AssemblyLevel::LEGACY);
112 0 : pi.Assemble(skip_zeros_interp);
113 0 : pi.Finalize(skip_zeros_interp);
114 0 : ParOperator RAP_Pi(std::make_unique<hypre::HypreCSRMatrix>(pi.SpMat()), h1d_fespace,
115 0 : nd_fespace, true);
116 : Pi = RAP_Pi.StealParallelAssemble();
117 0 : if (cycle_type >= 10)
118 : {
119 : // Get blocks of Pi corresponding to each component, and free Pi.
120 0 : mfem::Array2D<mfem::HypreParMatrix *> Pi_blocks(1, space_dim);
121 0 : Pi->GetBlocks(Pi_blocks, false, true);
122 0 : Pix.reset(Pi_blocks(0, 0));
123 0 : if (space_dim > 1)
124 : {
125 0 : Piy.reset(Pi_blocks(0, 1));
126 : }
127 0 : if (space_dim > 2)
128 : {
129 0 : Piz.reset(Pi_blocks(0, 2));
130 : }
131 : Pi.reset();
132 : }
133 0 : }
134 0 : }
135 :
136 0 : void HypreAmsSolver::InitializeSolver()
137 : {
138 : // Create the Hypre solver object.
139 0 : HYPRE_AMSCreate(&ams);
140 0 : HYPRE_AMSSetDimension(ams, space_dim);
141 0 : HYPRE_AMSSetCycleType(ams, cycle_type);
142 :
143 : // Control printing and number of iterations for use as a preconditioner.
144 0 : HYPRE_AMSSetPrintLevel(ams, print);
145 0 : HYPRE_AMSSetMaxIter(ams, ams_it);
146 : // HYPRE_AMSSetTol(ams, 1.0e-16); // Avoid issues with zero RHS
147 :
148 : // Set this option when solving a curl-curl problem with zero mass term.
149 0 : if (ams_singular)
150 : {
151 0 : HYPRE_AMSSetBetaPoissonMatrix(ams, nullptr);
152 : }
153 :
154 : // Set additional AMS options.
155 : int coarsen_type = 10; // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
156 0 : int amg_agg_levels = agg_coarsen ? 1 : 0; // Number of aggressive coarsening levels
157 : double theta = 0.5; // AMG strength parameter = 0.25 is 2D optimal (0.5-0.8 for 3D)
158 : int amg_relax_type = 8; // 3 = GS, 6 = symm. GS, 8 = l1-symm. GS, 13 = l1-GS,
159 : // 18 = l1-Jacobi, 16 = Chebyshev
160 : int interp_type = 6; // 6 = Extended+i, 0 = Classical, 13 = FF1
161 : int Pmax = 4; // Interpolation width
162 : int relax_type = 2; // 2 = l1-SSOR, 4 = trunc. l1-SSOR, 1 = l1-Jacobi, 16 = Chebyshev
163 : double weight = 1.0;
164 : double omega = 1.0;
165 0 : if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK))
166 : {
167 : // Modify options for GPU-supported features.
168 : coarsen_type = 8;
169 : amg_agg_levels = 0;
170 : amg_relax_type = 18;
171 : relax_type = 1;
172 : }
173 :
174 0 : HYPRE_AMSSetSmoothingOptions(ams, relax_type, ams_smooth_it, weight, omega);
175 0 : HYPRE_AMSSetAlphaAMGOptions(ams, coarsen_type, amg_agg_levels, amg_relax_type, theta,
176 : interp_type, Pmax);
177 0 : HYPRE_AMSSetBetaAMGOptions(ams, coarsen_type, amg_agg_levels, amg_relax_type, theta,
178 : interp_type, Pmax);
179 :
180 : // int coarse_relax_type = 8; // Default, l1-symm. GS
181 : // HYPRE_AMSSetAlphaAMGCoarseRelaxType(ams, coarse_relax_type);
182 : // HYPRE_AMSSetBetaAMGCoarseRelaxType(ams, coarse_relax_type);
183 :
184 : // Set the discrete gradient matrix.
185 0 : HYPRE_AMSSetDiscreteGradient(ams, (HYPRE_ParCSRMatrix)*G);
186 :
187 : // Set the mesh vertex coordinates or Nedelec interpolation matrix or matrices.
188 0 : HYPRE_ParVector HY_X = (x) ? (HYPRE_ParVector)*x : nullptr;
189 0 : HYPRE_ParVector HY_Y = (y) ? (HYPRE_ParVector)*y : nullptr;
190 0 : HYPRE_ParVector HY_Z = (z) ? (HYPRE_ParVector)*z : nullptr;
191 0 : HYPRE_AMSSetCoordinateVectors(ams, HY_X, HY_Y, HY_Z);
192 :
193 : HYPRE_ParCSRMatrix HY_Pi = (Pi) ? (HYPRE_ParCSRMatrix)*Pi : nullptr;
194 0 : HYPRE_ParCSRMatrix HY_Pix = (Pix) ? (HYPRE_ParCSRMatrix)*Pix : nullptr;
195 0 : HYPRE_ParCSRMatrix HY_Piy = (Piy) ? (HYPRE_ParCSRMatrix)*Piy : nullptr;
196 0 : HYPRE_ParCSRMatrix HY_Piz = (Piz) ? (HYPRE_ParCSRMatrix)*Piz : nullptr;
197 0 : HYPRE_AMSSetInterpolations(ams, HY_Pi, HY_Pix, HY_Piy, HY_Piz);
198 0 : }
199 :
200 0 : void HypreAmsSolver::SetOperator(const Operator &op)
201 : {
202 : // When the operator changes, we need to rebuild the AMS solver but can use the unchanged
203 : // auxiliary space matrices.
204 0 : if (A)
205 : {
206 0 : HYPRE_AMSDestroy(ams);
207 0 : InitializeSolver();
208 : }
209 0 : A = const_cast<mfem::HypreParMatrix *>(dynamic_cast<const mfem::HypreParMatrix *>(&op));
210 0 : MFEM_VERIFY(A, "HypreAmsSolver requires a HypreParMatrix operator!");
211 0 : height = A->Height();
212 0 : width = A->Width();
213 :
214 : // From mfem::HypreAMS: Update HypreSolver base class.
215 0 : setup_called = 0;
216 0 : delete X;
217 0 : delete B;
218 0 : B = X = nullptr;
219 0 : auxB.Delete();
220 : auxB.Reset();
221 0 : auxX.Delete();
222 : auxX.Reset();
223 0 : }
224 :
225 : } // namespace palace
|
