MFEM  v4.6.0
Finite element discretization library
lor_ads.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_ads.hpp"
13 #include "../../general/forall.hpp"
14 #include "../../fem/pbilinearform.hpp"
15 
16 namespace mfem
17 {
18 
19 #ifdef MFEM_USE_MPI
20 
22  const Vector &X_vert)
23  : face_fes(pfes_ho_),
24  order(face_fes.GetMaxElementOrder()),
25  edge_fec(order, dim),
26  edge_fes(face_fes.GetParMesh(), &edge_fec),
27  ams(edge_fes, X_vert)
28 {
29  MFEM_VERIFY(face_fes.GetParMesh()->Dimension() == dim, "Bad dimension.")
31 }
32 
34 {
35  const int o = order;
36  const int op1 = o + 1;
37  const int nface = dim*o*o*op1;
38 
39  face2edge.SetSize(4*nface);
40  auto f2e = Reshape(face2edge.HostWrite(), 4, nface);
41 
42  for (int c=0; c<dim; ++c)
43  {
44  const int nx = (c == 0) ? op1 : o;
45  const int ny = (c == 1) ? op1 : o;
46 
47  const int c1 = (c+1)%dim;
48  const int c2 = (c+2)%dim;
49 
50  const int nx_e1 = (c1 == 0) ? o : op1;
51  const int ny_e1 = (c1 == 1) ? o : op1;
52  const int nx_e2 = (c2 == 0) ? o : op1;
53  const int ny_e2 = (c2 == 1) ? o : op1;
54 
55  for (int i=0; i<o*o*op1; ++i)
56  {
57  int i1[dim], i2[dim];
58  const int ix = i1[0] = i2[0] = i%nx;
59  const int iy = i1[1] = i2[1] = (i/nx)%ny;
60  const int iz = i1[2] = i2[2] = i/nx/ny;
61 
62  const int iface = ix + iy*nx + iz*nx*ny + c*o*o*op1;
63 
64  ++i1[c2];
65  ++i2[c1];
66 
67  const int ie0 = c1*o*op1*op1 + ix + iy*nx_e1 + iz*nx_e1*ny_e1;
68  const int ie1 = c1*o*op1*op1 + i1[0] + i1[1]*nx_e1 + i1[2]*nx_e1*ny_e1;
69 
70  const int ie2 = c2*o*op1*op1 + ix + iy*nx_e2 + iz*nx_e2*ny_e2;
71  const int ie3 = c2*o*op1*op1 + i2[0] + i2[1]*nx_e2 + i2[2]*nx_e2*ny_e2;
72 
73  f2e(0, iface) = ie0;
74  f2e(1, iface) = ie1;
75  f2e(2, iface) = ie2;
76  f2e(3, iface) = ie3;
77  }
78  }
79 }
80 
82 {
83  // The curl matrix maps from LOR edges to LOR faces. Given a quadrilateral
84  // face (defined by its four edges) f_i = (e_j1, e_j2, e_j3, e_j4), the
85  // matrix has nonzeros A(i, jk), so there are always exactly four nonzeros
86  // per row.
87  const int nface_dof = face_fes.GetNDofs();
88  const int nedge_dof = edge_fes.GetNDofs();
89 
90  SparseMatrix C_local;
91  C_local.OverrideSize(nface_dof, nedge_dof);
92 
93  C_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
94  // Each row always has four nonzeros
95  const int nnz = 4*nedge_dof;
96  auto I = C_local.WriteI();
97  mfem::forall(nedge_dof+1, [=] MFEM_HOST_DEVICE (int i) { I[i] = 4*i; });
98 
99  // face2edge is a mapping of size (4, nface_per_el), such that with a macro
100  // element, face i (in lexicographic ordering) has four edges given by the
101  // entries (k, i), for k=1,2,3,4.
102  Array<int> face2edge;
103  Form3DFaceToEdge(face2edge);
104 
106  const auto *R_f = dynamic_cast<const ElementRestriction*>(
107  face_fes.GetElementRestriction(ordering));
108  const auto *R_e = dynamic_cast<const ElementRestriction*>(
109  edge_fes.GetElementRestriction(ordering));
110  MFEM_VERIFY(R_f != NULL && R_e != NULL, "");
111 
112  const int nel_ho = edge_fes.GetNE();
113  const int nedge_per_el = dim*order*(order+1)*(order+1);
114  const int nface_per_el = dim*order*order*(order+1);
115 
116  const auto offsets_f = R_f->Offsets().Read();
117  const auto indices_f = R_f->Indices().Read();
118  const auto gather_e = Reshape(R_e->GatherMap().Read(), nedge_per_el, nel_ho);
119 
120  const auto f2e = Reshape(face2edge.Read(), 4, nface_per_el);
121 
122  // Fill J and data
123  C_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
125 
126  auto J = C_local.WriteJ();
127  auto V = C_local.WriteData();
128 
129  // Loop over Raviart-Thomas L-DOFs
130  mfem::forall(nface_dof, [=] MFEM_HOST_DEVICE (int i)
131  {
132  const int sj = indices_f[offsets_f[i]]; // signed
133  const int j = (sj >= 0) ? sj : -1 - sj;
134  const int sgn_f = (sj >= 0) ? 1 : -1;
135  const int j_loc = j%nface_per_el;
136  const int j_el = j/nface_per_el;
137 
138  for (int k=0; k<4; ++k)
139  {
140  const int je_loc = f2e(k, j_loc);
141  const int sje = gather_e(je_loc, j_el); // signed
142  const int je = (sje >= 0) ? sje : -1 - sje;
143  const int sgn_e = (sje >= 0) ? 1 : -1;
144  const int sgn = (k == 1 || k == 2) ? -1 : 1;
145  J[i*4 + k] = je;
146  V[i*4 + k] = sgn*sgn_f*sgn_e;
147  }
148  });
149 
150  // Create a block diagonal parallel matrix
152  C_diag.MakeRectangularBlockDiag(edge_fes.GetComm(),
157  &C_local);
158 
159  // Assemble the parallel gradient matrix, must be deleted by the caller
161  {
162  C = C_diag.As<HypreParMatrix>();
163  C_diag.SetOperatorOwner(false);
164  HypreStealOwnership(*C, C_local);
165  }
166  else
167  {
175  Rt.As<SparseMatrix>());
176  C = RAP(Rt_diag.As<HypreParMatrix>(),
177  C_diag.As<HypreParMatrix>(),
179  }
180  C->CopyRowStarts();
181  C->CopyColStarts();
182 }
183 
185 {
186  return StealPointer(C);
187 }
188 
190 {
191  delete C;
192 }
193 
195  ParBilinearForm &a_ho, const Array<int> &ess_tdof_list, int ref_type)
196 {
197  MFEM_VERIFY(a_ho.FESpace()->GetMesh()->Dimension() == 3,
198  "The ADS solver is only valid in 3D.");
200  {
201  ParFiniteElementSpace &pfes = *a_ho.ParFESpace();
202  BatchedLORAssembly batched_lor(pfes);
203  BatchedLOR_ADS lor_ads(pfes, batched_lor.GetLORVertexCoordinates());
204  BatchedLOR_AMS &lor_ams = lor_ads.GetAMS();
205  batched_lor.Assemble(a_ho, ess_tdof_list, A);
206  xyz = lor_ams.StealCoordinateVector();
207  solver = new HypreADS(*A.As<HypreParMatrix>(),
208  lor_ads.StealCurlMatrix(),
209  lor_ams.StealGradientMatrix(),
210  lor_ams.StealXCoordinate(),
211  lor_ams.StealYCoordinate(),
212  lor_ams.StealZCoordinate());
213  }
214  else
215  {
216  ParLORDiscretization lor(a_ho, ess_tdof_list, ref_type);
217  // Assume ownership of the system matrix so that `lor` can be safely
218  // deleted
219  A.Reset(lor.GetAssembledSystem().Ptr());
221  solver = new HypreADS(lor.GetAssembledMatrix(), &lor.GetParFESpace());
222  }
223  width = solver->Width();
224  height = solver->Height();
225 }
226 
228 {
229  solver->SetOperator(op);
230 }
231 
232 void LORSolver<HypreADS>::Mult(const Vector &x, Vector &y) const
233 {
234  solver->Mult(x, y);
235 }
236 
238 
239 const HypreADS &LORSolver<HypreADS>::GetSolver() const { return *solver; }
240 
242 {
243  delete solver;
244  delete xyz;
245 }
246 
247 #endif
248 
249 } // namespace mfem
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
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1152
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1820
LORBase * lor
Definition: lor.hpp:206
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
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
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
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
HypreParMatrix * StealCurlMatrix()
Steal ownership of the discrete curl matrix.
Definition: lor_ads.cpp:184
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:269
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:1302
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:253
ParFiniteElementSpace edge_fes
The associated Nedelec space.
Definition: lor_ads.hpp:37
const OperatorHandle & GetAssembledSystem() const
Returns the assembled LOR system (const version).
Definition: lor.cpp:270
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
Make A_hyp steal ownership of its diagonal part A_diag.
Definition: hypre.cpp:2744
Data type sparse matrix.
Definition: sparsemat.hpp:50
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
SolverType solver
Definition: lor.hpp:208
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:241
HypreParMatrix * C
The discrete curl matrix.
Definition: lor_ads.hpp:39
ParFiniteElementSpace & face_fes
The RT space.
Definition: lor_ads.hpp:33
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3232
static constexpr int dim
Spatial dimension, always 3.
Definition: lor_ads.hpp:34
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:736
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
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
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
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:37
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:279
void FormCurlMatrix()
Form the discrete curl matrix (not part of the public API).
Definition: lor_ads.cpp:81
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void Assemble(BilinearForm &a, const Array< int > ess_dofs, OperatorHandle &A)
Assemble the given form as a matrix and place the result in A.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
BatchedLOR_ADS(ParFiniteElementSpace &pfes_ho_, const Vector &X_vert)
Construct the BatchedLOR_AMS object associated with the 3D RT space pfes_ho_.
Definition: lor_ads.cpp:21
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:277
T * HostWrite()
Shortcut for mfem::Write(a.GetMemory(), a.Size(), false).
Definition: array.hpp:319
const int order
Polynomial degree.
Definition: lor_ads.hpp:35
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
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
int dim
Definition: ex24.cpp:53
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Definition: sparsemat.cpp:297
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:269
Vector data type.
Definition: vector.hpp:58
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:317
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
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:278
void Form3DFaceToEdge(Array< int > &face2edge)
Form the local elementwise discrete curl matrix.
Definition: lor_ads.cpp:33
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