Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/nonlinear_implicit_system.h"
22 : #include "libmesh/diff_solver.h"
23 : #include "libmesh/equation_systems.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/nonlinear_solver.h"
26 : #include "libmesh/sparse_matrix.h"
27 : #include "libmesh/static_condensation.h"
28 : #include "libmesh/static_condensation_preconditioner.h"
29 :
30 : namespace libMesh
31 : {
32 :
33 1470 : NonlinearImplicitSystem::NonlinearImplicitSystem (EquationSystems & es,
34 : const std::string & name_in,
35 1470 : const unsigned int number_in) :
36 :
37 : Parent (es, name_in, number_in),
38 1386 : nonlinear_solver (NonlinearSolver<Number>::build(*this)),
39 1386 : diff_solver (),
40 1386 : _n_nonlinear_iterations (0),
41 1386 : _final_nonlinear_residual (1.e20),
42 1470 : _operator_matrix (nullptr)
43 : {
44 : // Set default parameters
45 : // These were chosen to match the Petsc defaults
46 1470 : es.parameters.set<Real> ("linear solver tolerance") = 1e-5;
47 1470 : es.parameters.set<Real> ("linear solver minimum tolerance") = 1e-5;
48 1470 : es.parameters.set<unsigned int>("linear solver maximum iterations") = 10000;
49 :
50 1470 : es.parameters.set<unsigned int>("nonlinear solver maximum iterations") = 50;
51 1470 : es.parameters.set<unsigned int>("nonlinear solver maximum function evaluations") = 10000;
52 :
53 1470 : es.parameters.set<Real>("nonlinear solver absolute residual tolerance") = 1e-35;
54 1470 : es.parameters.set<Real>("nonlinear solver relative residual tolerance") = 1e-8;
55 1470 : es.parameters.set<Real>("nonlinear solver divergence tolerance") = 1e+4;
56 1470 : es.parameters.set<Real>("nonlinear solver absolute step tolerance") = 1e-8;
57 1470 : es.parameters.set<Real>("nonlinear solver relative step tolerance") = 1e-8;
58 :
59 1470 : es.parameters.set<bool>("reuse preconditioner") = false;
60 1470 : es.parameters.set<unsigned int>("reuse preconditioner maximum linear iterations") = 1;
61 :
62 1470 : if (this->has_static_condensation())
63 280 : this->setup_static_condensation_preconditioner(*nonlinear_solver);
64 1470 : }
65 :
66 :
67 :
68 2508 : NonlinearImplicitSystem::~NonlinearImplicitSystem () = default;
69 :
70 :
71 :
72 0 : void NonlinearImplicitSystem::create_static_condensation()
73 : {
74 0 : Parent::create_static_condensation();
75 0 : this->setup_static_condensation_preconditioner(*nonlinear_solver);
76 0 : }
77 :
78 :
79 :
80 280 : void NonlinearImplicitSystem::clear ()
81 : {
82 : // clear the nonlinear solver
83 280 : nonlinear_solver->clear();
84 :
85 : // FIXME - this is necessary for petsc_auto_fieldsplit
86 : // nonlinear_solver->init_names(*this);
87 :
88 : // clear the parent data
89 280 : Parent::clear();
90 :
91 : // And restore any StaticCondensation to defaults
92 280 : if (this->has_static_condensation())
93 0 : this->setup_static_condensation_preconditioner(*nonlinear_solver);
94 280 : }
95 :
96 :
97 :
98 0 : void NonlinearImplicitSystem::reinit ()
99 : {
100 : // re-initialize the nonlinear solver interface
101 0 : nonlinear_solver->clear();
102 :
103 : // force the solver to get a new preconditioner, in
104 : // case reuse was set
105 0 : nonlinear_solver->force_new_preconditioner();
106 :
107 : // FIXME - this is necessary for petsc_auto_fieldsplit
108 : // nonlinear_solver->init_names(*this);
109 :
110 0 : if (diff_solver.get())
111 0 : diff_solver->reinit();
112 :
113 : // initialize parent data
114 0 : Parent::reinit();
115 0 : }
116 :
117 :
118 :
119 29400 : void NonlinearImplicitSystem::set_solver_parameters ()
120 : {
121 : // Get a reference to the EquationSystems
122 : const EquationSystems & es =
123 1680 : this->get_equation_systems();
124 :
125 : // Get the user-specified nonlinear solver tolerances
126 57960 : const unsigned int maxits = parameters.have_parameter<unsigned int>("nonlinear solver maximum iterations") ?
127 0 : parameters.get<unsigned int>("nonlinear solver maximum iterations") :
128 57960 : es.parameters.get<unsigned int>("nonlinear solver maximum iterations");
129 :
130 57960 : const unsigned int maxfuncs = parameters.have_parameter<unsigned int>("nonlinear solver maximum function evaluations") ?
131 0 : parameters.get<unsigned int>("nonlinear solver maximum function evaluations") :
132 57960 : es.parameters.get<unsigned int>("nonlinear solver maximum function evaluations");
133 :
134 57960 : const double abs_resid_tol = parameters.have_parameter<Real>("nonlinear solver absolute residual tolerance") ?
135 27720 : double(parameters.get<Real>("nonlinear solver absolute residual tolerance")) :
136 57960 : double(es.parameters.get<Real>("nonlinear solver absolute residual tolerance"));
137 :
138 57960 : const double rel_resid_tol = parameters.have_parameter<Real>("nonlinear solver relative residual tolerance") ?
139 0 : double(parameters.get<Real>("nonlinear solver relative residual tolerance")) :
140 57960 : double(es.parameters.get<Real>("nonlinear solver relative residual tolerance"));
141 :
142 57960 : const double div_tol = parameters.have_parameter<Real>("nonlinear solver divergence tolerance") ?
143 0 : double(parameters.get<Real>("nonlinear solver divergence tolerance")) :
144 57960 : double(es.parameters.get<Real>("nonlinear solver divergence tolerance"));
145 :
146 57960 : const double abs_step_tol = parameters.have_parameter<Real>("nonlinear solver absolute step tolerance") ?
147 0 : double(parameters.get<Real>("nonlinear solver absolute step tolerance")) :
148 57960 : double(es.parameters.get<Real>("nonlinear solver absolute step tolerance"));
149 :
150 57960 : const double rel_step_tol = parameters.have_parameter<Real>("nonlinear solver relative step tolerance")?
151 0 : double(parameters.get<Real>("nonlinear solver relative step tolerance")) :
152 29400 : double(es.parameters.get<Real>("nonlinear solver relative step tolerance"));
153 :
154 : // Get the user-specified linear solver tolerances
155 29400 : const auto [maxlinearits, linear_tol] = this->Parent::get_linear_solve_parameters();
156 :
157 57960 : const double linear_min_tol = parameters.have_parameter<Real>("linear solver minimum tolerance") ?
158 0 : double(parameters.get<Real>("linear solver minimum tolerance")) :
159 57960 : double(es.parameters.get<Real>("linear solver minimum tolerance"));
160 :
161 29400 : const bool reuse_preconditioner = parameters.have_parameter<unsigned int>("reuse preconditioner") ?
162 0 : parameters.get<unsigned int>("reuse preconditioner") :
163 29400 : es.parameters.get<bool>("reuse preconditioner");
164 : const unsigned int reuse_preconditioner_max_linear_its =
165 57960 : parameters.have_parameter<unsigned int>("reuse preconditioner maximum linear iterations") ?
166 0 : parameters.get<unsigned int>("reuse preconditioner maximum linear iterations") :
167 57960 : es.parameters.get<unsigned int>("reuse preconditioner maximum linear iterations");
168 :
169 : // Set all the parameters on the NonlinearSolver
170 29400 : nonlinear_solver->max_nonlinear_iterations = maxits;
171 29400 : nonlinear_solver->max_function_evaluations = maxfuncs;
172 29400 : nonlinear_solver->absolute_residual_tolerance = abs_resid_tol;
173 29400 : nonlinear_solver->relative_residual_tolerance = rel_resid_tol;
174 29400 : nonlinear_solver->divergence_tolerance = div_tol;
175 29400 : nonlinear_solver->absolute_step_tolerance = abs_step_tol;
176 29400 : nonlinear_solver->relative_step_tolerance = rel_step_tol;
177 29400 : nonlinear_solver->max_linear_iterations = maxlinearits;
178 29400 : nonlinear_solver->initial_linear_tolerance = linear_tol;
179 29400 : nonlinear_solver->minimum_linear_tolerance = linear_min_tol;
180 29400 : nonlinear_solver->set_reuse_preconditioner(reuse_preconditioner);
181 29400 : nonlinear_solver->set_reuse_preconditioner_max_linear_its(reuse_preconditioner_max_linear_its);
182 :
183 29400 : if (diff_solver.get())
184 : {
185 0 : diff_solver->max_nonlinear_iterations = maxits;
186 0 : diff_solver->absolute_residual_tolerance = abs_resid_tol;
187 0 : diff_solver->relative_residual_tolerance = rel_resid_tol;
188 0 : diff_solver->absolute_step_tolerance = abs_step_tol;
189 0 : diff_solver->relative_step_tolerance = rel_step_tol;
190 0 : diff_solver->max_linear_iterations = maxlinearits;
191 0 : diff_solver->initial_linear_tolerance = linear_tol;
192 0 : diff_solver->minimum_linear_tolerance = linear_min_tol;
193 : }
194 29400 : }
195 :
196 :
197 :
198 29400 : void NonlinearImplicitSystem::solve ()
199 : {
200 : // Log how long the nonlinear solve takes.
201 1680 : LOG_SCOPE("solve()", "System");
202 :
203 29400 : this->set_solver_parameters();
204 :
205 29400 : if (diff_solver.get())
206 : {
207 0 : diff_solver->solve();
208 :
209 : // Store the number of nonlinear iterations required to
210 : // solve and the final residual.
211 0 : _n_nonlinear_iterations = diff_solver->total_outer_iterations();
212 0 : _final_nonlinear_residual = 0.; // FIXME - support this!
213 : }
214 : else
215 : {
216 29400 : if (this->prefix_with_name())
217 0 : nonlinear_solver->init(this->prefix().c_str());
218 : else
219 29400 : nonlinear_solver->init();
220 :
221 : // FIXME - this is necessary for petsc_auto_fieldsplit
222 : // nonlinear_solver->init_names(*this);
223 :
224 : // Solve the nonlinear system.
225 : // Store the number of nonlinear iterations required to
226 : // solve and the final residual.
227 : //
228 : // If a distinct Jacobian operator matrix has been registered (see
229 : // set_operator_matrix()), use it as Amat while *matrix remains the preconditioning matrix
230 : // (Pmat); otherwise use the ordinary single-matrix solve.
231 29400 : if (_operator_matrix)
232 0 : std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
233 0 : nonlinear_solver->solve (*_operator_matrix, *matrix, *solution, *rhs,
234 0 : nonlinear_solver->relative_residual_tolerance,
235 0 : nonlinear_solver->max_linear_iterations);
236 : else
237 29400 : std::tie(_n_nonlinear_iterations, _final_nonlinear_residual) =
238 30240 : nonlinear_solver->solve (*matrix, *solution, *rhs,
239 840 : nonlinear_solver->relative_residual_tolerance,
240 3360 : nonlinear_solver->max_linear_iterations);
241 : }
242 :
243 : // Update the system after the solve
244 29400 : this->update();
245 29400 : }
246 :
247 :
248 :
249 0 : std::pair<unsigned int, Real> NonlinearImplicitSystem::get_linear_solve_parameters() const
250 : {
251 0 : if (diff_solver.get())
252 0 : return std::make_pair(this->diff_solver->max_linear_iterations,
253 0 : this->diff_solver->relative_residual_tolerance);
254 0 : return std::make_pair(this->nonlinear_solver->max_linear_iterations,
255 0 : this->nonlinear_solver->relative_residual_tolerance);
256 : }
257 :
258 :
259 :
260 0 : void NonlinearImplicitSystem::assembly(bool get_residual,
261 : bool get_jacobian,
262 : bool /*apply_heterogeneous_constraints*/,
263 : bool /*apply_no_constraints*/)
264 : {
265 : // Get current_local_solution in sync
266 0 : this->update();
267 :
268 : //-----------------------------------------------------------------------------
269 : // if the user has provided both function pointers and objects only the pointer
270 : // will be used, so catch that as an error
271 0 : libmesh_error_msg_if(nonlinear_solver->jacobian && nonlinear_solver->jacobian_object,
272 : "ERROR: cannot specify both a function and object to compute the Jacobian!");
273 :
274 0 : libmesh_error_msg_if(nonlinear_solver->residual && nonlinear_solver->residual_object,
275 : "ERROR: cannot specify both a function and object to compute the Residual!");
276 :
277 0 : libmesh_error_msg_if(nonlinear_solver->matvec && nonlinear_solver->residual_and_jacobian_object,
278 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
279 :
280 :
281 0 : if (get_jacobian)
282 : {
283 0 : if (nonlinear_solver->jacobian != nullptr)
284 0 : nonlinear_solver->jacobian (*current_local_solution.get(), *matrix, *this);
285 :
286 0 : else if (nonlinear_solver->jacobian_object != nullptr)
287 0 : nonlinear_solver->jacobian_object->jacobian (*current_local_solution.get(), *matrix, *this);
288 :
289 0 : else if (nonlinear_solver->matvec != nullptr)
290 0 : nonlinear_solver->matvec (*current_local_solution.get(), get_residual ? rhs : nullptr, matrix, *this);
291 :
292 0 : else if (nonlinear_solver->residual_and_jacobian_object != nullptr)
293 0 : nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), get_residual ? rhs : nullptr, matrix, *this);
294 :
295 : else
296 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
297 : }
298 :
299 0 : if (get_residual)
300 : {
301 0 : if (nonlinear_solver->residual != nullptr)
302 0 : nonlinear_solver->residual (*current_local_solution.get(), *rhs, *this);
303 :
304 0 : else if (nonlinear_solver->residual_object != nullptr)
305 0 : nonlinear_solver->residual_object->residual (*current_local_solution.get(), *rhs, *this);
306 :
307 0 : else if (nonlinear_solver->matvec != nullptr)
308 : {
309 : // we might have already grabbed the residual and jacobian together
310 0 : if (!get_jacobian)
311 0 : nonlinear_solver->matvec (*current_local_solution.get(), rhs, nullptr, *this);
312 : }
313 :
314 0 : else if (nonlinear_solver->residual_and_jacobian_object != nullptr)
315 : {
316 : // we might have already grabbed the residual and jacobian together
317 0 : if (!get_jacobian)
318 0 : nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian (*current_local_solution.get(), rhs, nullptr, *this);
319 : }
320 :
321 : else
322 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
323 : }
324 : else
325 0 : libmesh_assert(get_jacobian); // I can't believe you really wanted to assemble *nothing*
326 0 : }
327 :
328 :
329 :
330 :
331 0 : unsigned NonlinearImplicitSystem::get_current_nonlinear_iteration_number() const
332 : {
333 0 : return nonlinear_solver->get_current_nonlinear_iteration_number();
334 : }
335 :
336 :
337 :
338 : } // namespace libMesh
|
