MFEM  v4.3.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-2021, 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  // Initialize MPI.
53  int nprocs, myrank;
54  MPI_Init(&argc, &argv);
55  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
56  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
57 
58  // Parse command-line options.
59  const char *mesh_file = "../../data/star.mesh";
60  int order = 1;
61  bool static_cond = false;
62  int ser_ref_levels = 1;
63  int par_ref_levels = 1;
64  double newton_rel_tol = 1e-7;
65  double newton_abs_tol = 1e-12;
66  int newton_iter = 10;
67  int print_level = 1;
68  bool visualization = false;
69 
70  mfem::OptionsParser args(argc, argv);
71  args.AddOption(&mesh_file, "-m", "--mesh",
72  "Mesh file to use.");
73  args.AddOption(&ser_ref_levels,
74  "-rs",
75  "--refine-serial",
76  "Number of times to refine the mesh uniformly in serial.");
77  args.AddOption(&par_ref_levels,
78  "-rp",
79  "--refine-parallel",
80  "Number of times to refine the mesh uniformly in parallel.");
81  args.AddOption(&order, "-o", "--order",
82  "Finite element order (polynomial degree) or -1 for"
83  " isoparametric space.");
84  args.AddOption(&visualization,
85  "-vis",
86  "--visualization",
87  "-no-vis",
88  "--no-visualization",
89  "Enable or disable GLVis visualization.");
90  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
91  "--no-static-condensation", "Enable static condensation.");
92  args.AddOption(&newton_rel_tol,
93  "-rel",
94  "--relative-tolerance",
95  "Relative tolerance for the Newton solve.");
96  args.AddOption(&newton_abs_tol,
97  "-abs",
98  "--absolute-tolerance",
99  "Absolute tolerance for the Newton solve.");
100  args.AddOption(&newton_iter,
101  "-it",
102  "--newton-iterations",
103  "Maximum iterations for the Newton solve.");
104  args.Parse();
105  if (!args.Good())
106  {
107  if (myrank == 0)
108  {
109  args.PrintUsage(std::cout);
110  }
111  MPI_Finalize();
112  return 1;
113  }
114 
115  if (myrank == 0)
116  {
117  args.PrintOptions(std::cout);
118  }
119 
120  // Read the (serial) mesh from the given mesh file on all processors. We
121  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
122  // and volume meshes with the same code.
123  mfem::Mesh mesh(mesh_file, 1, 1);
124  int dim = mesh.Dimension();
125 
126  // Refine the serial mesh on all processors to increase the resolution. In
127  // this example we do 'ref_levels' of uniform refinement. We choose
128  // 'ref_levels' to be the largest number that gives a final mesh with no
129  // more than 10,000 elements.
130  {
131  int ref_levels =
132  (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
133  for (int l = 0; l < ref_levels; l++)
134  {
135  mesh.UniformRefinement();
136  }
137  }
138 
139  // Define a parallel mesh by a partitioning of the serial mesh. Refine
140  // this mesh further in parallel to increase the resolution. Once the
141  // parallel mesh is defined, the serial mesh can be deleted.
142  mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
143  mesh.Clear();
144  {
145  for (int l = 0; l < par_ref_levels; l++)
146  {
147  pmesh.UniformRefinement();
148  }
149  }
150 
151  // Define the Diffusion coefficient.
153  // Define the Heat source.
155  // Define the q-function.
156  mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
157  1e-7,4.0,0.5);
158 
159  // Define FE collection and space for the state solution.
160  mfem::H1_FECollection sfec(order, dim);
162  1);
163  // Define FE collection and space for the density field.
164  mfem::L2_FECollection pfec(order, dim);
166  1);
167 
168  // Define the arrays for the nonlinear form.
171 
172  asfes.Append(sfes);
173  apfes.Append(pfes);
174 
175  // Define parametric block nonlinear form using single scalar H1 field
176  // and L2 scalar density field.
178  // Add a parametric integrator.
180 
181  // Define true block vectors for state, adjoint, resudual.
182  mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
183  mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
184  mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
185  // Define true block vectors for parametric field and gradients.
187  prmbv=0.0;
189  grdbv=0.0;
190 
191  // Set the BCs for the physics.
192  mfem::Array<mfem::Array<int> *> ess_bdr;
194  ess_bdr.Append(new mfem::Array<int>(pmesh.bdr_attributes.Max()));
195  ess_rhs.Append(nullptr);
196  (*ess_bdr[0]) = 1;
197  nf->SetEssentialBC(ess_bdr,ess_rhs);
198  delete ess_bdr[0];
199 
200  // Set the density field to 0.5.
201  prmbv=0.5;
202  // Set the density as parametric field in the parametric BNLForm.
203  nf->SetParamFields(prmbv); //set the density
204 
205  // Compute the stiffness/tangent matrix for density prmbv=0.5.
206  mfem::BlockOperator *A = &nf->GetGradient(solbv);
208  prec->SetPrintLevel(print_level);
209  // Use only block (0,0) as in this case we have a single field.
210  prec->SetOperator(A->GetBlock(0,0));
211 
212  // Construct block preconditioner for the BNLForm.
214  nf->GetBlockTrueOffsets());
215  blpr->SetDiagonalBlock(0,prec);
216 
217  // Define the solvers.
218  mfem::GMRESSolver *gmres;
219  gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
220  gmres->SetAbsTol(newton_abs_tol/10);
221  gmres->SetRelTol(newton_rel_tol/10);
222  gmres->SetMaxIter(100);
223  gmres->SetPrintLevel(print_level);
224  gmres->SetPreconditioner(*blpr);
225  gmres->SetOperator(*A);
226 
227 
228  // Solve the problem.
229  solbv=0.0;
230  nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
231  gmres->Mult(resbv, solbv);
232 
233  // Compute the energy of the state system.
234  double energy = nf->GetEnergy(solbv);
235  if (myrank==0)
236  {
237  std::cout << "energy =" << energy << std::endl;
238  }
239 
240  // Define the block nonlinear form utilized for representing the objective -
241  // use the state array from the BNLForm.
243  // Add the integrator for the objective.
245 
246  // Compute the objective.
247  double obj=ob->GetEnergy(solbv);
248  if (myrank==0)
249  {
250  std::cout << "Objective =" << obj << std::endl;
251  }
252 
253  // Solve the adjoint.
254  {
255  mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets()); adjrhs=0.0;
256  // Compute the RHS for the adjoint, i.e., the gradients with respect to
257  // the parametric fields.
258  ob->Mult(solbv, adjrhs);
259  // Get the tangent matrix from the state problem. We do not need to
260  // transpose the operator for diffusion. Compute the adjoint solution.
261  gmres->Mult(adjrhs, adjbv);
262  }
263 
264  // Compute gradients.
265  // First set the adjoint field.
266  nf->SetAdjointFields(adjbv);
267  // Set the state field.
268  nf->SetStateFields(solbv);
269  // Call the parametric Mult.
270  nf->ParamMult(prmbv, grdbv);
271 
272  // Dump out the data.
273  if (visualization)
274  {
276  &pmesh);
277  mfem::ParGridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
278  mfem::ParGridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
279  // Define state grid function.
280  mfem::ParGridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
281  mfem::ParGridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
282 
283  dacol->SetLevelsOfDetail(order);
284  dacol->RegisterField("sol", &gfsol);
285  dacol->RegisterField("adj", &gfadj);
286  dacol->RegisterField("dns", &gfdns);
287  dacol->RegisterField("grd", &gfgrd);
288 
289  dacol->SetTime(1.0);
290  dacol->SetCycle(1);
291  dacol->Save();
292 
293  delete dacol;
294  }
295 
296  // FD check
297  {
298  mfem::BlockVector prtbv;
299  mfem::BlockVector tmpbv;
300  prtbv.Update(nf->ParamGetBlockTrueOffsets());
301  tmpbv.Update(nf->ParamGetBlockTrueOffsets());
302  prtbv.GetBlock(0).Randomize();
303  prtbv*=1.0;
304  double lsc=1.0;
305 
306  double gQoI=ob->GetEnergy(solbv);
307  double lQoI;
308 
309  double nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
310  double td=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,grdbv);
311  td=td/nd;
312 
313  for (int l = 0; l < 10; l++)
314  {
315  lsc/=10.0;
316  prtbv/=10.0;
317  add(prmbv,prtbv,tmpbv);
318  nf->SetParamFields(tmpbv);
319  // Solve the physics.
320  solbv=0.0;
321  nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
322  A = &nf->GetGradient(solbv);
323  prec->SetPrintLevel(0);
324  prec->SetOperator(A->GetBlock(0,0));
325  gmres->SetOperator(*A);
326  gmres->SetPrintLevel(0);
327  gmres->Mult(resbv,solbv);
328  // Compute the objective.
329  lQoI=ob->GetEnergy(solbv);
330  double ld=(lQoI-gQoI)/lsc;
331  if (myrank==0)
332  {
333  std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
334  << " adjoint gradient=" << td
335  << " err=" << std::fabs(ld/nd-td) << std::endl;
336  }
337  }
338  }
339 
340  delete ob;
341  delete gmres;
342  delete blpr;
343  delete prec;
344 
345  delete nf;
346  delete pfes;
347  delete sfes;
348 
349  delete qfun;
350  delete loadco;
351  delete diffco;
352 
353  MPI_Finalize();
354  return 0;
355 }
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 ...
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:846
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:763
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:291
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 Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:884
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:98
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:164
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
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:746
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1467
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
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:99
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
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:348
A class representing a general parallel block nonlinear operator defined on the Cartesian product of ...
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:316
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: hypre.cpp:4447
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:327
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
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:285
void Neg()
(*this) = -(*this)
Definition: vector.cpp:283
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150