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