libMesh
optimization_ex1.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Optimization Example 1 - Optimization of a quadratic objective function</h1>
21 // \author David Knezevic
22 // \date 2015
23 //
24 // In this example we demonstrate how to use OptimizationSystem to
25 // solve the optimization problem:
26 // min_U 0.5*U^T A U - U^T F
27 // We enforce Dirichlet constraints on U so that the minimization is
28 // well-posed. But note that we do not use OptimizationSystem's
29 // interface for imposing constraints in this case, so we can use an
30 // unconstrained solver (e.g. TAO's "Newton's method with line search"
31 // for unconstrained optimization is the default choice).
32 
33 // C++ include files that we need
34 #include <iostream>
35 
36 // libMesh includes
37 #include "libmesh/libmesh.h"
38 #include "libmesh/mesh.h"
39 #include "libmesh/mesh_generation.h"
40 #include "libmesh/exodusII_io.h"
41 #include "libmesh/equation_systems.h"
42 #include "libmesh/fe.h"
43 #include "libmesh/quadrature_gauss.h"
44 #include "libmesh/dof_map.h"
45 #include "libmesh/sparse_matrix.h"
46 #include "libmesh/numeric_vector.h"
47 #include "libmesh/dense_matrix.h"
48 #include "libmesh/dense_vector.h"
49 #include "libmesh/elem.h"
50 #include "libmesh/boundary_info.h"
51 #include "libmesh/string_to_enum.h"
52 #include "libmesh/getpot.h"
53 #include "libmesh/zero_function.h"
54 #include "libmesh/dirichlet_boundaries.h"
55 #include "libmesh/optimization_system.h"
56 #include "libmesh/optimization_solver.h"
57 
58 // Bring in everything from the libMesh namespace
59 using namespace libMesh;
60 
61 
70 {
71 private:
72 
77 
78 public:
79 
84 
90  void assemble_A_and_F();
91 
95  virtual Number objective (const NumericVector<Number> & soln,
96  OptimizationSystem & /*sys*/);
97 
101  virtual void gradient (const NumericVector<Number> & soln,
102  NumericVector<Number> & grad_f,
103  OptimizationSystem & /*sys*/);
104 
108  virtual void hessian (const NumericVector<Number> & soln,
109  SparseMatrix<Number> & H_f,
110  OptimizationSystem & /*sys*/);
111 
118 
124 };
125 
127  _sys(sys_in)
128 {}
129 
131 {
132  A_matrix->zero();
133  F_vector->zero();
134 
135  const MeshBase & mesh = _sys.get_mesh();
136 
137  const unsigned int dim = mesh.mesh_dimension();
138  const unsigned int u_var = _sys.variable_number ("u");
139 
140  const DofMap & dof_map = _sys.get_dof_map();
141  FEType fe_type = dof_map.variable_type(u_var);
142  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
143  QGauss qrule (dim, fe_type.default_quadrature_order());
144  fe->attach_quadrature_rule (&qrule);
145 
146  const std::vector<Real> & JxW = fe->get_JxW();
147  const std::vector<std::vector<Real>> & phi = fe->get_phi();
148  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
149 
150  std::vector<dof_id_type> dof_indices;
151 
154 
155  for (const auto & elem : mesh.active_local_element_ptr_range())
156  {
157  dof_map.dof_indices (elem, dof_indices);
158 
159  const unsigned int n_dofs = dof_indices.size();
160 
161  fe->reinit (elem);
162 
163  Ke.resize (n_dofs, n_dofs);
164  Fe.resize (n_dofs);
165 
166  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
167  {
168  for (unsigned int dof_i=0; dof_i<n_dofs; dof_i++)
169  {
170  for (unsigned int dof_j=0; dof_j<n_dofs; dof_j++)
171  {
172  Ke(dof_i, dof_j) += JxW[qp] * (dphi[dof_j][qp]* dphi[dof_i][qp]);
173  }
174  Fe(dof_i) += JxW[qp] * phi[dof_i][qp];
175  }
176  }
177 
178  // This will zero off-diagonal entries of Ke corresponding to
179  // Dirichlet dofs.
180  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
181 
182  // We want the diagonal of constrained dofs to be zero too
183  for (unsigned int local_dof_index=0; local_dof_index<n_dofs; local_dof_index++)
184  {
185  dof_id_type global_dof_index = dof_indices[local_dof_index];
186  if (dof_map.is_constrained_dof(global_dof_index))
187  {
188  Ke(local_dof_index, local_dof_index) = 0.;
189  }
190  }
191 
192  A_matrix->add_matrix (Ke, dof_indices);
193  F_vector->add_vector (Fe, dof_indices);
194  }
195 
196  A_matrix->close();
197  F_vector->close();
198 }
199 
201  OptimizationSystem & /*sys*/)
202 {
204 
205  A_matrix->vector_mult(*AxU, soln);
206 
207  Number
208  UTxAxU = AxU->dot(soln),
209  UTxF = F_vector->dot(soln);
210 
211  return 0.5 * UTxAxU - UTxF;
212 }
213 
215  NumericVector<Number> & grad_f,
216  OptimizationSystem & /*sys*/)
217 {
218  grad_f.zero();
219 
220  // Since we've enforced constraints on soln, A and F,
221  // this automatically sets grad_f to zero for constrained
222  // dofs.
223  A_matrix->vector_mult(grad_f, soln);
224  grad_f.add(-1, *F_vector);
225 }
226 
227 
229  SparseMatrix<Number> & H_f,
230  OptimizationSystem & /*sys*/)
231 {
232  H_f.zero();
233  H_f.add(1., *A_matrix);
234 }
235 
236 
237 int main (int argc, char ** argv)
238 {
239  LibMeshInit init (argc, argv);
240 
241 #ifndef LIBMESH_HAVE_PETSC_TAO
242 
243  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support");
244 
245 #elif !defined(LIBMESH_ENABLE_GHOSTED)
246 
247  libmesh_example_requires(false, "--enable-ghosted");
248 
249 #elif LIBMESH_USE_COMPLEX_NUMBERS
250 
251  // According to
252  // http://www.mcs.anl.gov/research/projects/tao/documentation/installation.html
253  // TAO & PETSc-complex are currently mutually exclusive
254  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support & real-numbers only");
255 
256 #endif
257 
258  // We use a 2D domain.
259  libmesh_example_requires(LIBMESH_DIM > 1, "--disable-1D-only");
260 
261  if (libMesh::on_command_line ("--use-eigen"))
262  {
263  libMesh::err << "This example requires an OptimizationSolver, and therefore does not "
264  << "support --use-eigen on the command line."
265  << std::endl;
266  return 0;
267  }
268 
269  GetPot infile("optimization_ex1.in");
270  const std::string approx_order = infile("approx_order", "FIRST");
271  const std::string fe_family = infile("fe_family", "LAGRANGE");
272  const unsigned int n_elem = infile("n_elem", 10);
273 
274  Mesh mesh(init.comm());
276  n_elem,
277  n_elem,
278  -1., 1.,
279  -1., 1.,
280  QUAD9);
281 
282  mesh.print_info();
283 
284  EquationSystems equation_systems (mesh);
285 
286  OptimizationSystem & system =
287  equation_systems.add_system<OptimizationSystem> ("Optimization");
288 
289  // The default is to use PETSc/Tao solvers, but let the user change
290  // the optimization solver package on the fly.
291  {
292  const std::string optimization_solver_type = infile("optimization_solver_type",
293  "PETSC_SOLVERS");
294  SolverPackage sp = Utility::string_to_enum<SolverPackage>(optimization_solver_type);
297  system.optimization_solver.reset(new_solver.release());
298  }
299 
300  // Set tolerances and maximum iteration counts directly on the optimization solver.
301  system.optimization_solver->max_objective_function_evaluations = 128;
302  system.optimization_solver->objective_function_relative_tolerance = 1.e-4;
303  system.optimization_solver->verbose = true;
304 
305  AssembleOptimization assemble_opt(system);
306 
307  unsigned int u_var = system.add_variable("u",
308  Utility::string_to_enum<Order> (approx_order),
309  Utility::string_to_enum<FEFamily>(fe_family));
310 
311  system.optimization_solver->objective_object = &assemble_opt;
312  system.optimization_solver->gradient_object = &assemble_opt;
313  system.optimization_solver->hessian_object = &assemble_opt;
314 
315  // system.matrix and system.rhs are used for the gradient and Hessian,
316  // so in this case we add an extra matrix and vector to store A and F.
317  // This makes it easy to write the code for evaluating the objective,
318  // gradient, and hessian.
319  system.add_matrix("A_matrix");
320  system.add_vector("F_vector");
321  assemble_opt.A_matrix = &system.get_matrix("A_matrix");
322  assemble_opt.F_vector = &system.get_vector("F_vector");
323 
324  // Apply Dirichlet constraints. This will be used to apply constraints
325  // to the objective function, gradient and Hessian.
326  std::set<boundary_id_type> boundary_ids;
327  boundary_ids.insert(0);
328  boundary_ids.insert(1);
329  boundary_ids.insert(2);
330  boundary_ids.insert(3);
331  std::vector<unsigned int> variables;
332  variables.push_back(u_var);
333  ZeroFunction<> zf;
334 
335  // Most DirichletBoundary users will want to supply a "locally
336  // indexed" functor
337  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
339  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
340 
341 
342  equation_systems.init();
343  equation_systems.print_info();
344 
345  assemble_opt.assemble_A_and_F();
346 
347  // We need to close the matrix so that we can use it to store the
348  // Hessian during the solve.
349  system.matrix->close();
350  system.solve();
351 
352  // Print convergence information
353  system.optimization_solver->print_converged_reason();
354 
355 #ifdef LIBMESH_HAVE_EXODUS_API
356  std::stringstream filename;
357  ExodusII_IO (mesh).write_equation_systems("optimization_soln.exo",
358  equation_systems);
359 #endif
360 
361  return 0;
362 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
OptimizationSystem & _sys
Keep a reference to the OptimizationSystem.
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
Inserts the entries of v in *this at the locations specified by v.
This is the EquationSystems class.
virtual Number objective(const NumericVector< Number > &soln, OptimizationSystem &)
Evaluate the objective function.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
Abstract base class to be used to calculate the objective function for optimization.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
AssembleOptimization(OptimizationSystem &sys_in)
Constructor.
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
SolverPackage
Defines an enum for various linear solver packages.
virtual void solve() libmesh_override
Solves the optimization problem.
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
MeshBase & mesh
UniquePtr< OptimizationSolver< Number > > optimization_solver
The OptimizationSolver that is used for performing the optimization.
NumericVector< Number > * F_vector
Vector for storing F.
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
int main(int argc, char **argv)
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
SparseMatrix< Number > * A_matrix
Sparse matrix for storing the matrix A.
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1734
This class encapsulate all functionality required for assembling the objective function, gradient, and hessian.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual void zero()=0
Set all entries to 0.
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Order default_quadrature_order() const
Definition: fe_type.h:332
This System subclass enables us to assemble an objective function, gradient, Hessian and bounds for o...
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
Abstract base class to be used to calculate the gradient of an objective function.
virtual void hessian(const NumericVector< Number > &soln, SparseMatrix< Number > &H_f, OptimizationSystem &)
Evaluate the Hessian.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
This base class can be inherited from to provide interfaces to optimization solvers from different pa...
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
virtual UniquePtr< NumericVector< T > > zero_clone() const =0
SparseMatrix< Number > * matrix
The system matrix.
virtual void gradient(const NumericVector< Number > &soln, NumericVector< Number > &grad_f, OptimizationSystem &)
Evaluate the gradient.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
Abstract base class to be used to calculate the Hessian of an objective function. ...
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual T dot(const NumericVector< T > &v) const =0
virtual void init()
Initialize all the systems.
void assemble_A_and_F()
The optimization problem we consider here is: min_U 0.5*U^T A U - U^T F.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917