libMesh
Classes | Functions
optimization_ex1.C File Reference

Go to the source code of this file.

Classes

class  AssembleOptimization
 This class encapsulate all functionality required for assembling the objective function, gradient, and hessian. More...
 

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 237 of file optimization_ex1.C.

References AssembleOptimization::A_matrix, libMesh::DofMap::add_dirichlet_boundary(), libMesh::ImplicitSystem::add_matrix(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::System::add_vector(), AssembleOptimization::assemble_A_and_F(), libMesh::MeshTools::Generation::build_square(), libMesh::SparseMatrix< T >::close(), libMesh::LibMeshInit::comm(), libMesh::err, AssembleOptimization::F_vector, libMesh::System::get_dof_map(), libMesh::ImplicitSystem::get_matrix(), libMesh::System::get_vector(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::NumericVector< T >::insert(), libMesh::LOCAL_VARIABLE_ORDER, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshTools::n_elem(), libMesh::on_command_line(), libMesh::OptimizationSystem::optimization_solver, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::OptimizationSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().

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 }
OStreamProxy err
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.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
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
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
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.
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
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
SparseMatrix< Number > * matrix
The system matrix.
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
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