MFEM  v4.5.2
Finite element discretization library
lor_ams.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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"
13 #include "../../general/forall.hpp"
14 #include "../../fem/pbilinearform.hpp"
15 
16 namespace 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(i, nedge_dof+1, 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*>(
167  vert_fes.GetElementRestriction(ordering));
168  const auto *R_e = dynamic_cast<const ElementRestriction*>(
169  edge_fes.GetElementRestriction(ordering));
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
183  G_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
185 
186  auto J = G_local.WriteJ();
187  auto V = G_local.WriteData();
188 
189  // Loop over Nedelec L-DOFs
190  MFEM_FORALL(i, nedge_dof,
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 
242 template <typename T>
243 static inline const T *HypreRead(const Memory<T> &mem)
244 {
245  return mem.Read(GetHypreMemoryClass(), mem.Capacity());
246 }
247 
248 template <typename T>
249 static 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 (dim,
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, dim)
264  // and then make the hypre vectors x, y, and z point to subvectors.
265  //
266  // In 2D, 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 = dim;
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*dim);
285 
286  auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, dim);
287  const auto xyz_e =
288  Reshape(HypreRead(X_vert.GetMemory()), dim, 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(i, ntdofs,
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
305  HYPRE_BigInt glob_size = vert_fes.GlobalTrueVSize();
307 
308  double *d_x_ptr = xyz_tv + 0*ntdofs;
309  x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
310  double *d_y_ptr = xyz_tv + 1*ntdofs;
311  y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
312  if (dim == 3)
313  {
314  double *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 
324 {
325  return StealPointer(G);
326 }
327 
329 {
330  return StealPointer(xyz_tvec);
331 }
332 
334 {
335  return StealPointer(x);
336 }
337 
339 {
340  return StealPointer(y);
341 }
342 
344 {
345  return StealPointer(z);
346 }
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 
403 void LORSolver<HypreAMS>::Mult(const Vector &x, Vector &y) const
404 {
405  solver->Mult(x, y);
406 }
407 
409 
410 const HypreAMS &LORSolver<HypreAMS>::GetSolver() const { return *solver; }
411 
413 {
414  delete solver;
415  delete xyz;
416 }
417 
418 #endif
419 
420 } // namespace mfem
HypreParVector * z
Definition: lor_ams.hpp:41
const T * Read(MemoryClass mc, int size) const
Get read-only access to the memory with the given MemoryClass.
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
Create and assemble a low-order refined version of a ParBilinearForm.
Definition: lor.hpp:165
HypreParVector * StealZCoordinate()
Definition: lor_ams.cpp:343
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
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
LORBase * lor
Definition: lor.hpp:206
ParFiniteElementSpace & edge_fes
The Nedelec space.
Definition: lor_ams.hpp:31
Device memory; using CUDA or HIP *Malloc and *Free.
HypreParVector * x
Definition: lor_ams.hpp:41
Efficient batched assembly of LOR discretizations on device.
Definition: lor_batched.hpp:33
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void Form2DEdgeToVertex_ND(Array< int > &edge2vert)
Definition: lor_ams.cpp:38
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: lor.hpp:255
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void FormGradientMatrix()
Construct the discrete gradient matrix (not part of the public API).
Definition: lor_ams.cpp:139
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:389
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:269
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:230
static bool FormIsSupported(BilinearForm &a)
Returns true if the form a supports batched assembly, false otherwise.
Definition: lor_batched.cpp:49
Memory< int > & GetMemoryI()
Definition: sparsemat.hpp:237
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1259
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:253
const OperatorHandle & GetAssembledSystem() const
Returns the assembled LOR system (const version).
Definition: lor.cpp:270
ParFiniteElementSpace vert_fes
The corresponding H1 space.
Definition: lor_ams.hpp:35
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
Make A_hyp steal ownership of its diagonal part A_diag.
Definition: hypre.cpp:2741
Vector * StealCoordinateVector()
Definition: lor_ams.cpp:328
Vector * xyz_tvec
Mesh vertex coordinates in true-vector format.
Definition: lor_ams.hpp:36
Data type sparse matrix.
Definition: sparsemat.hpp:50
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
void Form2DEdgeToVertex(Array< int > &edge2vert)
Definition: lor_ams.cpp:21
SolverType solver
Definition: lor.hpp:208
HypreParMatrix * G
Discrete gradient matrix.
Definition: lor_ams.hpp:37
const int dim
Spatial dimension.
Definition: lor_ams.hpp:32
void Form3DEdgeToVertex(Array< int > &edge2vert)
Definition: lor_ams.cpp:105
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:241
void FormCoordinateVectors(const Vector &X_vert)
Construct the mesh coordinate vectors (not part of the public API).
Definition: lor_ams.cpp:254
int Capacity() const
Return the size of the allocated memory.
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: lor.hpp:248
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
HypreParVector * StealYCoordinate()
Definition: lor_ams.cpp:338
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
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
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
const Vector & GetLORVertexCoordinates()
Return the vertices of the LOR mesh in E-vector format.
Definition: lor_batched.hpp:74
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:720
SolverType & GetSolver()
Access the underlying solver.
Definition: lor.hpp:258
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:37
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:283
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
void Assemble(BilinearForm &a, const Array< int > ess_dofs, OperatorHandle &A)
Assemble the given form as a matrix and place the result in A.
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
const int order
Polynomial degree.
Definition: lor_ams.hpp:33
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
HypreParVector * y
Definition: lor_ams.hpp:41
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:257
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:281
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition: array.hpp:319
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
constexpr MemoryClass GetHypreMemoryClass()
The MemoryClass used by Hypre objects.
Definition: hypre.hpp:137
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
HypreParVector * StealXCoordinate()
Definition: lor_ams.cpp:333
int dim
Definition: ex24.cpp:53
void Form2DEdgeToVertex_RT(Array< int > &edge2vert)
Definition: lor_ams.cpp:72
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Definition: sparsemat.cpp:297
Class used by MFEM to store pointers to host and/or device memory.
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
Lexicographic ordering for tensor-product FiniteElements.
Class for parallel bilinear form.
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
Vector data type.
Definition: vector.hpp:60
ID for class HypreParMatrix.
Definition: operator.hpp:287
T * StealPointer(T *&ptr)
Definition: lor_ams.hpp:95
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
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
HypreParMatrix * StealGradientMatrix()
Definition: lor_ams.cpp:323
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:282
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:73
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
double * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:273