optimization_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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
25 #include "libmesh/sparse_matrix.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/dof_map.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  numeric_index_type n_eq_constraints =
94  cast_int<numeric_index_type>(constraint_jac_sparsity.size());
95 
96  // Assign rows to each processor as evenly as possible
97  unsigned int n_procs = comm().size();
98  numeric_index_type n_local_rows = n_eq_constraints / n_procs;
99  if (comm().rank() < (n_eq_constraints % n_procs))
100  n_local_rows++;
101 
102  C_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
103  lambda_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
104 
105  // Get the maximum number of non-zeros per row
106  numeric_index_type max_nnz = 0;
107  for (numeric_index_type i=0; i<n_eq_constraints; i++)
108  {
109  numeric_index_type nnz =
110  cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());
111  if (nnz > max_nnz)
112  max_nnz = nnz;
113  }
114 
115  C_eq_jac->init(n_eq_constraints,
116  get_dof_map().n_dofs(),
117  n_local_rows,
119  max_nnz,
120  max_nnz);
121 
122  eq_constraint_jac_sparsity = constraint_jac_sparsity;
123 }
124 
125 
127 initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
128 {
129  numeric_index_type n_ineq_constraints =
130  cast_int<numeric_index_type>(constraint_jac_sparsity.size());
131 
132  // Assign rows to each processor as evenly as possible
133  unsigned int n_procs = comm().size();
134  numeric_index_type n_local_rows = n_ineq_constraints / n_procs;
135  if (comm().rank() < (n_ineq_constraints % n_procs))
136  n_local_rows++;
137 
138  C_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
139  lambda_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
140 
141  // Get the maximum number of non-zeros per row
142  numeric_index_type max_nnz = 0;
143  for (numeric_index_type i=0; i<n_ineq_constraints; i++)
144  {
145  numeric_index_type nnz =
146  cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());
147  if (nnz > max_nnz)
148  max_nnz = nnz;
149  }
150 
151  C_ineq_jac->init(n_ineq_constraints,
152  get_dof_map().n_dofs(),
153  n_local_rows,
155  max_nnz,
156  max_nnz);
157 
158  ineq_constraint_jac_sparsity = constraint_jac_sparsity;
159 }
160 
161 
163 {
164  START_LOG("solve()", "OptimizationSystem");
165 
166  optimization_solver->init();
167  optimization_solver->solve ();
168 
169  STOP_LOG("solve()", "OptimizationSystem");
170 
171  this->update();
172 }
173 
174 
175 } // namespace libMesh
virtual void init_data() override
Manages multiples systems of equations.
void initialize_inequality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
std::unique_ptr< NumericVector< Number > > lambda_ineq
processor_id_type size() const
Definition: communicator.h:175
virtual void solve() override
virtual void clear() override
virtual void init_data() override
virtual void reinit() override
std::unique_ptr< NumericVector< Number > > C_eq
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
const Parallel::Communicator & comm() const
dof_id_type n_local_dofs() const
Definition: system.C:187
std::unique_ptr< SparseMatrix< Number > > C_ineq_jac
dof_id_type n_dofs() const
Definition: system.C:150
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:661
dof_id_type numeric_index_type
Definition: id_types.h:92
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
virtual void reinit() override
std::vector< std::set< numeric_index_type > > ineq_constraint_jac_sparsity
std::vector< std::set< numeric_index_type > > eq_constraint_jac_sparsity
std::unique_ptr< SparseMatrix< Number > > C_eq_jac
std::unique_ptr< OptimizationSolver< Number > > optimization_solver
std::unique_ptr< NumericVector< Number > > lambda_eq
virtual void update()
Definition: system.C:408
void initialize_equality_constraints_storage(const std::vector< std::set< numeric_index_type >> &constraint_jac_sparsity)
virtual void clear() override
OptimizationSystem(EquationSystems &es, const std::string &name, const unsigned int number)
const DofMap & get_dof_map() const
Definition: system.h:2049
std::unique_ptr< NumericVector< Number > > C_ineq