MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nonlinearform_ext.cpp
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 // Implementations of classes FABilinearFormExtension, EABilinearFormExtension,
13 // PABilinearFormExtension and MFBilinearFormExtension.
14 
15 #include "nonlinearform.hpp"
16 #include "ceed/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(fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC)),
29  Grad(*this)
30 {
31  // TODO: optimize for the case when 'elemR' is identity
34  ye.UseDevice(true);
35 }
36 
38 {
39  double energy = 0.0;
40 
41  elemR->Mult(x, xe);
42  for (int i = 0; i < dnfi.Size(); i++)
43  {
44  energy += dnfi[i]->GetLocalStateEnergyPA(xe);
45  }
46  return energy;
47 }
48 
50 {
51  MFEM_VERIFY(nlf->GetInteriorFaceIntegrators().Size() == 0 &&
52  nlf->GetBdrFaceIntegrators().Size() == 0,
53  "face integrators are not supported yet");
54 
55  for (int i = 0; i < dnfi.Size(); ++i) { dnfi[i]->AssemblePA(fes); }
56 }
57 
59 {
60  if (!DeviceCanUseCeed())
61  {
62  ye = 0.0;
63  elemR->Mult(x, xe);
64  for (int i = 0; i < dnfi.Size(); ++i) { dnfi[i]->AddMultPA(xe, ye); }
65  elemR->MultTranspose(ye, y);
66  }
67  else
68  {
69  y.UseDevice(true); // typically this is a large vector, so store on device
70  y = 0.0;
71  for (int i = 0; i < dnfi.Size(); ++i)
72  {
73  dnfi[i]->AddMultPA(x, y);
74  }
75  }
76 }
77 
79 {
80  Grad.AssembleGrad(x);
81  return Grad;
82 }
83 
85 {
86  height = width = fes.GetVSize();
88  xe.SetSize(elemR->Height());
89  ye.SetSize(elemR->Height());
90  Grad.Update();
91 }
92 
93 PANonlinearFormExtension::Gradient::Gradient(const PANonlinearFormExtension &e):
94  Operator(e.Height()), ext(e)
95 { }
96 
97 void PANonlinearFormExtension::Gradient::AssembleGrad(const Vector &g)
98 {
99  ext.elemR->Mult(g, ext.xe);
100  for (int i = 0; i < ext.dnfi.Size(); ++i)
101  {
102  ext.dnfi[i]->AssembleGradPA(ext.xe, ext.fes);
103  }
104 }
105 
106 void PANonlinearFormExtension::Gradient::Mult(const Vector &x, Vector &y) const
107 {
108  ext.ye = 0.0;
109  ext.elemR->Mult(x, ext.xe);
110  for (int i = 0; i < ext.dnfi.Size(); ++i)
111  {
112  ext.dnfi[i]->AddMultGradPA(ext.xe, ext.ye);
113  }
114  ext.elemR->MultTranspose(ext.ye, y);
115 }
116 
117 void PANonlinearFormExtension::Gradient::AssembleDiagonal(Vector &diag) const
118 {
119  MFEM_ASSERT(diag.Size() == Height(),
120  "Vector for holding diagonal has wrong size!");
121  ext.ye = 0.0;
122  for (int i = 0; i < ext.dnfi.Size(); ++i)
123  {
124  ext.dnfi[i]->AssembleGradDiagonalPA(ext.ye);
125  }
126  ext.elemR->MultTranspose(ext.ye, diag);
127 }
128 
130 {
131  height = width = ext.Height();
132 }
133 
134 
136  NonlinearFormExtension(form), fes(*form->FESpace())
137 {
140  if (elem_restrict_lex) // replace with a check for not identity
141  {
144  localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
145  }
146 }
147 
149 {
150  const Array<NonlinearFormIntegrator*> &integrators = *nlf->GetDNFI();
151  const int Ni = integrators.Size();
152  for (int i = 0; i < Ni; ++i)
153  {
154  integrators[i]->AssembleMF(fes);
155  }
156 }
157 
159 {
160  const Array<NonlinearFormIntegrator*> &integrators = *nlf->GetDNFI();
161  const int iSz = integrators.Size();
162  // replace the check 'elem_restrict_lex' with a check for not identity
163  if (elem_restrict_lex && !DeviceCanUseCeed())
164  {
166  localY = 0.0;
167  for (int i = 0; i < iSz; ++i)
168  {
169  integrators[i]->AddMultMF(localX, localY);
170  }
172  }
173  else
174  {
175  y.UseDevice(true); // typically this is a large vector, so store on device
176  y = 0.0;
177  for (int i = 0; i < iSz; ++i)
178  {
179  integrators[i]->AddMultMF(x, y);
180  }
181  }
182 }
183 
185 {
186  height = width = fes.GetVSize();
189  if (elem_restrict_lex) // replace with a check for not identity
190  {
193  }
194 }
195 
196 } // namespace mfem
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:563
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
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.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
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 Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1261
const NonlinearForm * nlf
Not owned.
void Update() override
Called by NonlinearForm::Update() to reflect changes in the FE space.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
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.
const Array< NonlinearFormIntegrator * > & GetBdrFaceIntegrators() const
Access all boundary face integrators added with AddBdrFaceIntegrator().
double GetGridFunctionEnergy(const Vector &x) const override
Compute the local (to the MPI rank) energy of the L-vector state x.
const Array< NonlinearFormIntegrator * > & GetInteriorFaceIntegrators() const
Access all interior face integrators added with AddInteriorFaceIntegrator().
void Update(Vector &x, int k, DenseMatrix &h, Vector &s, Array< Vector * > &v)
Definition: solvers.cpp:950
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:66
PANonlinearFormExtension(const NonlinearForm *nlf)
Lexicographic ordering for tensor-product FiniteElements.
Vector data type.
Definition: vector.hpp:60
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