MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
lor_ads.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_ads.hpp"
15
16namespace 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*>(
108 const auto *R_e = dynamic_cast<const ElementRestriction*>(
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
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
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
233{
234 solver->Mult(x, y);
235}
236
238
239const HypreADS &LORSolver<HypreADS>::GetSolver() const { return *solver; }
240
242{
243 delete solver;
244 delete xyz;
245}
246
247#endif
248
249} // 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.
ParFiniteElementSpace & face_fes
The RT space.
Definition lor_ads.hpp:33
HypreParMatrix * StealCurlMatrix()
Steal ownership of the discrete curl matrix.
Definition lor_ads.cpp:184
BatchedLOR_AMS & GetAMS()
Return the associated BatchedLOR_AMS object.
Definition lor_ads.hpp:61
HypreParMatrix * C
The discrete curl matrix.
Definition lor_ads.hpp:39
void FormCurlMatrix()
Form the discrete curl matrix (not part of the public API).
Definition lor_ads.cpp:81
ParFiniteElementSpace edge_fes
The associated Nedelec space.
Definition lor_ads.hpp:37
const int order
Polynomial degree.
Definition lor_ads.hpp:35
static constexpr int dim
Spatial dimension, always 3.
Definition lor_ads.hpp:34
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
void Form3DFaceToEdge(Array< int > &face2edge)
Form the local elementwise discrete curl matrix.
Definition lor_ads.cpp:33
HypreParVector * StealZCoordinate()
Definition lor_ams.cpp:343
HypreParVector * StealYCoordinate()
Definition lor_ams.cpp:338
HypreParMatrix * StealGradientMatrix()
Definition lor_ams.cpp:323
Vector * StealCoordinateVector()
Definition lor_ams.cpp:328
HypreParVector * StealXCoordinate()
Definition lor_ams.cpp:333
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
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.
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
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
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
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
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
@ 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
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
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
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
int dim
Definition ex24.cpp:53
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
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
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