MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
parheat.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
48#include "mtop_integrators.hpp"
49
50using namespace mfem;
51
52int main(int argc, char *argv[])
53{
54 // 1. Initialize MPI and HYPRE.
55 mfem::Mpi::Init(argc, argv);
56 int myrank = mfem::Mpi::WorldRank();
58
59 // Parse command-line options.
60 const char *mesh_file = "../../data/star.mesh";
61 int order = 1;
62 bool static_cond = false;
63 int ser_ref_levels = 1;
64 int par_ref_levels = 1;
65 real_t newton_rel_tol = 1e-7;
66 real_t newton_abs_tol = 1e-12;
67 int newton_iter = 10;
68 int print_level = 1;
69 bool visualization = false;
70
71 mfem::OptionsParser args(argc, argv);
72 args.AddOption(&mesh_file, "-m", "--mesh",
73 "Mesh file to use.");
74 args.AddOption(&ser_ref_levels,
75 "-rs",
76 "--refine-serial",
77 "Number of times to refine the mesh uniformly in serial.");
78 args.AddOption(&par_ref_levels,
79 "-rp",
80 "--refine-parallel",
81 "Number of times to refine the mesh uniformly in parallel.");
82 args.AddOption(&order, "-o", "--order",
83 "Finite element order (polynomial degree) or -1 for"
84 " isoparametric space.");
85 args.AddOption(&visualization,
86 "-vis",
87 "--visualization",
88 "-no-vis",
89 "--no-visualization",
90 "Enable or disable GLVis visualization.");
91 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
92 "--no-static-condensation", "Enable static condensation.");
93 args.AddOption(&newton_rel_tol,
94 "-rel",
95 "--relative-tolerance",
96 "Relative tolerance for the Newton solve.");
97 args.AddOption(&newton_abs_tol,
98 "-abs",
99 "--absolute-tolerance",
100 "Absolute tolerance for the Newton solve.");
101 args.AddOption(&newton_iter,
102 "-it",
103 "--newton-iterations",
104 "Maximum iterations for the Newton solve.");
105 args.Parse();
106 if (!args.Good())
107 {
108 if (myrank == 0)
109 {
110 args.PrintUsage(std::cout);
111 }
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.
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 real_t 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 real_t 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 real_t lsc=1.0;
305
306 real_t gQoI=ob->GetEnergy(solbv);
307 real_t lQoI;
308
309 real_t nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
310 real_t 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 real_t 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 return 0;
354}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
A class to handle Block systems in a matrix-free implementation.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
A class to handle Vectors in a block fashion.
void Update(real_t *data, const Array< int > &bOffsets)
Update method.
Vector & GetBlock(int i)
Get the i-th vector in the block.
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
GMRES method.
Definition solvers.hpp:547
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition hypre.cpp:5092
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
A class representing a general parallel block nonlinear operator defined on the Cartesian product of ...
virtual real_t GetEnergy(const Vector &x) const
Computes the energy of the system.
virtual void Mult(const Vector &x, Vector &y) const
Block T-Vector to Block T-Vector.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
Class for parallel meshes.
Definition pmesh.hpp:34
A class representing a general parametric parallel block nonlinear operator defined on the Cartesian ...
virtual BlockOperator & GetGradient(const Vector &x) const
Return the block gradient matrix for the given true-dof vector x.
virtual void SetParamFields(const Vector &dv) const
Set the parameters/design fields.
virtual void SetAdjointFields(const Vector &av) const
Set the adjoint fields.
virtual real_t GetEnergy(const Vector &x) const
Computes the energy of the system.
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!
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!...
virtual void SetStateFields(const Vector &xv) const
Set the state 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 ...
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
const Array< int > & GetBlockTrueOffsets() const
Return the true-dof offsets.
const Array< int > & ParamGetBlockTrueOffsets() const
Return the true-dof offsets for the parameters.
void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
Adds new Domain Integrator.
void Randomize(int seed=0)
Set random values in the vector.
Definition vector.cpp:816
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
int dim
Definition ex24.cpp:53
int main()
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:439
float real_t
Definition config.hpp:43