libMesh
optimization_system.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 // C++ includes
21 
22 // Local includes
23 #include "libmesh/equation_systems.h"
24 #include "libmesh/libmesh_logging.h"
25 #include "libmesh/sparse_matrix.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/dof_map.h"
28 #include "libmesh/optimization_solver.h"
29 #include "libmesh/optimization_system.h"
30 
31 namespace libMesh
32 {
33 
34 // ------------------------------------------------------------
35 // OptimizationSystem implementation
37  const std::string & name_in,
38  const unsigned int number_in) :
39 
40  Parent(es, name_in, number_in),
41  optimization_solver(OptimizationSolver<Number>::build(*this)),
42  C_eq(NumericVector<Number>::build(this->comm())),
43  C_eq_jac(SparseMatrix<Number>::build(this->comm())),
44  C_ineq(NumericVector<Number>::build(this->comm())),
45  C_ineq_jac(SparseMatrix<Number>::build(this->comm())),
46  lambda_eq(NumericVector<Number>::build(this->comm())),
47  lambda_ineq(NumericVector<Number>::build(this->comm()))
48 {
49 }
50 
51 
52 
54 {
55  // Clear data
56  this->clear();
57 }
58 
59 
60 
62 {
63  // clear the optimization solver
64  optimization_solver->clear();
65 
66  // clear the parent data
67  Parent::clear();
68 }
69 
70 
72 {
73  this->add_vector("lower_bounds");
74  this->add_vector("upper_bounds");
75 
77 
78  optimization_solver->clear();
79 }
80 
81 
83 {
84  optimization_solver->clear();
85 
87 }
88 
89 
91 initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
92 {
93  unsigned int n_eq_constraints = constraint_jac_sparsity.size();
94 
95  // Assign rows to each processor as evenly as possible
96  unsigned int n_procs = comm().size();
97  unsigned int n_local_rows = n_eq_constraints / n_procs;
98  if (comm().rank() < (n_eq_constraints % n_procs))
99  n_local_rows++;
100 
101  C_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
102  lambda_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
103 
104  // Get the maximum number of non-zeros per row
105  unsigned int max_nnz = 0;
106  for (unsigned int i=0; i<n_eq_constraints; i++)
107  {
108  unsigned int nnz = constraint_jac_sparsity[i].size();
109  if (nnz > max_nnz)
110  max_nnz = nnz;
111  }
112 
113  C_eq_jac->init(n_eq_constraints,
114  get_dof_map().n_dofs(),
115  n_local_rows,
117  max_nnz,
118  max_nnz);
119 
120  eq_constraint_jac_sparsity = constraint_jac_sparsity;
121 }
122 
123 
125 initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
126 {
127  unsigned int n_ineq_constraints = constraint_jac_sparsity.size();
128 
129  // Assign rows to each processor as evenly as possible
130  unsigned int n_procs = comm().size();
131  unsigned int n_local_rows = n_ineq_constraints / n_procs;
132  if (comm().rank() < (n_ineq_constraints % n_procs))
133  n_local_rows++;
134 
135  C_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
136  lambda_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
137 
138  // Get the maximum number of non-zeros per row
139  unsigned int max_nnz = 0;
140  for (unsigned int i=0; i<n_ineq_constraints; i++)
141  {
142  unsigned int nnz = constraint_jac_sparsity[i].size();
143  if (nnz > max_nnz)
144  max_nnz = nnz;
145  }
146 
147  C_ineq_jac->init(n_ineq_constraints,
148  get_dof_map().n_dofs(),
149  n_local_rows,
151  max_nnz,
152  max_nnz);
153 
154  ineq_constraint_jac_sparsity = constraint_jac_sparsity;
155 }
156 
157 
159 {
160  START_LOG("solve()", "OptimizationSystem");
161 
162  optimization_solver->init();
163  optimization_solver->solve ();
164 
165  STOP_LOG("solve()", "OptimizationSystem");
166 
167  this->update();
168 }
169 
170 
171 } // namespace libMesh
This is the EquationSystems class.
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
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.
unsigned int size() const
Definition: parallel.h:726
UniquePtr< NumericVector< Number > > lambda_ineq
virtual void solve() libmesh_override
Solves the optimization problem.
UniquePtr< OptimizationSolver< Number > > optimization_solver
The OptimizationSolver that is used for performing the optimization.
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
virtual void init_data() libmesh_override
Initializes new data members of the system.
UniquePtr< NumericVector< Number > > lambda_eq
Vectors to store the dual variables associated with equality and inequality constraints.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
Generic sparse matrix.
Definition: dof_map.h:65
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 is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
UniquePtr< SparseMatrix< Number > > C_eq_jac
The sparse matrix that stores the Jacobian of C_eq.
std::vector< std::set< numeric_index_type > > ineq_constraint_jac_sparsity
std::vector< std::set< numeric_index_type > > eq_constraint_jac_sparsity
A copy of the equality and inequality constraint Jacobian sparsity patterns.
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
dof_id_type n_local_dofs() const
Definition: system.C:185
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
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.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
virtual ~OptimizationSystem()
Destructor.
This base class can be inherited from to provide interfaces to optimization solvers from different pa...
const Parallel::Communicator & comm() const
dof_id_type n_dofs() const
Definition: system.C:148
UniquePtr< NumericVector< Number > > C_eq
The vector that stores equality constraints.
OptimizationSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
UniquePtr< SparseMatrix< Number > > C_ineq_jac
The sparse matrix that stores the Jacobian of C_ineq.
UniquePtr< NumericVector< Number > > C_ineq
The vector that stores inequality constraints.