MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
parheat.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 // ----------------------------------------------------------------
13 // ParHeat Miniapp: Gradients of PDE constrained objective function
14 // ----------------------------------------------------------------
15 // (Parallel Version)
16 //
17 // The following example computes the gradients of a specified objective
18 // function with respect to parametric fields. The objective function is having
19 // the following form f(u(\rho)) where u(\rho) is a solution of a specific state
20 // problem (in the example that is the diffusion equation), and \rho is a
21 // parametric field discretized by finite elements. The parametric field (also
22 // called density in topology optimization) controls the coefficients of the
23 // state equation. For the considered case, the density controls the diffusion
24 // coefficient within the computational domain.
25 //
26 // For more information, the users are referred to:
27 //
28 // Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
29 // Optimization with PDE Constraints
30 // Springer Netherlands, 2009
31 //
32 // Bendsøe, M. P. & Sigmund, O.
33 // Topology Optimization - Theory, Methods and Applications
34 // Springer Verlag, Berlin Heidelberg, 2003
35 //
36 // Compile with: make parheat
37 //
38 // Sample runs:
39 //
40 // mpirun -np 4 parheat --visualization
41 // mpirun -np 4 parheat --visualization -m ../../data/beam-quad.mesh
42 
43 #include "mfem.hpp"
44 #include <fstream>
45 #include <iostream>
46 
47 #include "pparamnonlinearform.hpp"
48 #include "mtop_integrators.hpp"
49 
50 int main(int argc, char *argv[])
51 {
52  // 1. Initialize MPI and HYPRE.
53  mfem::Mpi::Init(argc, argv);
54  int myrank = mfem::Mpi::WorldRank();
56 
57  // Parse command-line options.
58  const char *mesh_file = "../../data/star.mesh";
59  int order = 1;
60  bool static_cond = false;
61  int ser_ref_levels = 1;
62  int par_ref_levels = 1;
63  double newton_rel_tol = 1e-7;
64  double newton_abs_tol = 1e-12;
65  int newton_iter = 10;
66  int print_level = 1;
67  bool visualization = false;
68 
69  mfem::OptionsParser args(argc, argv);
70  args.AddOption(&mesh_file, "-m", "--mesh",
71  "Mesh file to use.");
72  args.AddOption(&ser_ref_levels,
73  "-rs",
74  "--refine-serial",
75  "Number of times to refine the mesh uniformly in serial.");
76  args.AddOption(&par_ref_levels,
77  "-rp",
78  "--refine-parallel",
79  "Number of times to refine the mesh uniformly in parallel.");
80  args.AddOption(&order, "-o", "--order",
81  "Finite element order (polynomial degree) or -1 for"
82  " isoparametric space.");
83  args.AddOption(&visualization,
84  "-vis",
85  "--visualization",
86  "-no-vis",
87  "--no-visualization",
88  "Enable or disable GLVis visualization.");
89  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
90  "--no-static-condensation", "Enable static condensation.");
91  args.AddOption(&newton_rel_tol,
92  "-rel",
93  "--relative-tolerance",
94  "Relative tolerance for the Newton solve.");
95  args.AddOption(&newton_abs_tol,
96  "-abs",
97  "--absolute-tolerance",
98  "Absolute tolerance for the Newton solve.");
99  args.AddOption(&newton_iter,
100  "-it",
101  "--newton-iterations",
102  "Maximum iterations for the Newton solve.");
103  args.Parse();
104  if (!args.Good())
105  {
106  if (myrank == 0)
107  {
108  args.PrintUsage(std::cout);
109  }
110  return 1;
111  }
112 
113  if (myrank == 0)
114  {
115  args.PrintOptions(std::cout);
116  }
117 
118  // Read the (serial) mesh from the given mesh file on all processors. We
119  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
120  // and volume meshes with the same code.
121  mfem::Mesh mesh(mesh_file, 1, 1);
122  int dim = mesh.Dimension();
123 
124  // Refine the serial mesh on all processors to increase the resolution. In
125  // this example we do 'ref_levels' of uniform refinement. We choose
126  // 'ref_levels' to be the largest number that gives a final mesh with no
127  // more than 10,000 elements.
128  {
129  int ref_levels =
130  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
131  for (int l = 0; l < ref_levels; l++)
132  {
133  mesh.UniformRefinement();
134  }
135  }
136 
137  // Define a parallel mesh by a partitioning of the serial mesh. Refine
138  // this mesh further in parallel to increase the resolution. Once the
139  // parallel mesh is defined, the serial mesh can be deleted.
140  mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
141  mesh.Clear();
142  {
143  for (int l = 0; l < par_ref_levels; l++)
144  {
145  pmesh.UniformRefinement();
146  }
147  }
148 
149  // Define the Diffusion coefficient.
151  // Define the Heat source.
153  // Define the q-function.
154  mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
155  1e-7,4.0,0.5);
156 
157  // Define FE collection and space for the state solution.
158  mfem::H1_FECollection sfec(order, dim);
160  1);
161  // Define FE collection and space for the density field.
162  mfem::L2_FECollection pfec(order, dim);
164  1);
165 
166  // Define the arrays for the nonlinear form.
169 
170  asfes.Append(sfes);
171  apfes.Append(pfes);
172 
173  // Define parametric block nonlinear form using single scalar H1 field
174  // and L2 scalar density field.
176  // Add a parametric integrator.
178 
179  // Define true block vectors for state, adjoint, resudual.
180  mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
181  mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
182  mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
183  // Define true block vectors for parametric field and gradients.
185  prmbv=0.0;
187  grdbv=0.0;
188 
189  // Set the BCs for the physics.
190  mfem::Array<mfem::Array<int> *> ess_bdr;
192  ess_bdr.Append(new mfem::Array<int>(pmesh.bdr_attributes.Max()));
193  ess_rhs.Append(nullptr);
194  (*ess_bdr[0]) = 1;
195  nf->SetEssentialBC(ess_bdr,ess_rhs);
196  delete ess_bdr[0];
197 
198  // Set the density field to 0.5.
199  prmbv=0.5;
200  // Set the density as parametric field in the parametric BNLForm.
201  nf->SetParamFields(prmbv); //set the density
202 
203  // Compute the stiffness/tangent matrix for density prmbv=0.5.
204  mfem::BlockOperator *A = &nf->GetGradient(solbv);
206  prec->SetPrintLevel(print_level);
207  // Use only block (0,0) as in this case we have a single field.
208  prec->SetOperator(A->GetBlock(0,0));
209 
210  // Construct block preconditioner for the BNLForm.
212  nf->GetBlockTrueOffsets());
213  blpr->SetDiagonalBlock(0,prec);
214 
215  // Define the solvers.
216  mfem::GMRESSolver *gmres;
217  gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
218  gmres->SetAbsTol(newton_abs_tol/10);
219  gmres->SetRelTol(newton_rel_tol/10);
220  gmres->SetMaxIter(100);
221  gmres->SetPrintLevel(print_level);
222  gmres->SetPreconditioner(*blpr);
223  gmres->SetOperator(*A);
224 
225 
226  // Solve the problem.
227  solbv=0.0;
228  nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
229  gmres->Mult(resbv, solbv);
230 
231  // Compute the energy of the state system.
232  double energy = nf->GetEnergy(solbv);
233  if (myrank==0)
234  {
235  std::cout << "energy =" << energy << std::endl;
236  }
237 
238  // Define the block nonlinear form utilized for representing the objective -
239  // use the state array from the BNLForm.
241  // Add the integrator for the objective.
243 
244  // Compute the objective.
245  double obj=ob->GetEnergy(solbv);
246  if (myrank==0)
247  {
248  std::cout << "Objective =" << obj << std::endl;
249  }
250 
251  // Solve the adjoint.
252  {
253  mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets()); adjrhs=0.0;
254  // Compute the RHS for the adjoint, i.e., the gradients with respect to
255  // the parametric fields.
256  ob->Mult(solbv, adjrhs);
257  // Get the tangent matrix from the state problem. We do not need to
258  // transpose the operator for diffusion. Compute the adjoint solution.
259  gmres->Mult(adjrhs, adjbv);
260  }
261 
262  // Compute gradients.
263  // First set the adjoint field.
264  nf->SetAdjointFields(adjbv);
265  // Set the state field.
266  nf->SetStateFields(solbv);
267  // Call the parametric Mult.
268  nf->ParamMult(prmbv, grdbv);
269 
270  // Dump out the data.
271  if (visualization)
272  {
274  &pmesh);
275  mfem::ParGridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
276  mfem::ParGridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
277  // Define state grid function.
278  mfem::ParGridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
279  mfem::ParGridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
280 
281  dacol->SetLevelsOfDetail(order);
282  dacol->RegisterField("sol", &gfsol);
283  dacol->RegisterField("adj", &gfadj);
284  dacol->RegisterField("dns", &gfdns);
285  dacol->RegisterField("grd", &gfgrd);
286 
287  dacol->SetTime(1.0);
288  dacol->SetCycle(1);
289  dacol->Save();
290 
291  delete dacol;
292  }
293 
294  // FD check
295  {
296  mfem::BlockVector prtbv;
297  mfem::BlockVector tmpbv;
298  prtbv.Update(nf->ParamGetBlockTrueOffsets());
299  tmpbv.Update(nf->ParamGetBlockTrueOffsets());
300  prtbv.GetBlock(0).Randomize();
301  prtbv*=1.0;
302  double lsc=1.0;
303 
304  double gQoI=ob->GetEnergy(solbv);
305  double lQoI;
306 
307  double nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
308  double td=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,grdbv);
309  td=td/nd;
310 
311  for (int l = 0; l < 10; l++)
312  {
313  lsc/=10.0;
314  prtbv/=10.0;
315  add(prmbv,prtbv,tmpbv);
316  nf->SetParamFields(tmpbv);
317  // Solve the physics.
318  solbv=0.0;
319  nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
320  A = &nf->GetGradient(solbv);
321  prec->SetPrintLevel(0);
322  prec->SetOperator(A->GetBlock(0,0));
323  gmres->SetOperator(*A);
324  gmres->SetPrintLevel(0);
325  gmres->Mult(resbv,solbv);
326  // Compute the objective.
327  lQoI=ob->GetEnergy(solbv);
328  double ld=(lQoI-gQoI)/lsc;
329  if (myrank==0)
330  {
331  std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
332  << " adjoint gradient=" << td
333  << " err=" << std::fabs(ld/nd-td) << std::endl;
334  }
335  }
336  }
337 
338  delete ob;
339  delete gmres;
340  delete blpr;
341  delete prec;
342 
343  delete nf;
344  delete pfes;
345  delete sfes;
346 
347  delete qfun;
348  delete loadco;
349  delete diffco;
350 
351  return 0;
352 }
virtual void SetParamFields(const Vector &dv) const
Set the parameters/design fields.
virtual void ParamMult(const Vector &x, Vector &y) const
Calculates the product of the adjoint field and the derivative of the state residual with respect to ...
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:67
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
Set the state essential BCs. Here, rhs is a true dof vector!
A class representing a general parametric parallel block nonlinear operator defined on the Cartesian ...
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
Helper class for ParaView visualization data.
virtual void Mult(const Vector &x, Vector &y) const
Calculates the residual for a state input given by block T-Vector. The result is Block T-Vector! The ...
void Update(double *data, const Array< int > &bOffsets)
Update method.
Definition: blockvector.cpp:85
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:772
virtual void Mult(const Vector &x, Vector &y) const
Block T-Vector to Block T-Vector.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
const Array< int > & GetBlockTrueOffsets() const
Return the true-dof offsets.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:971
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1443
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1523
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
static void Init()
Singleton creation with Mpi::Init();.
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double GetEnergy(const Vector &x) const
Computes the energy of the system.
virtual void SetStateFields(const Vector &xv) const
Set the state fields.
const Array< int > & ParamGetBlockTrueOffsets() const
Return the true-dof offsets for the parameters.
void SetAbsTol(double atol)
Definition: solvers.hpp:199
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void SetRelTol(double rtol)
Definition: solvers.hpp:198
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
virtual void SetAdjointFields(const Vector &av) const
Set the adjoint fields.
GMRES method.
Definition: solvers.hpp:497
A class representing a general parallel block nonlinear operator defined on the Cartesian product of ...
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:426
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4653
virtual double GetEnergy(const Vector &x) const
Computes the energy of the system.
int dim
Definition: ex24.cpp:53
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
virtual void Save() override
Class for parallel grid function.
Definition: pgridfunc.hpp:32
A class to handle Block systems in a matrix-free implementation.
Class for parallel meshes.
Definition: pmesh.hpp:32
void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
Adds new Domain Integrator.
int main()
virtual BlockOperator & GetGradient(const Vector &x) const
Return the block gradient matrix for the given true-dof vector x.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
void Neg()
(*this) = -(*this)
Definition: vector.cpp:292
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150