MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pnonlinearform.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 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "fem.hpp"
17 #include "../general/forall.hpp"
18 
19 namespace mfem
20 {
21 
23  : NonlinearForm(pf), pGrad(Operator::Hypre_ParCSR)
24 {
25  X.MakeRef(pf, NULL);
26  Y.MakeRef(pf, NULL);
27  MFEM_VERIFY(!Serial(), "internal MFEM error");
28 }
29 
31 {
32  double loc_energy, glob_energy;
33 
34  loc_energy = GetGridFunctionEnergy(x);
35 
36  if (fnfi.Size())
37  {
38  MFEM_ABORT("TODO: add energy contribution from shared faces");
39  }
40 
41  MPI_Allreduce(&loc_energy, &glob_energy, 1, MPI_DOUBLE, MPI_SUM,
42  ParFESpace()->GetComm());
43 
44  return glob_energy;
45 }
46 
47 void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
48 {
49  NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
50 
51  if (fnfi.Size())
52  {
53  MFEM_VERIFY(!NonlinearForm::ext, "Not implemented (extensions + faces");
54  // Terms over shared interior faces in parallel.
56  ParMesh *pmesh = pfes->GetParMesh();
58  const FiniteElement *fe1, *fe2;
59  Array<int> vdofs1, vdofs2;
60  Vector el_x, el_y;
61 
63  X.MakeRef(aux1, 0); // aux1 contains P.x
65  const int n_shared_faces = pmesh->GetNSharedFaces();
66  for (int i = 0; i < n_shared_faces; i++)
67  {
68  tr = pmesh->GetSharedFaceTransformations(i, true);
69  int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
70 
71  fe1 = pfes->GetFE(tr->Elem1No);
72  fe2 = pfes->GetFaceNbrFE(Elem2NbrNo);
73 
74  pfes->GetElementVDofs(tr->Elem1No, vdofs1);
75  pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
76 
77  el_x.SetSize(vdofs1.Size() + vdofs2.Size());
78  X.GetSubVector(vdofs1, el_x.GetData());
79  X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
80 
81  for (int k = 0; k < fnfi.Size(); k++)
82  {
83  fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
84  aux2.AddElementVector(vdofs1, el_y.GetData());
85  }
86  }
87  }
88 
89  P->MultTranspose(aux2, y);
90 
91  const int N = ess_tdof_list.Size();
92  const auto idx = ess_tdof_list.Read();
93  auto Y_RW = y.ReadWrite();
94  MFEM_FORALL(i, N, Y_RW[idx[i]] = 0.0; );
95 }
96 
98 {
99  MFEM_VERIFY(NonlinearForm::ext == nullptr,
100  "this method is not supported yet with partial assembly");
101 
102  NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
103 
104  return *Grad;
105 }
106 
108 {
110 
112 
113  pGrad.Clear();
114 
115  NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
116 
117  OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
118 
119  if (fnfi.Size() == 0)
120  {
121  dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
122  pfes->GetDofOffsets(), Grad);
123  }
124  else
125  {
126  MFEM_ABORT("TODO: assemble contributions from shared face terms");
127  }
128 
129  // RAP the local gradient dA.
130  // TODO - construct Dof_TrueDof_Matrix directly in the pGrad format
131  Ph.ConvertFrom(pfes->Dof_TrueDof_Matrix());
132  pGrad.MakePtAP(dA, Ph);
133 
134  // Impose b.c. on pGrad
135  OperatorHandle pGrad_e;
137 
138  return *pGrad.Ptr();
139 }
140 
142 {
143  Y.MakeRef(ParFESpace(), NULL);
144  X.MakeRef(ParFESpace(), NULL);
145  pGrad.Clear();
147 }
148 
149 
152 {
153  pBlockGrad = NULL;
154  SetParSpaces(pf);
155 }
156 
158 {
159  delete pBlockGrad;
160  pBlockGrad = NULL;
161 
162  for (int s1=0; s1<fes.Size(); ++s1)
163  {
164  for (int s2=0; s2<fes.Size(); ++s2)
165  {
166  delete phBlockGrad(s1,s2);
167  }
168  }
169 
170  Array<FiniteElementSpace *> serialSpaces(pf.Size());
171 
172  for (int s=0; s<pf.Size(); s++)
173  {
174  serialSpaces[s] = (FiniteElementSpace *) pf[s];
175  }
176 
177  SetSpaces(serialSpaces);
178 
179  phBlockGrad.SetSize(fes.Size(), fes.Size());
180 
181  for (int s1=0; s1<fes.Size(); ++s1)
182  {
183  for (int s2=0; s2<fes.Size(); ++s2)
184  {
186  }
187  }
188 }
189 
191 {
192  return (ParFiniteElementSpace *)fes[k];
193 }
194 
196 {
197  return (const ParFiniteElementSpace *)fes[k];
198 }
199 
200 // Here, rhs is a true dof vector
202  Array<Array<int> *>&bdr_attr_is_ess,
203  Array<Vector *> &rhs)
204 {
205  Array<Vector *> nullarray(fes.Size());
206  nullarray = NULL;
207 
208  BlockNonlinearForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
209 
210  for (int s = 0; s < fes.Size(); ++s)
211  {
212  if (rhs[s])
213  {
214  rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
215  }
216  }
217 }
218 
220 {
221  // xs_true is not modified, so const_cast is okay
222  xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
224 
225  for (int s = 0; s < fes.Size(); ++s)
226  {
227  fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
228  }
229 
231  double englo = 0.0;
232 
233  MPI_Allreduce(&enloc, &englo, 1, MPI_DOUBLE, MPI_SUM,
234  ParFESpace(0)->GetComm());
235 
236  return englo;
237 }
238 
239 void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
240 {
241  // xs_true is not modified, so const_cast is okay
242  xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
246 
247  for (int s=0; s<fes.Size(); ++s)
248  {
249  fes[s]->GetProlongationMatrix()->Mult(
251  }
252 
254 
255  if (fnfi.Size() > 0)
256  {
257  MFEM_ABORT("TODO: assemble contributions from shared face terms");
258  }
259 
260  for (int s=0; s<fes.Size(); ++s)
261  {
262  fes[s]->GetProlongationMatrix()->MultTranspose(
264 
266  }
267 
269  y.SyncMemory(ys_true);
270 }
271 
272 /// Return the local gradient matrix for the given true-dof vector x
274  const Vector &x) const
275 {
276  // xs_true is not modified, so const_cast is okay
277  xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
279 
280  for (int s=0; s<fes.Size(); ++s)
281  {
282  fes[s]->GetProlongationMatrix()->Mult(
284  }
285 
286  // (re)assemble Grad without b.c. into 'Grads'
288 
289  delete BlockGrad;
291 
292  for (int i = 0; i < fes.Size(); ++i)
293  {
294  for (int j = 0; j < fes.Size(); ++j)
295  {
296  BlockGrad->SetBlock(i, j, Grads(i, j));
297  }
298  }
299  return *BlockGrad;
300 }
301 
302 // Set the operator type id for the parallel gradient matrix/operator.
304 {
305  for (int s1=0; s1<fes.Size(); ++s1)
306  {
307  for (int s2=0; s2<fes.Size(); ++s2)
308  {
309  phBlockGrad(s1,s2)->SetType(tid);
310  }
311  }
312 }
313 
315 {
316  if (pBlockGrad == NULL)
317  {
319  }
320 
322 
323  for (int s1=0; s1<fes.Size(); ++s1)
324  {
325  pfes[s1] = ParFESpace(s1);
326 
327  for (int s2=0; s2<fes.Size(); ++s2)
328  {
329  phBlockGrad(s1,s2)->Clear();
330  }
331  }
332 
333  GetLocalGradient(x); // gradients are stored in 'Grads'
334 
335  if (fnfi.Size() > 0)
336  {
337  MFEM_ABORT("TODO: assemble contributions from shared face terms");
338  }
339 
340  for (int s1=0; s1<fes.Size(); ++s1)
341  {
342  for (int s2=0; s2<fes.Size(); ++s2)
343  {
344  OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
345  Ph(phBlockGrad(s1,s2)->Type()),
346  Rh(phBlockGrad(s1,s2)->Type());
347 
348  if (s1 == s2)
349  {
350  dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
351  pfes[s1]->GetDofOffsets(), Grads(s1,s1));
352  Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
353  phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
354 
355  OperatorHandle Ae;
356  Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
357  }
358  else
359  {
360  dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
361  pfes[s1]->GlobalVSize(),
362  pfes[s2]->GlobalVSize(),
363  pfes[s1]->GetDofOffsets(),
364  pfes[s2]->GetDofOffsets(),
365  Grads(s1,s2));
366  Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
367  Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
368 
369  phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
370 
371  phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
372  phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
373  }
374 
375  pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
376  }
377  }
378 
379  return *pBlockGrad;
380 }
381 
383 {
384  delete pBlockGrad;
385  for (int s1=0; s1<fes.Size(); ++s1)
386  {
387  for (int s2=0; s2<fes.Size(); ++s2)
388  {
389  delete phBlockGrad(s1,s2);
390  }
391  }
392 }
393 
394 }
395 
396 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:235
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:200
ParFiniteElementSpace * ParFESpace() const
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3127
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:85
A class representing a general block nonlinear operator defined on the Cartesian product of multiple ...
virtual void Update()
Update the ParNonlinearForm to propagate updates of the associated parallel FE space.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
Abstract parallel finite element space.
Definition: pfespace.hpp:28
Array< int > ess_tdof_list
A list of all essential true dofs.
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
virtual void Mult(const Vector &x, Vector &y) const
Block T-Vector to Block T-Vector.
Array2D< SparseMatrix * > Grads
ParBlockNonlinearForm()
Construct an empty ParBlockNonlinearForm. Initialize with SetParSpaces().
void SetSpaces(Array< FiniteElementSpace * > &f)
(Re)initialize the BlockNonlinearForm.
void SetGradientType(Operator::Type tid)
Set the operator type id for the blocks of the parallel gradient matrix/operator. The default type is...
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:94
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:102
ParNonlinearForm(ParFiniteElementSpace *pf)
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:283
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1517
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
Array< FiniteElementSpace * > fes
FE spaces on which the form lives.
Data type sparse matrix.
Definition: sparsemat.hpp:50
Array< NonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
void EliminateRowsCols(OperatorHandle &A, const Array< int > &ess_dof_list)
Reset the OperatorHandle to be the eliminated part of A after elimination of the essential dofs ess_d...
Definition: handle.cpp:252
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3011
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
Vector aux1
Auxiliary Vectors.
Array< BlockNonlinearFormIntegrator * > fnfi
Set of interior face Integrators to be assembled (added).
const Operator * P
Pointer to the prolongation matrix of fes, may be NULL.
const BlockOperator & GetLocalGradient(const Vector &x) const
Return the local block gradient matrix for the given true-dof vector x.
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:281
void SyncMemory(const Vector &v) const
Update the memory location of the vector to match v.
Definition: vector.hpp:242
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:639
void SetParSpaces(Array< ParFiniteElementSpace * > &pf)
After a call to SetParSpaces(), the essential b.c. and the gradient-type (if different from the defau...
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
virtual double GetEnergy(const Vector &x) const
Computes the energy of the system.
double GetParGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Vector & FaceNbrData()
Definition: pgridfunc.hpp:208
bool Serial() const
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:256
virtual void Update()
Update the NonlinearForm to propagate updates of the associated FE space.
void MakeSquareBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt *row_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel square block-diagonal matrix using the currently set type...
Definition: handle.cpp:60
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
const SparseMatrix & GetLocalGradient(const Vector &x) const
Return the local gradient matrix for the given true-dof vector x.
NonlinearFormExtension * ext
Extension for supporting different AssemblyLevels.
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:124
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:465
double GetGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Array2D< OperatorHandle * > phBlockGrad
virtual Operator & GetGradient(const Vector &x) const
Compute the gradient Operator of the NonlinearForm corresponding to the state x.
double GetEnergyBlocked(const BlockVector &bx) const
Specialized version of GetEnergy() for BlockVectors.
void SyncFromBlocks() const
Synchronize the memory location flags (i.e. the memory validity flags) of the big/monolithic block-ve...
Vector data type.
Definition: vector.hpp:60
void ComputeGradientBlocked(const BlockVector &bx) const
Specialized version of GetGradient() for BlockVector.
SparseMatrix * Grad
ID for class HypreParMatrix.
Definition: operator.hpp:260
virtual BlockOperator & GetGradient(const Vector &x) const
virtual ~ParBlockNonlinearForm()
Destructor.
RefCoord s[3]
Array< Array< int > * > ess_tdofs
ParFiniteElementSpace * ParFESpace(int k)
Return the k-th parallel FE space of the ParBlockNonlinearForm.
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:32
void MultBlocked(const BlockVector &bx, BlockVector &by) const
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:469
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:123
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1464
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)