MFEM  v4.6.0
Finite element discretization library
nonlinearform_ext.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 // Implementations of classes FABilinearFormExtension, EABilinearFormExtension,
13 // PABilinearFormExtension and MFBilinearFormExtension.
14 
15 #include "nonlinearform.hpp"
16 #include "ceed/interface/util.hpp"
17 
18 namespace mfem
19 {
20 
22  : Operator(nlf->FESpace()->GetVSize()), nlf(nlf) { }
23 
26  fes(*nlf->FESpace()),
27  dnfi(*nlf->GetDNFI()),
28  elemR(nullptr),
29  Grad(*this)
30 {
31  if (!DeviceCanUseCeed())
32  {
34  // TODO: optimize for the case when 'elemR' is identity
37  }
38  ye.UseDevice(true);
39 }
40 
42 {
43  double energy = 0.0;
44 
45  elemR->Mult(x, xe);
46  for (int i = 0; i < dnfi.Size(); i++)
47  {
48  energy += dnfi[i]->GetLocalStateEnergyPA(xe);
49  }
50  return energy;
51 }
52 
54 {
55  MFEM_VERIFY(nlf->GetInteriorFaceIntegrators().Size() == 0 &&
56  nlf->GetBdrFaceIntegrators().Size() == 0,
57  "face integrators are not supported yet");
58 
59  for (int i = 0; i < dnfi.Size(); ++i) { dnfi[i]->AssemblePA(fes); }
60 }
61 
63 {
64  if (!DeviceCanUseCeed())
65  {
66  ye = 0.0;
67  elemR->Mult(x, xe);
68  for (int i = 0; i < dnfi.Size(); ++i) { dnfi[i]->AddMultPA(xe, ye); }
69  elemR->MultTranspose(ye, y);
70  }
71  else
72  {
73  y.UseDevice(true); // typically this is a large vector, so store on device
74  y = 0.0;
75  for (int i = 0; i < dnfi.Size(); ++i)
76  {
77  dnfi[i]->AddMultPA(x, y);
78  }
79  }
80 }
81 
83 {
84  Grad.AssembleGrad(x);
85  return Grad;
86 }
87 
89 {
90  height = width = fes.GetVSize();
92  xe.SetSize(elemR->Height());
93  ye.SetSize(elemR->Height());
94  Grad.Update();
95 }
96 
97 PANonlinearFormExtension::Gradient::Gradient(const PANonlinearFormExtension &e):
98  Operator(e.Height()), ext(e)
99 { }
100 
101 void PANonlinearFormExtension::Gradient::AssembleGrad(const Vector &g)
102 {
103  ext.elemR->Mult(g, ext.xe);
104  for (int i = 0; i < ext.dnfi.Size(); ++i)
105  {
106  ext.dnfi[i]->AssembleGradPA(ext.xe, ext.fes);
107  }
108 }
109 
110 void PANonlinearFormExtension::Gradient::Mult(const Vector &x, Vector &y) const
111 {
112  ext.ye = 0.0;
113  ext.elemR->Mult(x, ext.xe);
114  for (int i = 0; i < ext.dnfi.Size(); ++i)
115  {
116  ext.dnfi[i]->AddMultGradPA(ext.xe, ext.ye);
117  }
118  ext.elemR->MultTranspose(ext.ye, y);
119 }
120 
121 void PANonlinearFormExtension::Gradient::AssembleDiagonal(Vector &diag) const
122 {
123  MFEM_ASSERT(diag.Size() == Height(),
124  "Vector for holding diagonal has wrong size!");
125  ext.ye = 0.0;
126  for (int i = 0; i < ext.dnfi.Size(); ++i)
127  {
128  ext.dnfi[i]->AssembleGradDiagonalPA(ext.ye);
129  }
130  ext.elemR->MultTranspose(ext.ye, diag);
131 }
132 
134 {
135  height = width = ext.Height();
136 }
137 
138 
140  NonlinearFormExtension(form), fes(*form->FESpace())
141 {
142  if (!DeviceCanUseCeed())
143  {
146  if (elem_restrict_lex) // replace with a check for not identity
147  {
150  localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
151  }
152  }
153 }
154 
156 {
157  const Array<NonlinearFormIntegrator*> &integrators = *nlf->GetDNFI();
158  const int Ni = integrators.Size();
159  for (int i = 0; i < Ni; ++i)
160  {
161  integrators[i]->AssembleMF(fes);
162  }
163 }
164 
166 {
167  const Array<NonlinearFormIntegrator*> &integrators = *nlf->GetDNFI();
168  const int iSz = integrators.Size();
169  // replace the check 'elem_restrict_lex' with a check for not identity
171  {
173  localY = 0.0;
174  for (int i = 0; i < iSz; ++i)
175  {
176  integrators[i]->AddMultMF(localX, localY);
177  }
179  }
180  else
181  {
182  y.UseDevice(true); // typically this is a large vector, so store on device
183  y = 0.0;
184  for (int i = 0; i < iSz; ++i)
185  {
186  integrators[i]->AddMultMF(x, y);
187  }
188  }
189 }
190 
192 {
193  height = width = fes.GetVSize();
196  if (elem_restrict_lex) // replace with a check for not identity
197  {
200  }
201 }
202 
203 } // namespace mfem
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
const Array< NonlinearFormIntegrator * > & GetBdrFaceIntegrators() const
Access all boundary face integrators added with AddBdrFaceIntegrator().
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
const Array< NonlinearFormIntegrator * > & dnfi
NonlinearFormExtension(const NonlinearForm *)
void Mult(const Vector &x, Vector &y) const override
Perform the action of the MFNonlinearFormExtension.
void Assemble() override
Prepare the PANonlinearFormExtension for evaluation with Mult().
const FiniteElementSpace & fes
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Data and methods for partially-assembled nonlinear forms.
void Mult(const Vector &x, Vector &y) const override
Perform the action of the PANonlinearFormExtension.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1302
const NonlinearForm * nlf
Not owned.
void Update() override
Called by NonlinearForm::Update() to reflect changes in the FE space.
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector *> &v)
Definition: solvers.cpp:959
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition: device.hpp:277
void Update() override
Called by NonlinearForm::Update() to reflect changes in the FE space.
Class extending the NonlinearForm class to support the different AssemblyLevels.
double GetGridFunctionEnergy(const Vector &x) const override
Compute the local (to the MPI rank) energy of the L-vector state x.
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
Definition: util.cpp:33
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void Assemble() override
Prepare the MFNonlinearFormExtension for evaluation with Mult().
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
const FiniteElementSpace & fes
MFNonlinearFormExtension(const NonlinearForm *)
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
const Array< NonlinearFormIntegrator * > & GetInteriorFaceIntegrators() const
Access all interior face integrators added with AddInteriorFaceIntegrator().
PANonlinearFormExtension(const NonlinearForm *nlf)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
Lexicographic ordering for tensor-product FiniteElements.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:58
Operator & GetGradient(const Vector &x) const override
Return the gradient as an L-to-L Operator. The input x must be an L-vector (i.e. GridFunction-size ve...
Abstract operator.
Definition: operator.hpp:24
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28