MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
lor_ams.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "lor_ams.hpp"
15
16namespace mfem
17{
18
19#ifdef MFEM_USE_MPI
20
22{
24 if (dynamic_cast<const ND_FECollection*>(fec))
25 {
26 Form2DEdgeToVertex_ND(edge2vert);
27 }
28 else if (dynamic_cast<const RT_FECollection*>(fec))
29 {
30 Form2DEdgeToVertex_RT(edge2vert);
31 }
32 else
33 {
34 MFEM_ABORT("Bad finite element type.")
35 }
36}
37
39{
40 const int o = order;
41 const int op1 = o + 1;
42 const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
43
44 edge2vert.SetSize(2*nedge);
45 auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
46
47 for (int c=0; c<dim; ++c)
48 {
49 const int nx = (c == 0) ? o : op1;
50 for (int i2=0; i2<op1; ++i2)
51 {
52 for (int i1=0; i1<o; ++i1)
53 {
54 const int ix = (c == 0) ? i1 : i2;
55 const int iy = (c == 0) ? i2 : i1;
56
57 const int iedge = ix + iy*nx + c*o*op1;
58
59 const int ix1 = (c == 0) ? ix + 1 : ix;
60 const int iy1 = (c == 1) ? iy + 1 : iy;
61
62 const int iv0 = ix + iy*op1;
63 const int iv1 = ix1 + iy1*op1;
64
65 e2v(0, iedge) = iv0;
66 e2v(1, iedge) = iv1;
67 }
68 }
69 }
70}
71
73{
74 const int o = order;
75 const int op1 = o + 1;
76 const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
77
78 edge2vert.SetSize(2*nedge);
79 auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
80
81 for (int c=0; c<dim; ++c)
82 {
83 const int nx = (c == 0) ? op1 : o;
84 for (int i=0; i<o*op1; ++i)
85 {
86 const int ix = i%nx;
87 const int iy = i/nx;
88
89 const int iedge = ix + iy*nx + c*o*op1;
90
91 const int ix1 = (c == 0) ? ix : ix + 1;
92 const int iy1 = (c == 1) ? iy : iy + 1;
93
94 const int iv0 = ix + iy*op1;
95 const int iv1 = ix1 + iy1*op1;
96
97 // Rotated gradient in 2D (-dy, dx), so flip the sign for the first
98 // component (c == 0).
99 e2v(0, iedge) = (c == 1) ? iv0 : iv1;
100 e2v(1, iedge) = (c == 1) ? iv1 : iv0;
101 }
102 }
103}
104
106{
107 const int o = order;
108 const int op1 = o + 1;
109 const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
110
111 edge2vert.SetSize(2*nedge);
112 auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
113
114 for (int c=0; c<dim; ++c)
115 {
116 const int nx = (c == 0) ? o : op1;
117 const int ny = (c == 1) ? o : op1;
118 for (int i=0; i<o*op1*op1; ++i)
119 {
120 const int ix = i%nx;
121 const int iy = (i/nx)%ny;
122 const int iz = i/nx/ny;
123
124 const int iedge = ix + iy*nx + iz*nx*ny + c*o*op1*op1;
125
126 const int ix1 = (c == 0) ? ix + 1 : ix;
127 const int iy1 = (c == 1) ? iy + 1 : iy;
128 const int iz1 = (c == 2) ? iz + 1 : iz;
129
130 const int iv0 = ix + iy*op1 + iz*op1*op1;
131 const int iv1 = ix1 + iy1*op1 + iz1*op1*op1;
132
133 e2v(0, iedge) = iv0;
134 e2v(1, iedge) = iv1;
135 }
136 }
137}
138
140{
141 // The gradient matrix maps from LOR vertices to LOR edges. Given an edge
142 // (defined by its two vertices) e_i = (v_j1, v_j2), the matrix has nonzeros
143 // A(i, j1) = -1 and A(i, j2) = 1, so there are always exactly two nonzeros
144 // per row.
145 const int nedge_dof = edge_fes.GetNDofs();
146 const int nvert_dof = vert_fes.GetNDofs();
147
148 SparseMatrix G_local;
149 G_local.OverrideSize(nedge_dof, nvert_dof);
150
151 G_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
152 // Each row always has two nonzeros
153 const int nnz = 2*nedge_dof;
154 auto I = G_local.WriteI();
155 mfem::forall(nedge_dof+1, [=] MFEM_HOST_DEVICE (int i) { I[i] = 2*i; });
156
157 // edge2vertex is a mapping of size (2, nedge_per_el), such that with a macro
158 // element, edge i (in lexicographic ordering) has vertices (also in
159 // lexicographic ordering) given by the entries (0, i) and (1, i) of the
160 // matrix.
161 Array<int> edge2vertex;
162 if (dim == 2) { Form2DEdgeToVertex(edge2vertex); }
163 else { Form3DEdgeToVertex(edge2vertex); }
164
166 const auto *R_v = dynamic_cast<const ElementRestriction*>(
168 const auto *R_e = dynamic_cast<const ElementRestriction*>(
170 MFEM_VERIFY(R_v != NULL && R_e != NULL, "");
171
172 const int nel_ho = edge_fes.GetNE();
173 const int nedge_per_el = static_cast<int>(dim*order*pow(order + 1, dim - 1));
174 const int nvert_per_el = static_cast<int>(pow(order + 1, dim));
175
176 const auto offsets_e = R_e->Offsets().Read();
177 const auto indices_e = R_e->Indices().Read();
178 const auto gather_v = Reshape(R_v->GatherMap().Read(), nvert_per_el, nel_ho);
179
180 const auto e2v = Reshape(edge2vertex.Read(), 2, nedge_per_el);
181
182 // Fill J and data
185
186 auto J = G_local.WriteJ();
187 auto V = G_local.WriteData();
188
189 // Loop over Nedelec L-DOFs
190 mfem::forall(nedge_dof, [=] MFEM_HOST_DEVICE (int i)
191 {
192 const int sj = indices_e[offsets_e[i]]; // signed
193 const int j = (sj >= 0) ? sj : -1 - sj;
194 const int sgn = (sj >= 0) ? 1 : -1;
195 const int j_loc = j%nedge_per_el;
196 const int j_el = j/nedge_per_el;
197
198 const int jv0_loc = e2v(0, j_loc);
199 const int jv1_loc = e2v(1, j_loc);
200
201 J[i*2 + 0] = gather_v(jv0_loc, j_el);
202 J[i*2 + 1] = gather_v(jv1_loc, j_el);
203
204 V[i*2 + 0] = -sgn;
205 V[i*2 + 1] = sgn;
206 });
207
208 // Create a block diagonal parallel matrix
210 G_diag.MakeRectangularBlockDiag(vert_fes.GetComm(),
215 &G_local);
216
217 // Assemble the parallel gradient matrix, must be deleted by the caller
219 {
220 G = G_diag.As<HypreParMatrix>();
221 G_diag.SetOperatorOwner(false);
222 HypreStealOwnership(*G, G_local);
223 }
224 else
225 {
233 Rt.As<SparseMatrix>());
234 G = RAP(Rt_diag.As<HypreParMatrix>(),
235 G_diag.As<HypreParMatrix>(),
237 }
238 G->CopyRowStarts();
239 G->CopyColStarts();
240}
241
242template <typename T>
243static inline const T *HypreRead(const Memory<T> &mem)
244{
245 return mem.Read(GetHypreMemoryClass(), mem.Capacity());
246}
247
248template <typename T>
249static inline T *HypreWrite(Memory<T> &mem)
250{
251 return mem.Write(GetHypreMemoryClass(), mem.Capacity());
252}
253
255{
256 // Create true-DOF vectors x, y, and z that contain the coordinates of the
257 // vertices of the LOR mesh. The vertex coordinates are already computed in
258 // E-vector format and passed in in X_vert.
259 //
260 // In this function, we need to convert X_vert (which has the shape (sdim,
261 // ndof_per_el, nel_ho)) to T-DOF format.
262 //
263 // We place the results in the vector xyz_tvec, which has shape (ntdofs, sdim)
264 // and then make the hypre vectors x, y, and z point to subvectors.
265 //
266 // When the space dimension is 2, z is NULL.
267
268 // Create the H1 vertex space and get the element restriction
270 const Operator *op = vert_fes.GetElementRestriction(ordering);
271 const auto *el_restr = dynamic_cast<const ElementRestriction*>(op);
272 MFEM_VERIFY(el_restr != NULL, "");
274
275 const int nel_ho = vert_fes.GetNE();
276 const int ndp1 = order + 1;
277 const int ndof_per_el = static_cast<int>(pow(ndp1, dim));
278 const int sdim = vert_fes.GetMesh()->SpaceDimension();
279 const int ntdofs = R->Height();
280
281 const MemoryClass mc = GetHypreMemoryClass();
282 bool dev = (mc == MemoryClass::DEVICE);
283
284 xyz_tvec = new Vector(ntdofs*sdim);
285
286 auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, sdim);
287 const auto xyz_e =
288 Reshape(HypreRead(X_vert.GetMemory()), sdim, ndof_per_el, nel_ho);
289 const auto d_offsets = HypreRead(el_restr->Offsets().GetMemory());
290 const auto d_indices = HypreRead(el_restr->Indices().GetMemory());
291 const auto ltdof_ldof = HypreRead(R->GetMemoryJ());
292
293 // Go from E-vector format directly to T-vector format
294 mfem::hypre_forall(ntdofs, [=] MFEM_HOST_DEVICE (int i)
295 {
296 const int j = d_offsets[ltdof_ldof[i]];
297 for (int c = 0; c < sdim; ++c)
298 {
299 const int idx_j = d_indices[j];
300 xyz_tv(i,c) = xyz_e(c, idx_j%ndof_per_el, idx_j/ndof_per_el);
301 }
302 });
303
304 // Make x, y, z HypreParVectors point to T-vector data
307
308 real_t *d_x_ptr = xyz_tv + 0*ntdofs;
309 x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
310 real_t *d_y_ptr = xyz_tv + 1*ntdofs;
311 y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
312 if (sdim == 3)
313 {
314 real_t *d_z_ptr = xyz_tv + 2*ntdofs;
315 z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);
316 }
317 else
318 {
319 z = NULL;
320 }
321}
322
327
332
337
342
347
349{
350 delete x;
351 delete y;
352 delete z;
353 delete xyz_tvec;
354 delete G;
355}
356
358 const Vector &X_vert)
359 : edge_fes(pfes_ho_),
360 dim(edge_fes.GetParMesh()->Dimension()),
361 order(edge_fes.GetMaxElementOrder()),
362 vert_fec(order, dim),
363 vert_fes(edge_fes.GetParMesh(), &vert_fec)
364{
365 FormCoordinateVectors(X_vert);
367}
368
370 ParBilinearForm &a_ho, const Array<int> &ess_tdof_list, int ref_type)
371{
373 {
374 ParFiniteElementSpace &pfes = *a_ho.ParFESpace();
375 BatchedLORAssembly batched_lor(pfes);
376 batched_lor.Assemble(a_ho, ess_tdof_list, A);
377 BatchedLOR_AMS lor_ams(pfes, batched_lor.GetLORVertexCoordinates());
378 xyz = lor_ams.StealCoordinateVector();
379 solver = new HypreAMS(*A.As<HypreParMatrix>(),
380 lor_ams.StealGradientMatrix(),
381 lor_ams.StealXCoordinate(),
382 lor_ams.StealYCoordinate(),
383 lor_ams.StealZCoordinate());
384 }
385 else
386 {
387 ParLORDiscretization lor(a_ho, ess_tdof_list, ref_type);
388 // Assume ownership of the system matrix so that `lor` can be safely
389 // deleted
390 A.Reset(lor.GetAssembledSystem().Ptr());
392 solver = new HypreAMS(lor.GetAssembledMatrix(), &lor.GetParFESpace());
393 }
394 width = solver->Width();
395 height = solver->Height();
396}
397
399{
400 solver->SetOperator(op);
401}
402
404{
405 solver->Mult(x, y);
406}
407
409
410const HypreAMS &LORSolver<HypreAMS>::GetSolver() const { return *solver; }
411
413{
414 delete solver;
415 delete xyz;
416}
417
418#endif
419
420} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition array.hpp:329
Efficient batched assembly of LOR discretizations on device.
const Vector & GetLORVertexCoordinates()
Return the vertices of the LOR mesh in E-vector format.
static bool FormIsSupported(BilinearForm &a)
Returns true if the form a supports batched assembly, false otherwise.
void Assemble(BilinearForm &a, const Array< int > ess_dofs, OperatorHandle &A)
Assemble the given form as a matrix and place the result in A.
HypreParVector * StealZCoordinate()
Definition lor_ams.cpp:343
HypreParVector * StealYCoordinate()
Definition lor_ams.cpp:338
HypreParMatrix * StealGradientMatrix()
Definition lor_ams.cpp:323
const int order
Polynomial degree.
Definition lor_ams.hpp:33
Vector * StealCoordinateVector()
Definition lor_ams.cpp:328
void FormGradientMatrix()
Construct the discrete gradient matrix (not part of the public API).
Definition lor_ams.cpp:139
HypreParVector * StealXCoordinate()
Definition lor_ams.cpp:333
HypreParVector * x
Definition lor_ams.hpp:41
HypreParMatrix * G
Discrete gradient matrix.
Definition lor_ams.hpp:37
void Form2DEdgeToVertex_RT(Array< int > &edge2vert)
Definition lor_ams.cpp:72
ParFiniteElementSpace & edge_fes
The Nedelec space.
Definition lor_ams.hpp:31
const int dim
Spatial dimension.
Definition lor_ams.hpp:32
BatchedLOR_AMS(ParFiniteElementSpace &pfes_ho_, const Vector &X_vert)
Construct the BatchedLOR_AMS object associated with the Nedelec space (or RT in 2D) pfes_ho_.
Definition lor_ams.cpp:357
Vector * xyz_tvec
Mesh vertex coordinates in true-vector format.
Definition lor_ams.hpp:36
ParFiniteElementSpace vert_fes
The corresponding H1 space.
Definition lor_ams.hpp:35
void FormCoordinateVectors(const Vector &X_vert)
Construct the mesh coordinate vectors (not part of the public API).
Definition lor_ams.cpp:254
HypreParVector * z
Definition lor_ams.hpp:41
void Form2DEdgeToVertex(Array< int > &edge2vert)
Definition lor_ams.cpp:21
void Form3DEdgeToVertex(Array< int > &edge2vert)
Definition lor_ams.cpp:105
void Form2DEdgeToVertex_ND(Array< int > &edge2vert)
Definition lor_ams.cpp:38
HypreParVector * y
Definition lor_ams.hpp:41
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:274
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:710
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1336
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
const OperatorHandle & GetAssembledSystem() const
Returns the assembled LOR system (const version).
Definition lor.cpp:270
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition lor.hpp:255
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition lor.hpp:248
SolverType & GetSolver()
Access the underlying solver.
Definition lor.hpp:258
LORSolver(BilinearForm &a_ho, const Array< int > &ess_tdof_list, int ref_type=BasisType::GaussLobatto)
Create a solver of type SolverType, formed using the assembled SparseMatrix of the LOR version of a_h...
Definition lor.hpp:212
Class used by MFEM to store pointers to host and/or device memory.
int Capacity() const
Return the size of the allocated memory.
const T * Read(MemoryClass mc, int size) const
Get read-only access to the memory with the given MemoryClass.
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition handle.hpp:120
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
void MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_num_rows, HYPRE_BigInt glob_num_cols, HYPRE_BigInt *row_starts, HYPRE_BigInt *col_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel rectangular block-diagonal matrix using the currently set...
Definition handle.cpp:91
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
@ Hypre_ParCSR
ID for class HypreParMatrix.
Definition operator.hpp:287
Class for parallel bilinear form.
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
HYPRE_BigInt * GetTrueDofOffsets() const
Definition pfespace.hpp:282
HYPRE_BigInt GlobalVSize() const
Definition pfespace.hpp:283
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
HYPRE_BigInt * GetDofOffsets() const
Definition pfespace.hpp:281
const Operator * GetProlongationMatrix() const override
The returned Operator is owned by the FiniteElementSpace.
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition pfespace.hpp:327
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
Definition pfespace.hpp:389
Create and assemble a low-order refined version of a ParBilinearForm.
Definition lor.hpp:166
ParFiniteElementSpace & GetParFESpace() const
Return the LOR ParFiniteElementSpace.
Definition lor.cpp:528
HypreParMatrix & GetAssembledMatrix() const
Return the assembled LOR operator as a HypreParMatrix.
Definition lor.cpp:522
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Data type sparse matrix.
Definition sparsemat.hpp:51
int * WriteJ(bool on_dev=true)
Memory< int > & GetMemoryI()
int * WriteI(bool on_dev=true)
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Memory< int > & GetMemoryJ()
Memory< real_t > & GetMemoryData()
real_t * WriteData(bool on_dev=true)
Vector data type.
Definition vector.hpp:80
Memory< real_t > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition vector.hpp:249
int dim
Definition ex24.cpp:53
HYPRE_Int HYPRE_BigInt
T * StealPointer(T *&ptr)
Definition lor_ams.hpp:95
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
MemoryClass
Memory classes identify sets of memory types.
MemoryClass GetHypreMemoryClass()
The MemoryClass used by Hypre objects.
Definition hypre.hpp:154
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
bool IsIdentityProlongation(const Operator *P)
Definition operator.hpp:720
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
Make A_hyp steal ownership of its diagonal part A_diag.
Definition hypre.cpp:2840
void hypre_forall(int N, lambda &&body)
Definition forall.hpp:823
float real_t
Definition config.hpp:43
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition fespace.hpp:75
@ LEXICOGRAPHIC
Lexicographic ordering for tensor-product FiniteElements.
void forall(int N, lambda &&body)
Definition forall.hpp:754