Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local includes
21 : #include "libmesh/equation_systems.h"
22 : #include "libmesh/libmesh_logging.h"
23 : #include "libmesh/sparse_matrix.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/dof_map.h"
26 : #include "libmesh/optimization_solver.h"
27 : #include "libmesh/optimization_system.h"
28 :
29 : namespace libMesh
30 : {
31 :
32 71 : OptimizationSystem::OptimizationSystem (EquationSystems & es,
33 : const std::string & name_in,
34 71 : const unsigned int number_in) :
35 :
36 : Parent(es, name_in, number_in),
37 67 : optimization_solver(OptimizationSolver<Number>::build(*this)),
38 67 : C_eq(NumericVector<Number>::build(this->comm())),
39 67 : C_eq_jac(SparseMatrix<Number>::build(this->comm())),
40 67 : C_ineq(NumericVector<Number>::build(this->comm())),
41 67 : C_ineq_jac(SparseMatrix<Number>::build(this->comm())),
42 67 : lambda_eq(NumericVector<Number>::build(this->comm())),
43 138 : lambda_ineq(NumericVector<Number>::build(this->comm()))
44 : {
45 71 : }
46 :
47 :
48 :
49 134 : OptimizationSystem::~OptimizationSystem () = default;
50 :
51 :
52 :
53 0 : void OptimizationSystem::clear ()
54 : {
55 : // clear the optimization solver
56 0 : optimization_solver->clear();
57 :
58 : // clear the parent data
59 0 : Parent::clear();
60 0 : }
61 :
62 :
63 71 : void OptimizationSystem::init_data ()
64 : {
65 71 : this->add_vector("lower_bounds");
66 71 : this->add_vector("upper_bounds");
67 :
68 71 : Parent::init_data();
69 :
70 71 : optimization_solver->clear();
71 71 : }
72 :
73 :
74 0 : void OptimizationSystem::reinit ()
75 : {
76 0 : optimization_solver->clear();
77 :
78 0 : Parent::reinit();
79 0 : }
80 :
81 :
82 1 : void OptimizationSystem::
83 : initialize_equality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
84 : {
85 : numeric_index_type n_eq_constraints =
86 0 : cast_int<numeric_index_type>(constraint_jac_sparsity.size());
87 :
88 : // Assign rows to each processor as evenly as possible
89 0 : unsigned int n_procs = comm().size();
90 1 : numeric_index_type n_local_rows = n_eq_constraints / n_procs;
91 1 : if (comm().rank() < (n_eq_constraints % n_procs))
92 0 : n_local_rows++;
93 :
94 1 : C_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
95 1 : lambda_eq->init(n_eq_constraints, n_local_rows, false, PARALLEL);
96 :
97 : // Get the maximum number of non-zeros per row
98 0 : numeric_index_type max_nnz = 0;
99 4 : for (numeric_index_type i=0; i<n_eq_constraints; i++)
100 : {
101 : numeric_index_type nnz =
102 0 : cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());
103 0 : if (nnz > max_nnz)
104 0 : max_nnz = nnz;
105 : }
106 :
107 1 : C_eq_jac->init(n_eq_constraints,
108 0 : get_dof_map().n_dofs(),
109 : n_local_rows,
110 0 : get_dof_map().n_local_dofs(),
111 : max_nnz,
112 0 : max_nnz);
113 :
114 1 : eq_constraint_jac_sparsity = constraint_jac_sparsity;
115 1 : }
116 :
117 :
118 1 : void OptimizationSystem::
119 : initialize_inequality_constraints_storage(const std::vector<std::set<numeric_index_type>> & constraint_jac_sparsity)
120 : {
121 : numeric_index_type n_ineq_constraints =
122 0 : cast_int<numeric_index_type>(constraint_jac_sparsity.size());
123 :
124 : // Assign rows to each processor as evenly as possible
125 0 : unsigned int n_procs = comm().size();
126 1 : numeric_index_type n_local_rows = n_ineq_constraints / n_procs;
127 1 : if (comm().rank() < (n_ineq_constraints % n_procs))
128 0 : n_local_rows++;
129 :
130 1 : C_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
131 1 : lambda_ineq->init(n_ineq_constraints, n_local_rows, false, PARALLEL);
132 :
133 : // Get the maximum number of non-zeros per row
134 0 : numeric_index_type max_nnz = 0;
135 2 : for (numeric_index_type i=0; i<n_ineq_constraints; i++)
136 : {
137 : numeric_index_type nnz =
138 0 : cast_int<numeric_index_type>(constraint_jac_sparsity[i].size());
139 0 : if (nnz > max_nnz)
140 0 : max_nnz = nnz;
141 : }
142 :
143 1 : C_ineq_jac->init(n_ineq_constraints,
144 0 : get_dof_map().n_dofs(),
145 : n_local_rows,
146 0 : get_dof_map().n_local_dofs(),
147 : max_nnz,
148 0 : max_nnz);
149 :
150 1 : ineq_constraint_jac_sparsity = constraint_jac_sparsity;
151 1 : }
152 :
153 :
154 71 : void OptimizationSystem::solve ()
155 : {
156 4 : LOG_SCOPE("solve()", "OptimizationSystem");
157 :
158 71 : optimization_solver->init();
159 71 : optimization_solver->solve ();
160 :
161 71 : this->update();
162 71 : }
163 :
164 :
165 : } // namespace libMesh
|
