libMesh
Classes | Functions
optimization_ex2.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 395 of file optimization_ex2.C.

References AssembleOptimization::A_matrix, 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(), AssembleOptimization::F_vector, libMesh::ImplicitSystem::get_matrix(), libMesh::System::get_vector(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshTools::n_elem(), libMesh::OptimizationSystem::optimization_solver, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::Parallel::Communicator::size(), libMesh::OptimizationSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().

396 {
397  LibMeshInit init (argc, argv);
398 
399 #ifndef LIBMESH_HAVE_PETSC_TAO
400 
401  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support");
402 
403 #elif LIBMESH_USE_COMPLEX_NUMBERS
404 
405  // According to
406  // http://www.mcs.anl.gov/research/projects/tao/documentation/installation.html
407  // TAO & PETSc-complex are currently mutually exclusive
408  libmesh_example_requires(false, "PETSc >= 3.5.0 with built-in TAO support & real-numbers only");
409 
410 #endif
411 
412  // We use a 2D domain.
413  libmesh_example_requires(LIBMESH_DIM > 1, "--disable-1D-only");
414 
415  // TAO is giving us problems in parallel?
416  if (init.comm().size() != 1)
417  {
418  libMesh::out << "This example can currently only be run in serial." << std::endl;
419  return 77;
420  }
421 
422  GetPot infile("optimization_ex2.in");
423  const std::string approx_order = infile("approx_order", "FIRST");
424  const std::string fe_family = infile("fe_family", "LAGRANGE");
425  const unsigned int n_elem = infile("n_elem", 10);
426 
427  Mesh mesh(init.comm());
429  n_elem,
430  n_elem,
431  -1., 1.,
432  -1., 1.,
433  QUAD9);
434 
435  mesh.print_info();
436 
437  EquationSystems equation_systems (mesh);
438 
439  OptimizationSystem & system =
440  equation_systems.add_system<OptimizationSystem> ("Optimization");
441 
442  // The default is to use PETSc/Tao solvers, but let the user change
443  // the optimization solver package on the fly.
444  {
445  const std::string optimization_solver_type = infile("optimization_solver_type",
446  "PETSC_SOLVERS");
447  SolverPackage sp = Utility::string_to_enum<SolverPackage>(optimization_solver_type);
450  system.optimization_solver.reset(new_solver.release());
451  }
452 
453  // Set tolerances and maximum iteration counts directly on the optimization solver.
454  system.optimization_solver->max_objective_function_evaluations = 128;
455  system.optimization_solver->objective_function_relative_tolerance = 1.e-4;
456  system.optimization_solver->verbose = true;
457 
458  AssembleOptimization assemble_opt(system);
459 
460  system.add_variable("u",
461  Utility::string_to_enum<Order> (approx_order),
462  Utility::string_to_enum<FEFamily>(fe_family));
463 
464  system.optimization_solver->objective_object = &assemble_opt;
465  system.optimization_solver->gradient_object = &assemble_opt;
466  system.optimization_solver->hessian_object = &assemble_opt;
467  system.optimization_solver->equality_constraints_object = &assemble_opt;
468  system.optimization_solver->equality_constraints_jacobian_object = &assemble_opt;
469  system.optimization_solver->inequality_constraints_object = &assemble_opt;
470  system.optimization_solver->inequality_constraints_jacobian_object = &assemble_opt;
471  system.optimization_solver->lower_and_upper_bounds_object = &assemble_opt;
472 
473  // system.matrix and system.rhs are used for the gradient and Hessian,
474  // so in this case we add an extra matrix and vector to store A and F.
475  // This makes it easy to write the code for evaluating the objective,
476  // gradient, and hessian.
477  system.add_matrix("A_matrix");
478  system.add_vector("F_vector");
479  assemble_opt.A_matrix = &system.get_matrix("A_matrix");
480  assemble_opt.F_vector = &system.get_vector("F_vector");
481 
482 
483  equation_systems.init();
484  equation_systems.print_info();
485 
486  assemble_opt.assemble_A_and_F();
487 
488  {
489  std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
490  std::set<numeric_index_type> sparsity_row;
491  sparsity_row.insert(17);
492  constraint_jac_sparsity.push_back(sparsity_row);
493  sparsity_row.clear();
494 
495  sparsity_row.insert(23);
496  constraint_jac_sparsity.push_back(sparsity_row);
497  sparsity_row.clear();
498 
499  sparsity_row.insert(98);
500  sparsity_row.insert(185);
501  constraint_jac_sparsity.push_back(sparsity_row);
502 
503  system.initialize_equality_constraints_storage(constraint_jac_sparsity);
504  }
505 
506  {
507  std::vector<std::set<numeric_index_type>> constraint_jac_sparsity;
508  std::set<numeric_index_type> sparsity_row;
509  sparsity_row.insert(200);
510  sparsity_row.insert(201);
511  constraint_jac_sparsity.push_back(sparsity_row);
512 
513  system.initialize_inequality_constraints_storage(constraint_jac_sparsity);
514  }
515 
516  // We need to close the matrix so that we can use it to store the
517  // Hessian during the solve.
518  system.matrix->close();
519  system.solve();
520 
521  // Print convergence information
522  system.optimization_solver->print_converged_reason();
523 
524  // Write the solution to file if the optimization solver converged
525  if (system.optimization_solver->get_converged_reason() > 0)
526  {
527  std::stringstream filename;
528  ExodusII_IO (mesh).write_equation_systems("optimization_soln.exo",
529  equation_systems);
530  }
531 
532  return 0;
533 }
This is the EquationSystems class.
void initialize_inequality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
Initialize storage for the inequality constraints, as per initialize_equality_constraints_storage.
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
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.
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
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...
void initialize_equality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
Initialize storage for the equality constraints, and the corresponding Jacobian.
This base class can be inherited from to provide interfaces to optimization solvers from different pa...
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
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