MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lor_batched.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #ifndef MFEM_LOR_BATCHED
13 #define MFEM_LOR_BATCHED
14 
15 #include "lor.hpp"
16 #include "../qspace.hpp"
17 
18 namespace mfem
19 {
20 
21 /// @brief Efficient batched assembly of LOR discretizations on device.
22 ///
23 /// This class should typically be used by the user-facing classes
24 /// LORDiscretization and ParLORDiscretization. Only certain bilinear forms are
25 /// supported, currently:
26 ///
27 /// - H1 diffusion + mass
28 /// - ND curl-curl + mass
29 /// - RT div-div + mass
30 ///
31 /// Whether a form is supported can be checked with the static member function
32 /// BatchedLORAssembly::FormIsSupported.
34 {
35 protected:
36  FiniteElementSpace &fes_ho; ///< The high-order space.
37 
38  Vector X_vert; ///< LOR vertex coordinates.
39 
40  /// @brief The elementwise LOR matrices in a sparse "ij" format.
41  ///
42  /// This is interpreted to have shape (nnz_per_row, ndof_per_el, nel_ho). For
43  /// index (i, j, k), this represents row @a j of the @a kth element matrix.
44  /// The column index is given by sparse_mapping(i, j).
46 
47  /// @brief The sparsity pattern of the element matrices.
48  ///
49  /// This array should be interpreted as having shape (nnz_per_row,
50  /// ndof_per_el). For local DOF index @a j, sparse_mapping(i, j) is the
51  /// column index of the @a ith nonzero in the @a jth row. If the index is
52  /// negative, that entry should be skipped (there is no corresponding
53  /// nonzero).
55 
56 public:
57  /// Construct the batched assembly object corresponding to @a fes_ho_.
59 
60  /// Returns true if the form @a a supports batched assembly, false otherwise.
61  static bool FormIsSupported(BilinearForm &a);
62 
63  /// @brief Assemble the given form as a matrix and place the result in @a A.
64  ///
65  /// In serial, the result will be a SparseMatrix. In parallel, the result
66  /// will be a HypreParMatrix.
67  void Assemble(BilinearForm &a, const Array<int> ess_dofs, OperatorHandle &A);
68 
69  /// Compute the vertices of the LOR mesh and place the result in @a X_vert.
71  Vector &X_vert);
72 
73  /// Return the vertices of the LOR mesh in E-vector format
74  const Vector &GetLORVertexCoordinates() { return X_vert; }
75 
76 protected:
77  /// After assembling the "sparse IJ" format, convert it to CSR.
78  void SparseIJToCSR(OperatorHandle &A) const;
79 
80  /// Assemble the system without eliminating essential DOFs.
82 
83  /// @brief Fill in @a sparse_ij and @a sparse_mapping using one of the
84  /// specialized LOR assembly kernel classes.
85  ///
86  /// @sa Specialization classes: BatchedLOR_H1, BatchedLOR_ND, BatchedLOR_RT
87  template <typename LOR_KERNEL> void AssemblyKernel(BilinearForm &a);
88 
89 public:
90  /// @name GPU kernel functions
91  /// These functions should be considered protected, but they contain
92  /// MFEM_FORALL kernels, and so they must be public (this is a compiler
93  /// limitation).
94  ///@{
95 
96  /// @brief Fill the I array of the sparse matrix @a A.
97  ///
98  /// @note AssemblyKernel must be called first to populate @a sparse_mapping.
99  int FillI(SparseMatrix &A) const;
100 
101  /// @brief Fill the J and data arrays of the sparse matrix @a A.
102  ///
103  /// @note AssemblyKernel must be called first to populate @a sparse_mapping
104  /// and @a sparse_ij.
105  void FillJAndData(SparseMatrix &A) const;
106 
107 #ifdef MFEM_USE_MPI
108  /// Assemble the system in parallel and place the result in @a A.
109  void ParAssemble(BilinearForm &a, const Array<int> &ess_dofs,
110  OperatorHandle &A);
111 #endif
112  ///@}
113 };
114 
115 /// @brief Ensure that @a mem has at least capacity @a capacity.
116 ///
117 /// If the capacity of @a mem is not large enough, delete it and allocate new
118 /// memory with size @a capacity.
119 template <typename T>
120 void EnsureCapacity(Memory<T> &mem, int capacity)
121 {
122  if (mem.Capacity() < capacity)
123  {
124  mem.Delete();
125  mem.New(capacity, mem.GetMemoryType());
126  }
127 }
128 
129 /// Return the first domain integrator in the form @a i of type @a T.
130 template <typename T>
131 static T *GetIntegrator(BilinearForm &a)
132 {
133  Array<BilinearFormIntegrator*> *integs = a.GetDBFI();
134  if (integs != NULL)
135  {
136  for (auto *i : *integs)
137  {
138  if (auto *ti = dynamic_cast<T*>(i))
139  {
140  return ti;
141  }
142  }
143  }
144  return nullptr;
145 }
146 
147 IntegrationRule GetCollocatedIntRule(FiniteElementSpace &fes);
148 
149 template <typename INTEGRATOR>
151 {
152  INTEGRATOR *i = GetIntegrator<INTEGRATOR>(a);
153  if (i)
154  {
155  // const_cast since Coefficient::Eval is not const...
156  auto *coeff = const_cast<Coefficient*>(i->GetCoefficient());
157  if (coeff) { coeff_vector.Project(*coeff); }
158  else { coeff_vector.SetConstant(1.0); }
159  }
160  else
161  {
162  coeff_vector.SetConstant(0.0);
163  }
164 }
165 
166 /// Abstract base class for the batched LOR assembly kernels.
168 {
169 protected:
170  FiniteElementSpace &fes_ho; ///< The associated high-order space.
171  Vector &X_vert; ///< Mesh coordinate vector.
172  Vector &sparse_ij; ///< Local element sparsity matrix data.
173  Array<int> &sparse_mapping; ///< Local element sparsity pattern.
174  IntegrationRule ir; ///< Collocated integration rule.
175  QuadratureSpace qs; ///< Quadrature space for coefficients.
176  CoefficientVector c1; ///< Coefficient of first integrator.
177  CoefficientVector c2; ///< Coefficient of second integrator.
179  Vector &X_vert_,
180  Vector &sparse_ij_,
181  Array<int> &sparse_mapping_)
182  : fes_ho(fes_ho_), X_vert(X_vert_), sparse_ij(sparse_ij_),
183  sparse_mapping(sparse_mapping_), ir(GetCollocatedIntRule(fes_ho)),
186  { }
187 };
188 
189 }
190 
191 #endif
void ParAssemble(BilinearForm &a, const Array< int > &ess_dofs, OperatorHandle &A)
Assemble the system in parallel and place the result in A.
FiniteElementSpace & fes_ho
The high-order space.
Definition: lor_batched.hpp:36
CoefficientStorage
Flags that determine what storage optimizations to use in CoefficientVector.
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Abstract base class for the batched LOR assembly kernels.
Array< int > & sparse_mapping
Local element sparsity pattern.
void ProjectLORCoefficient(BilinearForm &a, CoefficientVector &coeff_vector)
void Delete()
Delete the owned pointers and reset the Memory object.
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
BatchedLORKernel(FiniteElementSpace &fes_ho_, Vector &X_vert_, Vector &sparse_ij_, Array< int > &sparse_mapping_)
void SetConstant(double constant)
Set this vector to the given constant.
CoefficientVector c2
Coefficient of second integrator.
int Capacity() const
Return the size of the allocated memory.
static bool FormIsSupported(BilinearForm &a)
Returns true if the form a supports batched assembly, false otherwise.
Definition: lor_batched.cpp:49
Class to represent a coefficient evaluated at quadrature points.
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
static void FormLORVertexCoordinates(FiniteElementSpace &fes_ho, Vector &X_vert)
Compute the vertices of the LOR mesh and place the result in X_vert.
Definition: lor_batched.cpp:72
void Project(Coefficient &coeff)
Evaluate the given Coefficient at the quadrature points defined by qs.
Data type sparse matrix.
Definition: sparsemat.hpp:50
IntegrationRule ir
Collocated integration rule.
Vector sparse_ij
The elementwise LOR matrices in a sparse &quot;ij&quot; format.
Definition: lor_batched.hpp:45
Enable all above compressions.
Vector & sparse_ij
Local element sparsity matrix data.
Vector X_vert
LOR vertex coordinates.
Definition: lor_batched.hpp:38
Array< int > sparse_mapping
The sparsity pattern of the element matrices.
Definition: lor_batched.hpp:54
FiniteElementSpace & fes_ho
The associated high-order space.
int FillI(SparseMatrix &A) const
Fill the I array of the sparse matrix A.
void FillJAndData(SparseMatrix &A) const
Fill the J and data arrays of the sparse matrix A.
QuadratureSpace qs
Quadrature space for coefficients.
const Vector & GetLORVertexCoordinates()
Return the vertices of the LOR mesh in E-vector format.
Definition: lor_batched.hpp:74
Vector & X_vert
Mesh coordinate vector.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void Assemble(BilinearForm &a, const Array< int > ess_dofs, OperatorHandle &A)
Assemble the given form as a matrix and place the result in A.
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
void EnsureCapacity(Memory< T > &mem, int capacity)
Ensure that mem has at least capacity capacity.
Class used by MFEM to store pointers to host and/or device memory.
void AssembleWithoutBC(BilinearForm &a, OperatorHandle &A)
Assemble the system without eliminating essential DOFs.
void SparseIJToCSR(OperatorHandle &A) const
After assembling the &quot;sparse IJ&quot; format, convert it to CSR.
Vector data type.
Definition: vector.hpp:60
BatchedLORAssembly(FiniteElementSpace &fes_ho_)
Construct the batched assembly object corresponding to fes_ho_.
Class representing the storage layout of a QuadratureFunction.
Definition: qspace.hpp:92
CoefficientVector c1
Coefficient of first integrator.
IntegrationRule GetCollocatedIntRule(FiniteElementSpace &fes)
void AssemblyKernel(BilinearForm &a)
Fill in sparse_ij and sparse_mapping using one of the specialized LOR assembly kernel classes...