MFEM  v4.5.2
Finite element discretization library
tmop_pa.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 "../tmop.hpp"
13 #include "../linearform.hpp"
14 #include "../pgridfunc.hpp"
15 #include "../tmop_tools.hpp"
16 #include "../quadinterpolator.hpp"
17 #include "../../general/forall.hpp"
18 #include "../../linalg/kernels.hpp"
19 
20 namespace mfem
21 {
22 
24  const FiniteElementSpace &fes)
25 {
26  MFEM_VERIFY(PA.enabled, "PA extension setup has not been done!");
27  MFEM_VERIFY(PA.fes == &fes, "");
28  // TODO: we need a more robust way to check that the 'fes' used when
29  // AssemblePA() was called has not been modified or completely destroyed and
30  // a new object created at the same address.
31 
32  if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
33  {
35  PA.Jtr_debug_grad = true;
36  }
37 
38  if (PA.dim == 2)
39  {
41  if (lim_coeff) { AssembleGradPA_C0_2D(xe); }
42  }
43 
44  if (PA.dim == 3)
45  {
47  if (lim_coeff) { AssembleGradPA_C0_3D(xe); }
48  }
49 }
50 
52 {
53  const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
55  // Return immediately if limiting is not enabled
56  if (lim_coeff == nullptr) { return; }
57  MFEM_VERIFY(lim_nodes0, "internal error");
58 
59  MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
60  MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
61  MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func) ||
62  dynamic_cast<TMOP_ExponentialLimiter*>(lim_func),
63  "Only TMOP_QuadraticLimiter and TMOP_ExponentialLimiter are supported");
64 
65  const FiniteElementSpace *fes = PA.fes;
66  const int NE = PA.ne;
67  if (NE == 0) { return; } // Quick return for empty processors
68  const IntegrationRule &ir = *PA.ir;
69 
71 
72  // H0 for lim_coeff, (dim x dim) Q-vector
73  PA.H0.UseDevice(true);
74  PA.H0.SetSize(PA.dim * PA.dim * PA.nq * NE, mt);
75 
76  // lim_coeff -> PA.C0 (Q-vector)
77  PA.C0.UseDevice(true);
78  if (ConstantCoefficient* cQ =
79  dynamic_cast<ConstantCoefficient*>(lim_coeff))
80  {
81  PA.C0.SetSize(1, Device::GetMemoryType());
82  PA.C0.HostWrite();
83  PA.C0(0) = cQ->constant;
84  }
85  else
86  {
87  PA.C0.SetSize(PA.nq * PA.ne, Device::GetMemoryType());
88  auto C0 = Reshape(PA.C0.HostWrite(), PA.nq, PA.ne);
89  for (int e = 0; e < NE; ++e)
90  {
92  for (int q = 0; q < ir.GetNPoints(); ++q)
93  {
94  C0(q,e) = lim_coeff->Eval(T, ir.IntPoint(q));
95  }
96  }
97  }
98 
99  // lim_nodes0 -> PA.X0 (E-vector)
100  MFEM_VERIFY(lim_nodes0->FESpace() == fes, "");
101  const Operator *n0_R = fes->GetElementRestriction(ordering);
102  PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
103  PA.X0.UseDevice(true);
104  n0_R->Mult(*lim_nodes0, PA.X0);
105 
106  // Limiting distances: lim_dist -> PA.LD (E-vector)
107  // TODO: remove the hack for the case lim_dist == NULL.
108  const FiniteElementSpace *limfes = (lim_dist) ? lim_dist->FESpace() : fes;
109  const FiniteElement &lim_fe = *limfes->GetFE(0);
110  PA.maps_lim = &lim_fe.GetDofToQuad(ir, DofToQuad::TENSOR);
111  PA.LD.SetSize(NE*lim_fe.GetDof(), Device::GetMemoryType());
112  PA.LD.UseDevice(true);
113  if (lim_dist)
114  {
115  const Operator *ld_R = limfes->GetElementRestriction(ordering);
116  ld_R->Mult(*lim_dist, PA.LD);
117  }
118  else
119  {
120  PA.LD = 1.0;
121  }
122 }
123 
125  const IntegrationRule &ir,
126  const Vector &xe,
127  DenseTensor &Jtr) const
128 {
129  MFEM_VERIFY(Jtr.SizeI() == Jtr.SizeJ() && Jtr.SizeI() > 1, "");
130  const int dim = Jtr.SizeI();
131  bool done = false;
132  if (dim == 2) { done = ComputeAllElementTargets<2>(fes, ir, xe, Jtr); }
133  if (dim == 3) { done = ComputeAllElementTargets<3>(fes, ir, xe, Jtr); }
134 
135  if (!done) { ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr); }
136 }
137 
139  const IntegrationRule &ir,
140  const Vector &xe,
141  DenseTensor &Jtr) const
142 {
143  ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr);
144 }
145 
146 
147 // Code paths leading to ComputeElementTargets:
148 // - GetElementEnergy(elfun) which is done through GetGridFunctionEnergyPA(x)
149 // - AssembleElementVectorExact(elfun)
150 // - AssembleElementGradExact(elfun)
151 // - EnableNormalization(x) -> ComputeNormalizationEnergies(x)
152 // - (AssembleElementVectorFD(elfun))
153 // - (AssembleElementGradFD(elfun))
154 // ============================================================================
155 // - TargetConstructor():
156 // - IDEAL_SHAPE_UNIT_SIZE: Wideal
157 // - IDEAL_SHAPE_EQUAL_SIZE: α * Wideal
158 // - IDEAL_SHAPE_GIVEN_SIZE: β * Wideal
159 // - GIVEN_SHAPE_AND_SIZE: β * Wideal
160 // - AnalyticAdaptTC(elfun):
161 // - GIVEN_FULL: matrix_tspec->Eval(Jtr(elfun))
162 // - DiscreteAdaptTC():
163 // - IDEAL_SHAPE_GIVEN_SIZE: size^{1.0/dim} * Jtr(i) (size)
164 // - GIVEN_SHAPE_AND_SIZE: Jtr(i) *= D_rho (ratio)
165 // Jtr(i) *= Q_phi (skew)
166 // Jtr(i) *= R_theta (orientation)
168 {
169  PA.Jtr_needs_update = false;
170  PA.Jtr_debug_grad = false;
171  const FiniteElementSpace *fes = PA.fes;
172  if (PA.ne == 0) { return; } // Quick return for empty processors
173  const IntegrationRule &ir = *PA.ir;
174 
175  // Compute PA.Jtr for all elements
177 }
178 
180 {
181  const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
183  PA.enabled = true;
184  PA.fes = &fes;
185  Mesh *mesh = fes.GetMesh();
186  const int ne = PA.ne = mesh->GetNE();
187  if (ne == 0) { return; } // Quick return for empty processors
188  const int dim = PA.dim = mesh->Dimension();
189  MFEM_VERIFY(PA.dim == 2 || PA.dim == 3, "Not yet implemented!");
190  MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
191  "mixed meshes are not supported");
192  MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
193  const FiniteElement &fe = *fes.GetFE(0);
194  PA.ir = &EnergyIntegrationRule(fe);
195  const IntegrationRule &ir = *PA.ir;
196  MFEM_VERIFY(fes.GetOrdering() == Ordering::byNODES,
197  "PA Only supports Ordering::byNODES!");
198 
199  const int nq = PA.nq = ir.GetNPoints();
200  const DofToQuad::Mode mode = DofToQuad::TENSOR;
201  PA.maps = &fe.GetDofToQuad(ir, mode);
203 
204  // Energy vector, scalar Q-vector
205  PA.E.UseDevice(true);
206  PA.E.SetSize(ne*nq, Device::GetDeviceMemoryType());
207 
208  // H for Grad, (dim x dim) Q-vector
209  PA.H.UseDevice(true);
210  PA.H.SetSize(dim*dim * dim*dim * nq*ne, mt);
211 
212  // Scalar Q-vector of '1', used to compute sums via dot product
213  PA.O.SetSize(ne*nq, Device::GetDeviceMemoryType());
214  PA.O = 1.0;
215 
216  // Setup ref->target Jacobians, PA.Jtr, (dim x dim) Q-vector, DenseTensor
217  PA.Jtr.SetSize(dim, dim, PA.ne*PA.nq, mt);
218  PA.Jtr_needs_update = true;
219  PA.Jtr_debug_grad = false;
220 
221  // Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
222  if (lim_coeff) { AssemblePA_Limiting(); }
223 }
224 
226 {
227  // This method must be called after AssembleGradPA().
228 
229  MFEM_VERIFY(PA.Jtr_needs_update == false, "");
230 
232  {
233  MFEM_VERIFY(PA.Jtr_debug_grad == true, "AssembleGradPA() was not called"
234  " or Jtr was overwritten by another method!");
235  }
236 
237  if (PA.dim == 2)
238  {
240  if (lim_coeff) { AssembleDiagonalPA_C0_2D(de); }
241  }
242 
243  if (PA.dim == 3)
244  {
246  if (lim_coeff) { AssembleDiagonalPA_C0_3D(de); }
247  }
248 }
249 
250 void TMOP_Integrator::AddMultPA(const Vector &xe, Vector &ye) const
251 {
252  // This method must be called after AssemblePA().
253 
254  if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
255  {
257  }
258 
259  if (PA.dim == 2)
260  {
261  AddMultPA_2D(xe,ye);
262  if (lim_coeff) { AddMultPA_C0_2D(xe,ye); }
263  }
264 
265  if (PA.dim == 3)
266  {
267  AddMultPA_3D(xe,ye);
268  if (lim_coeff) { AddMultPA_C0_3D(xe,ye); }
269  }
270 }
271 
272 void TMOP_Integrator::AddMultGradPA(const Vector &re, Vector &ce) const
273 {
274  // This method must be called after AssembleGradPA().
275 
276  MFEM_VERIFY(PA.Jtr_needs_update == false, "");
277 
279  {
280  MFEM_VERIFY(PA.Jtr_debug_grad == true, "AssembleGradPA() was not called or "
281  "Jtr was overwritten by another method!");
282  }
283 
284  if (PA.dim == 2)
285  {
286  AddMultGradPA_2D(re,ce);
287  if (lim_coeff) { AddMultGradPA_C0_2D(re,ce); }
288  }
289 
290  if (PA.dim == 3)
291  {
292  AddMultGradPA_3D(re,ce);
293  if (lim_coeff) { AddMultGradPA_C0_3D(re,ce); }
294  }
295 }
296 
298 {
299  // This method must be called after AssemblePA().
300 
301  double energy = 0.0;
302 
303  if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
304  {
306  }
307 
308  if (PA.dim == 2)
309  {
310  energy = GetLocalStateEnergyPA_2D(xe);
311  if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
312  }
313 
314  if (PA.dim == 3)
315  {
316  energy = GetLocalStateEnergyPA_3D(xe);
317  if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
318  }
319 
320  return energy;
321 }
322 
323 } // namespace mfem
Abstract class for all finite elements.
Definition: fe_base.hpp:232
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:463
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:160
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:840
const GridFunction * lim_dist
Definition: tmop.hpp:1663
void AddMultGradPA_2D(const Vector &, Vector &) const
bool UsesPhysicalCoordinates() const
Return true if the methods ComputeElementTargets(), ComputeAllElementTargets(), and ComputeElementTar...
Definition: tmop.hpp:1327
struct mfem::TMOP_Integrator::@23 PA
int Dimension() const
Definition: mesh.hpp:1047
void AssembleGradPA_2D(const Vector &) const
void ComputeAllElementTargets(const Vector &xe=Vector()) const
Definition: tmop_pa.cpp:167
void AddMultPA_3D(const Vector &, Vector &) const
Definition: tmop_pa_p3.cpp:192
int SizeJ() const
Definition: densemat.hpp:1025
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const IntegrationRule * ir
Definition: tmop.hpp:1749
virtual double GetLocalStateEnergyPA(const Vector &) const
Compute the local (to the MPI rank) energy with partial assembly.
Definition: tmop_pa.cpp:297
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5921
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1259
const IntegrationRule & EnergyIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1810
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
Definition: tmop_pa.cpp:250
void AddMultPA_C0_2D(const Vector &, Vector &) const
bool ComputeAllElementTargets(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
double GetLocalStateEnergyPA_C0_3D(const Vector &) const
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:364
void AssembleGradPA_C0_2D(const Vector &) const
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
virtual void AssembleGradPA(const Vector &, const FiniteElementSpace &)
Prepare the integrator for partial assembly (PA) gradient evaluations on the given FE space fes at th...
Definition: tmop_pa.cpp:23
void AddMultGradPA_3D(const Vector &, Vector &) const
void AssembleDiagonalPA_3D(Vector &) const
void AddMultGradPA_C0_2D(const Vector &, Vector &) const
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition: device.hpp:277
const TargetConstructor * targetC
Definition: tmop.hpp:1646
const GridFunction * lim_nodes0
Definition: tmop.hpp:1660
virtual void AssemblePA(const FiniteElementSpace &)
Method defining partial assembly.
Definition: tmop_pa.cpp:179
void AssembleGradPA_C0_3D(const Vector &) const
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
void AssembleDiagonalPA_2D(Vector &) const
void AssemblePA_Limiting()
Definition: tmop_pa.cpp:51
double GetLocalStateEnergyPA_3D(const Vector &) const
Definition: tmop_pa_w3.cpp:155
Coefficient * lim_coeff
Definition: tmop.hpp:1661
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void AddMultPA_C0_3D(const Vector &, Vector &) const
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
void AssembleDiagonalPA_C0_2D(Vector &) const
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
virtual void AssembleGradDiagonalPA(Vector &) const
Method for computing the diagonal of the gradient with partial assembly.
Definition: tmop_pa.cpp:225
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
double GetLocalStateEnergyPA_C0_2D(const Vector &) const
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition: fe_base.hpp:148
const FiniteElementSpace * fes
Definition: tmop.hpp:1748
void AssembleDiagonalPA_C0_3D(Vector &) const
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
virtual void AddMultGradPA(const Vector &, Vector &) const
Method for partially assembled gradient action.
Definition: tmop_pa.cpp:272
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
int dim
Definition: ex24.cpp:53
int SizeI() const
Definition: densemat.hpp:1024
Lexicographic ordering for tensor-product FiniteElements.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void AssembleGradPA_3D(const Vector &) const
double GetLocalStateEnergyPA_2D(const Vector &) const
Definition: tmop_pa_w2.cpp:145
Vector data type.
Definition: vector.hpp:60
virtual void ComputeAllElementTargets(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Computes reference-to-target transformation Jacobians for all quadrature points in all elements...
Definition: tmop_pa.cpp:138
Abstract operator.
Definition: operator.hpp:24
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 AddMultPA_2D(const Vector &, Vector &) const
Definition: tmop_pa_p2.cpp:167
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:978
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:1665
void ComputeAllElementTargets_Fallback(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Definition: tmop.cpp:1468
void AddMultGradPA_C0_3D(const Vector &, Vector &) const
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647