1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
OptimizationSystem(EquationSystems & es,const std::string & name_in,const unsigned int number_in)36 OptimizationSystem::OptimizationSystem (EquationSystems & es,
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
~OptimizationSystem()53 OptimizationSystem::~OptimizationSystem ()
54 {
55 // Clear data
56 this->clear();
57 }
58
59
60
clear()61 void OptimizationSystem::clear ()
62 {
63 // clear the optimization solver
64 optimization_solver->clear();
65
66 // clear the parent data
67 Parent::clear();
68 }
69
70
init_data()71 void OptimizationSystem::init_data ()
72 {
73 this->add_vector("lower_bounds");
74 this->add_vector("upper_bounds");
75
76 Parent::init_data();
77
78 optimization_solver->clear();
79 }
80
81
reinit()82 void OptimizationSystem::reinit ()
83 {
84 optimization_solver->clear();
85
86 Parent::reinit();
87 }
88
89
90 void OptimizationSystem::
initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)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,
118 get_dof_map().n_local_dofs(),
119 max_nnz,
120 max_nnz);
121
122 eq_constraint_jac_sparsity = constraint_jac_sparsity;
123 }
124
125
126 void OptimizationSystem::
initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)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,
154 get_dof_map().n_local_dofs(),
155 max_nnz,
156 max_nnz);
157
158 ineq_constraint_jac_sparsity = constraint_jac_sparsity;
159 }
160
161
solve()162 void OptimizationSystem::solve ()
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
176