MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
seqheat.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 // SeqHeat Miniapp: Gradients of PDE constrained objective function
14 // ----------------------------------------------------------------
15 // (Sequential 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 seqheat
37 //
38 // Sample runs:
39 //
40 // seqheat -m ../../data/star-mixed.mesh
41 // seqheat --visualization
42 
43 #include "mfem.hpp"
44 #include <fstream>
45 #include <iostream>
46 
47 #include "mtop_integrators.hpp"
48 
49 int main(int argc, char *argv[])
50 {
51  const char *mesh_file = "../../data/star.vtk";
52  int ser_ref_levels = 1;
53  int order = 2;
54  bool visualization = false;
55  double newton_rel_tol = 1e-4;
56  double newton_abs_tol = 1e-6;
57  int newton_iter = 10;
58  int print_level = 0;
59 
60  mfem::OptionsParser args(argc, argv);
61  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
62  args.AddOption(&ser_ref_levels,
63  "-rs",
64  "--refine-serial",
65  "Number of times to refine the mesh uniformly in serial.");
66  args.AddOption(&order,
67  "-o",
68  "--order",
69  "Order (degree) of the finite elements.");
70  args.AddOption(&visualization,
71  "-vis",
72  "--visualization",
73  "-no-vis",
74  "--no-visualization",
75  "Enable or disable GLVis visualization.");
76  args.AddOption(&newton_rel_tol,
77  "-rel",
78  "--relative-tolerance",
79  "Relative tolerance for the Newton solve.");
80  args.AddOption(&newton_abs_tol,
81  "-abs",
82  "--absolute-tolerance",
83  "Absolute tolerance for the Newton solve.");
84  args.AddOption(&newton_iter,
85  "-it",
86  "--newton-iterations",
87  "Maximum iterations for the Newton solve.");
88  args.Parse();
89  if (!args.Good())
90  {
91  args.PrintUsage(std::cout);
92  return 1;
93  }
94  args.PrintOptions(std::cout);
95 
96  // Read the (serial) mesh from the given mesh file on all processors. We
97  // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
98  // with the same code.
99  mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
100  int dim = mesh->Dimension();
101 
102  // Refine the mesh in serial to increase the resolution. In this example
103  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
104  // a command-line parameter.
105  for (int lev = 0; lev < ser_ref_levels; lev++)
106  {
107  mesh->UniformRefinement();
108  }
109 
110  // Diffusion coefficient
112  // Heat source
114  // Define the q-function
115  mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
116  1e-7,4.0,0.5);
117 
118  // Define FE collection and space for the state solution
119  mfem::H1_FECollection sfec(order, dim);
120  mfem::FiniteElementSpace* sfes=new mfem::FiniteElementSpace(mesh,&sfec,1);
121  // Define FE collection and space for the density field
122  mfem::L2_FECollection pfec(order, dim);
123  mfem::FiniteElementSpace* pfes=new mfem::FiniteElementSpace(mesh,&pfec,1);
124 
125  // Define the arrays for the nonlinear form
128 
129  asfes.Append(sfes);
130  apfes.Append(pfes);
131  // Define parametric block nonlinear form using single scalar H1 field
132  // and L2 scalar density field
134  // Add the parametric integrator
136 
137  // Define true block vectors for state, adjoint, residual
138  mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
139  mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
140  mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
141  // Define true block vectors for parametric field and gradients
143  prmbv=0.0;
145  grdbv=0.0;
146 
147  // Set the BC for the physics
148  mfem::Array<mfem::Array<int> *> ess_bdr;
150  ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
151  ess_rhs.Append(nullptr);
152  (*ess_bdr[0]) = 1;
153  nf->SetEssentialBC(ess_bdr,ess_rhs);
154  delete ess_bdr[0];
155 
156  // Define the linear solvers
157  mfem::GMRESSolver *gmres;
158  gmres = new mfem::GMRESSolver();
159  gmres->SetAbsTol(newton_abs_tol/10);
160  gmres->SetRelTol(newton_rel_tol/10);
161  gmres->SetMaxIter(300);
162  gmres->SetPrintLevel(print_level);
163 
164  // Define the Newton solver
165  mfem::NewtonSolver *ns;
166  ns = new mfem::NewtonSolver();
167  ns->iterative_mode = true;
168  ns->SetSolver(*gmres);
169  ns->SetOperator(*nf);
170  ns->SetPrintLevel(print_level);
171  ns->SetRelTol(newton_rel_tol);
172  ns->SetAbsTol(newton_abs_tol);
173  ns->SetMaxIter(newton_iter);
174 
175  // Solve the problem
176  // Set the density to 0.5
177  prmbv=0.5;
178  nf->SetParamFields(prmbv); // Set the density
179  // Define the RHS
180  mfem::Vector b;
181  solbv=0.0;
182  // Newton solve
183  ns->Mult(b, solbv);
184 
185  // Compute the residual
186  nf->Mult(solbv,resbv);
187  std::cout<<"Norm residual="<<resbv.Norml2()<<std::endl;
188 
189  // Compute the energy of the state system
190  double energy = nf->GetEnergy(solbv);
191  std::cout<<"energy ="<< energy<<std::endl;
192 
193  // Define the block nonlinear form utilized for representing the
194  // objective. The input is the state array asfes defined earlier.
196 
197  // Add the integrator for the objective
199 
200  // Compute the objective
201  double obj=ob->GetEnergy(solbv);
202  std::cout<<"Objective ="<<obj<<std::endl;
203 
204  // Solve the adjoint
205  {
206  mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets()); adjrhs=0.0;
207  // Compute the RHS for the adjoint
208  ob->Mult(solbv, adjrhs);
209  // Get the tangent matrix from the state problem
210  mfem::BlockOperator& A=nf->GetGradient(solbv);
211  // We do not need to transpose the operator for diffusion
212  gmres->SetOperator(A.GetBlock(0,0));
213  // Compute the adjoint solution
214  gmres->Mult(adjrhs.GetBlock(0), adjbv.GetBlock(0));
215  }
216 
217  // Compute gradients
218  nf->SetAdjointFields(adjbv);
219  nf->SetStateFields(solbv);
220  nf->ParamMult(prmbv, grdbv);
221 
222  // Dump out the data
223  if (visualization)
224  {
226  mesh);
227  mfem::GridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
228  mfem::GridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
229  // Define state grid function
230  mfem::GridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
231  mfem::GridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
232 
233  dacol->SetLevelsOfDetail(order);
234  dacol->RegisterField("sol", &gfsol);
235  dacol->RegisterField("adj", &gfadj);
236  dacol->RegisterField("dns", &gfdns);
237  dacol->RegisterField("grd", &gfgrd);
238 
239  dacol->SetTime(1.0);
240  dacol->SetCycle(1);
241  dacol->Save();
242 
243  delete dacol;
244  }
245 
246  // FD check
247  {
248  // Perturbation vector
249  mfem::BlockVector prtbv;
250  mfem::BlockVector tmpbv;
251  prtbv.Update(nf->ParamGetBlockTrueOffsets());
252  tmpbv.Update(nf->ParamGetBlockTrueOffsets());
253  // Generate the perturbation
254  prtbv.GetBlock(0).Randomize();
255  prtbv*=1.0;
256  // Scaling parameter
257  double lsc=1.0;
258 
259  // Compute initial objective
260  double gQoI=ob->GetEnergy(solbv);
261  double lQoI;
262 
263  // Norm of the perturbation
264  double nd=mfem::InnerProduct(prtbv,prtbv);
265  // Projection of the adjoint gradient on the perturbation
266  double td=mfem::InnerProduct(prtbv,grdbv);
267  // Normalize the directional derivative
268  td=td/nd;
269 
270  for (int l = 0; l < 10; l++)
271  {
272  lsc/=10.0;
273  // Scale the perturbation
274  prtbv/=10.0;
275  // Add the perturbation to the original density
276  add(prmbv,prtbv,tmpbv);
277  nf->SetParamFields(tmpbv);
278  // Solve the physics
279  ns->Mult(b,solbv);
280  // Compute the objective
281  lQoI=ob->GetEnergy(solbv);
282  // FD approximation
283  double ld=(lQoI-gQoI)/lsc;
284  std::cout << "dx=" << lsc << " FD gradient=" << ld/nd
285  << " adjoint gradient=" << td
286  << " err=" << std::fabs(ld/nd-td) << std::endl;
287  }
288  }
289 
290  delete ob;
291 
292  delete ns;
293  delete gmres;
294 
295  delete nf;
296  delete pfes;
297  delete sfes;
298 
299  delete qfun;
300  delete loadco;
301  delete diffco;
302 
303  delete mesh;
304 
305  return 0;
306 }
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:783
A class representing a general parametric block nonlinear operator defined on the Cartesian product o...
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 SetAdjointFields(const Vector &av) const
Set the adjoint fields.
virtual void SetSolver(Solver &solver)
Set the linear solver for inverting the Jacobian.
Definition: solvers.hpp:512
virtual double GetEnergy(const Vector &x) const
Computes the energy for a state vector x.
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:763
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 double GetEnergy(const Vector &x) const
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:884
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
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:746
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
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
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:464
virtual void Mult(const Vector &x, Vector &y) const
virtual void SetEssentialBC(const Array< Array< int > * > &bdr_attr_is_ess, Array< Vector * > &rhs)
Set the essential boundary conditions.
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 void SetParamFields(const Vector &dv) const
Set the parameters/design fields.
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1653
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
virtual void SetStateFields(const Vector &xv) const
Set the state fields.
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
GMRES method.
Definition: solvers.hpp:348
virtual void ParamMult(const Vector &x, Vector &y) const
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:316
virtual void Mult(const Vector &x, Vector &y) const
virtual BlockOperator & GetGradient(const Vector &x) const
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_)
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
virtual void Save() override
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:347
A class to handle Block systems in a matrix-free implementation.
void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
Adds new Domain Integrator.
int main()
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1642
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
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150