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 : // libMesh includes
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/equation_systems.h"
23 : #include "libmesh/fe_base.h"
24 : #include "libmesh/fem_context.h"
25 : #include "libmesh/fem_system.h"
26 : #include "libmesh/libmesh_logging.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/numeric_vector.h"
29 : #include "libmesh/parallel_algebra.h"
30 : #include "libmesh/parallel_ghost_sync.h"
31 : #include "libmesh/quadrature.h"
32 : #include "libmesh/sparse_matrix.h"
33 : #include "libmesh/time_solver.h"
34 : #include "libmesh/unsteady_solver.h" // For eulerian_residual
35 : #include "libmesh/fe_interface.h"
36 :
37 : namespace {
38 : using namespace libMesh;
39 :
40 : // give this guy some scope since there
41 : // is underlying vector allocation upon
42 : // creation/deletion
43 : ConstElemRange elem_range;
44 :
45 : typedef Threads::spin_mutex femsystem_mutex;
46 : femsystem_mutex assembly_mutex;
47 :
48 21041685 : void assemble_unconstrained_element_system(const FEMSystem & _sys,
49 : const bool _get_jacobian,
50 : const bool _constrain_heterogeneously,
51 : FEMContext & _femcontext)
52 : {
53 21041685 : if (_sys.print_element_solutions)
54 : {
55 0 : std::streamsize old_precision = libMesh::out.precision();
56 0 : libMesh::out.precision(16);
57 0 : if (_femcontext.has_elem())
58 0 : libMesh::out << "U_elem " << _femcontext.get_elem().id();
59 : else
60 0 : libMesh::out << "U_scalar ";
61 0 : libMesh::out << " = " << _femcontext.get_elem_solution() << std::endl;
62 :
63 0 : if (_sys.use_fixed_solution)
64 : {
65 0 : if (_femcontext.has_elem())
66 0 : libMesh::out << "Ufixed_elem " << _femcontext.get_elem().id();
67 : else
68 0 : libMesh::out << "Ufixed_scalar ";
69 0 : libMesh::out << " = " << _femcontext.get_elem_fixed_solution() << std::endl;
70 0 : libMesh::out.precision(old_precision);
71 : }
72 : }
73 :
74 : // We need jacobians to do heterogeneous residual constraints
75 21041685 : const bool need_jacobian =
76 1885846 : (_get_jacobian || _constrain_heterogeneously);
77 :
78 : bool jacobian_computed =
79 21041685 : _sys.time_solver->element_residual(need_jacobian, _femcontext);
80 :
81 : // Compute a numeric jacobian if we have to
82 21041685 : if (need_jacobian && !jacobian_computed)
83 : {
84 : // Make sure we didn't compute a jacobian and lie about it
85 21556 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
86 : // Logging of numerical jacobians is done separately
87 227292 : _sys.numerical_elem_jacobian(_femcontext);
88 : }
89 :
90 : // Compute a numeric jacobian if we're asked to verify the
91 : // analytic jacobian we got
92 21041685 : if (need_jacobian && jacobian_computed &&
93 6144148 : _sys.verify_analytic_jacobians != 0.0)
94 : {
95 0 : DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
96 :
97 0 : _femcontext.get_elem_jacobian().zero();
98 : // Logging of numerical jacobians is done separately
99 0 : _sys.numerical_elem_jacobian(_femcontext);
100 :
101 0 : Real analytic_norm = analytic_jacobian.l1_norm();
102 0 : Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
103 :
104 : // If we can continue, we'll probably prefer the analytic jacobian
105 0 : analytic_jacobian.swap(_femcontext.get_elem_jacobian());
106 :
107 : // The matrix "analytic_jacobian" will now hold the error matrix
108 0 : analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
109 0 : Real error_norm = analytic_jacobian.l1_norm();
110 :
111 0 : Real relative_error = error_norm /
112 0 : std::max(analytic_norm, numerical_norm);
113 :
114 0 : if (relative_error > _sys.verify_analytic_jacobians)
115 : {
116 0 : libMesh::err << "Relative error " << relative_error
117 0 : << " detected in analytic jacobian on element "
118 0 : << _femcontext.get_elem().id() << '!' << std::endl;
119 :
120 0 : std::streamsize old_precision = libMesh::out.precision();
121 0 : libMesh::out.precision(16);
122 0 : libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
123 0 : << _femcontext.get_elem_jacobian() << std::endl;
124 0 : analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
125 0 : libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
126 0 : << analytic_jacobian << std::endl;
127 :
128 0 : libMesh::out.precision(old_precision);
129 :
130 0 : libmesh_error_msg("Relative error too large, exiting!");
131 : }
132 0 : }
133 :
134 21041685 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
135 105634793 : for (_femcontext.side = 0; _femcontext.side != n_sides;
136 84593108 : ++_femcontext.side)
137 : {
138 : // Don't compute on non-boundary sides unless requested
139 169186216 : if (!_sys.get_physics()->compute_internal_sides &&
140 84593108 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr)
141 82111062 : continue;
142 :
143 : // Any mesh movement has already been done (and restored,
144 : // if the TimeSolver isn't broken), but
145 : // reinitializing the side FE objects is still necessary
146 2482046 : _femcontext.side_fe_reinit();
147 :
148 2917658 : DenseMatrix<Number> old_jacobian;
149 : // If we're in DEBUG mode, we should always verify that the
150 : // user's side_residual function doesn't alter our existing
151 : // jacobian and then lie about it
152 : #ifndef DEBUG
153 : // Even if we're not in DEBUG mode, when we're verifying
154 : // analytic jacobians we'll want to verify each side's
155 : // jacobian contribution separately.
156 2264240 : if (_sys.verify_analytic_jacobians != 0.0 && need_jacobian)
157 : #endif // ifndef DEBUG
158 : {
159 217806 : old_jacobian = _femcontext.get_elem_jacobian();
160 217806 : _femcontext.get_elem_jacobian().zero();
161 : }
162 :
163 : jacobian_computed =
164 2482046 : _sys.time_solver->side_residual(need_jacobian, _femcontext);
165 :
166 : // Compute a numeric jacobian if we have to
167 2482046 : if (need_jacobian && !jacobian_computed)
168 : {
169 : // If we have already backed up old_jacobian,
170 : // we can make sure side_residual didn't compute a
171 : // jacobian and lie about it.
172 : //
173 : // If we haven't, then we need to, to let
174 : // numerical_side_jacobian work.
175 49166 : if (old_jacobian.m())
176 5620 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
177 : else
178 : {
179 43546 : old_jacobian = _femcontext.get_elem_jacobian();
180 0 : _femcontext.get_elem_jacobian().zero();
181 : }
182 :
183 : // Logging of numerical jacobians is done separately
184 49166 : _sys.numerical_side_jacobian(_femcontext);
185 :
186 : // Add back in element interior numerical Jacobian
187 43546 : _femcontext.get_elem_jacobian() += old_jacobian;
188 : }
189 :
190 : // Compute a numeric jacobian if we're asked to verify the
191 : // analytic jacobian we got
192 2432880 : else if (need_jacobian && jacobian_computed &&
193 797997 : _sys.verify_analytic_jacobians != 0.0)
194 : {
195 0 : DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
196 :
197 0 : _femcontext.get_elem_jacobian().zero();
198 : // Logging of numerical jacobians is done separately
199 0 : _sys.numerical_side_jacobian(_femcontext);
200 :
201 0 : Real analytic_norm = analytic_jacobian.l1_norm();
202 0 : Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
203 :
204 : // If we can continue, we'll probably prefer the analytic jacobian
205 0 : analytic_jacobian.swap(_femcontext.get_elem_jacobian());
206 :
207 : // The matrix "analytic_jacobian" will now hold the error matrix
208 0 : analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
209 0 : Real error_norm = analytic_jacobian.l1_norm();
210 :
211 0 : Real relative_error = error_norm /
212 0 : std::max(analytic_norm, numerical_norm);
213 :
214 0 : if (relative_error > _sys.verify_analytic_jacobians)
215 : {
216 0 : libMesh::err << "Relative error " << relative_error
217 0 : << " detected in analytic jacobian on element "
218 0 : << _femcontext.get_elem().id()
219 0 : << ", side "
220 0 : << static_cast<unsigned int>(_femcontext.side) << '!' << std::endl;
221 :
222 0 : std::streamsize old_precision = libMesh::out.precision();
223 0 : libMesh::out.precision(16);
224 0 : libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
225 0 : << _femcontext.get_elem_jacobian() << std::endl;
226 0 : analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
227 0 : libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
228 0 : << analytic_jacobian << std::endl;
229 0 : libMesh::out.precision(old_precision);
230 :
231 0 : libmesh_error_msg("Relative error too large, exiting!");
232 : }
233 : // Once we've verified a side, we'll want to add back the
234 : // rest of the accumulated jacobian
235 0 : _femcontext.get_elem_jacobian() += old_jacobian;
236 0 : }
237 :
238 : // In DEBUG mode, we've set elem_jacobian == 0, and we
239 : // may have yet to add the old jacobian back
240 : #ifdef DEBUG
241 : else
242 : {
243 212186 : _femcontext.get_elem_jacobian() += old_jacobian;
244 : }
245 : #endif // ifdef DEBUG
246 2046434 : }
247 21041685 : }
248 :
249 21041685 : void add_element_system(FEMSystem & _sys,
250 : const bool _get_residual,
251 : const bool _get_jacobian,
252 : const bool _constrain_heterogeneously,
253 : const bool _no_constraints,
254 : FEMContext & _femcontext)
255 : {
256 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
257 21041685 : if (_get_residual && _sys.print_element_residuals)
258 : {
259 0 : std::streamsize old_precision = libMesh::out.precision();
260 0 : libMesh::out.precision(16);
261 0 : if (_femcontext.has_elem())
262 0 : libMesh::out << "Rraw_elem " << _femcontext.get_elem().id();
263 : else
264 0 : libMesh::out << "Rraw_scalar ";
265 0 : libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
266 0 : libMesh::out.precision(old_precision);
267 : }
268 :
269 21041685 : if (_get_jacobian && _sys.print_element_jacobians)
270 : {
271 0 : std::streamsize old_precision = libMesh::out.precision();
272 0 : libMesh::out.precision(16);
273 0 : if (_femcontext.has_elem())
274 0 : libMesh::out << "Jraw_elem " << _femcontext.get_elem().id();
275 : else
276 0 : libMesh::out << "Jraw_scalar ";
277 0 : libMesh::out << " = " << _femcontext.get_elem_jacobian() << std::endl;
278 0 : libMesh::out.precision(old_precision);
279 : }
280 :
281 : // We turn off the asymmetric constraint application iff we expect
282 : // enforce_constraints_exactly() to be called in the solver
283 21041685 : const bool constrain_in_solver = _sys.get_constrain_in_solver();
284 :
285 21041685 : if (_get_residual && _get_jacobian)
286 : {
287 4193579 : if (constrain_in_solver)
288 : {
289 3846417 : if (_constrain_heterogeneously)
290 0 : _sys.get_dof_map().heterogenously_constrain_element_matrix_and_vector
291 0 : (_femcontext.get_elem_jacobian(),
292 : _femcontext.get_elem_residual(),
293 : _femcontext.get_dof_indices(), false);
294 3846417 : else if (!_no_constraints)
295 345242 : _sys.get_dof_map().constrain_element_matrix_and_vector
296 3846417 : (_femcontext.get_elem_jacobian(),
297 : _femcontext.get_elem_residual(),
298 : _femcontext.get_dof_indices(), false);
299 : }
300 1760 : else if (!_no_constraints)
301 160 : _sys.get_dof_map().heterogeneously_constrain_element_jacobian_and_residual
302 1760 : (_femcontext.get_elem_jacobian(),
303 : _femcontext.get_elem_residual(),
304 : _femcontext.get_dof_indices(),
305 160 : *_sys.current_local_solution);
306 : // Do nothing if (_no_constraints)
307 : }
308 17193508 : else if (_get_residual)
309 : {
310 14672005 : if (constrain_in_solver)
311 : {
312 14662589 : if (_constrain_heterogeneously)
313 0 : _sys.get_dof_map().heterogenously_constrain_element_vector
314 0 : (_femcontext.get_elem_jacobian(),
315 : _femcontext.get_elem_residual(),
316 : _femcontext.get_dof_indices(), false);
317 14662589 : else if (!_no_constraints)
318 1299222 : _sys.get_dof_map().constrain_element_vector
319 14522693 : (_femcontext.get_elem_residual(),
320 : _femcontext.get_dof_indices(), false);
321 : }
322 9416 : else if (!_no_constraints)
323 160 : _sys.get_dof_map().heterogeneously_constrain_element_residual
324 1760 : (_femcontext.get_elem_residual(),
325 : _femcontext.get_dof_indices(),
326 160 : *_sys.current_local_solution);
327 : // Do nothing if (_no_constraints)
328 : }
329 2521503 : else if (_get_jacobian)
330 : {
331 : // Heterogeneous and homogeneous constraints are the same on the
332 : // matrix
333 : // Only get these contribs if we are applying some constraints
334 2521503 : if (!_no_constraints)
335 2750037 : _sys.get_dof_map().constrain_element_matrix (_femcontext.get_elem_jacobian(),
336 : _femcontext.get_dof_indices(),
337 2521503 : !constrain_in_solver);
338 : }
339 : #else
340 : libmesh_ignore(_constrain_heterogeneously, _no_constraints);
341 : #endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS
342 :
343 21041685 : if (_get_residual && _sys.print_element_residuals)
344 : {
345 0 : std::streamsize old_precision = libMesh::out.precision();
346 0 : libMesh::out.precision(16);
347 0 : if (_femcontext.has_elem())
348 0 : libMesh::out << "R_elem " << _femcontext.get_elem().id();
349 : else
350 0 : libMesh::out << "R_scalar ";
351 0 : libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
352 0 : libMesh::out.precision(old_precision);
353 : }
354 :
355 21041685 : if (_get_jacobian && _sys.print_element_jacobians)
356 : {
357 0 : std::streamsize old_precision = libMesh::out.precision();
358 0 : libMesh::out.precision(16);
359 0 : if (_femcontext.has_elem())
360 0 : libMesh::out << "J_elem " << _femcontext.get_elem().id();
361 : else
362 0 : libMesh::out << "J_scalar ";
363 0 : libMesh::out << " = " << _femcontext.get_elem_jacobian() << std::endl;
364 0 : libMesh::out.precision(old_precision);
365 : }
366 :
367 : { // A lock is necessary around access to the global system
368 3771692 : femsystem_mutex::scoped_lock lock(assembly_mutex);
369 :
370 21041685 : if (_get_jacobian)
371 6943616 : _sys.get_system_matrix().add_matrix (_femcontext.get_elem_jacobian(),
372 1147872 : _femcontext.get_dof_indices());
373 21041685 : if (_get_residual)
374 18520182 : _sys.rhs->add_vector (_femcontext.get_elem_residual(),
375 1657312 : _femcontext.get_dof_indices());
376 : } // Scope for assembly mutex
377 21041685 : }
378 :
379 :
380 :
381 : class AssemblyContributions
382 : {
383 : public:
384 : /**
385 : * constructor to set context
386 : */
387 5900 : AssemblyContributions(FEMSystem & sys,
388 : bool get_residual,
389 : bool get_jacobian,
390 : bool constrain_heterogeneously,
391 204784 : bool no_constraints) :
392 192984 : _sys(sys),
393 192984 : _get_residual(get_residual),
394 192984 : _get_jacobian(get_jacobian),
395 192984 : _constrain_heterogeneously(constrain_heterogeneously),
396 204784 : _no_constraints(no_constraints) {}
397 :
398 : /**
399 : * operator() for use with Threads::parallel_for().
400 : */
401 220030 : void operator()(const ConstElemRange & range) const
402 : {
403 231042 : std::unique_ptr<DiffContext> con = _sys.build_context();
404 11012 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
405 220030 : _sys.init_context(_femcontext);
406 :
407 21261715 : for (const auto & elem : range)
408 : {
409 21041685 : _femcontext.pre_fe_reinit(_sys, elem);
410 21041685 : _femcontext.elem_fe_reinit();
411 :
412 : assemble_unconstrained_element_system
413 21041685 : (_sys, _get_jacobian, _constrain_heterogeneously, _femcontext);
414 :
415 : add_element_system
416 24813377 : (_sys, _get_residual, _get_jacobian,
417 21041685 : _constrain_heterogeneously, _no_constraints, _femcontext);
418 : }
419 220030 : }
420 :
421 : private:
422 :
423 : FEMSystem & _sys;
424 :
425 : const bool _get_residual, _get_jacobian, _constrain_heterogeneously, _no_constraints;
426 : };
427 :
428 : class PostprocessContributions
429 : {
430 : public:
431 : /**
432 : * constructor to set context
433 : */
434 : explicit
435 1976 : PostprocessContributions(FEMSystem & sys) : _sys(sys) {}
436 :
437 : /**
438 : * operator() for use with Threads::parallel_for().
439 : */
440 2048 : void operator()(const ConstElemRange & range) const
441 : {
442 2252 : std::unique_ptr<DiffContext> con = _sys.build_context();
443 204 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
444 2048 : _sys.init_context(_femcontext);
445 :
446 64788 : for (const auto & elem : range)
447 : {
448 62740 : _femcontext.pre_fe_reinit(_sys, elem);
449 :
450 : // Optionally initialize all the interior FE objects on elem.
451 62740 : if (_sys.fe_reinit_during_postprocess)
452 62740 : _femcontext.elem_fe_reinit();
453 :
454 62740 : _sys.element_postprocess(_femcontext);
455 :
456 62740 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
457 310884 : for (_femcontext.side = 0; _femcontext.side != n_sides;
458 248144 : ++_femcontext.side)
459 : {
460 : // Don't compute on non-boundary sides unless requested
461 269216 : if (!_sys.postprocess_sides ||
462 316992 : (!_sys.get_physics()->compute_internal_sides &&
463 169032 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
464 232531 : continue;
465 :
466 : // Optionally initialize all the FE objects on this side.
467 15613 : if (_sys.fe_reinit_during_postprocess)
468 15613 : _femcontext.side_fe_reinit();
469 :
470 15613 : _sys.side_postprocess(_femcontext);
471 : }
472 : }
473 2048 : }
474 :
475 : private:
476 :
477 : FEMSystem & _sys;
478 : };
479 :
480 : class QoIContributions
481 : {
482 : public:
483 : /**
484 : * constructor to set context
485 : */
486 : explicit
487 39155 : QoIContributions(FEMSystem & sys,
488 : DifferentiableQoI & diff_qoi,
489 39155 : const QoISet & qoi_indices) :
490 40273 : qoi(sys.n_qois(), 0.), _sys(sys), _diff_qoi(diff_qoi),_qoi_indices(qoi_indices) {}
491 :
492 : /**
493 : * splitting constructor
494 : */
495 2236 : QoIContributions(const QoIContributions & other,
496 2236 : Threads::split) :
497 2236 : qoi(other._sys.n_qois(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {}
498 :
499 : /**
500 : * operator() for use with Threads::parallel_reduce().
501 : */
502 42509 : void operator()(const ConstElemRange & range)
503 : {
504 44745 : std::unique_ptr<DiffContext> con = _sys.build_context();
505 2236 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
506 42509 : _diff_qoi.init_context(_femcontext);
507 :
508 2236 : bool have_some_heterogenous_qoi_bc = false;
509 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
510 46981 : std::vector<bool> have_heterogenous_qoi_bc(_sys.n_qois(), false);
511 103258 : for (auto q : make_range(_sys.n_qois()))
512 63945 : if (_qoi_indices.has_index(q) &&
513 6392 : _sys.get_dof_map().has_heterogenous_adjoint_constraints(q))
514 : {
515 0 : have_heterogenous_qoi_bc[q] = true;
516 0 : have_some_heterogenous_qoi_bc = true;
517 : }
518 : #endif
519 :
520 42509 : if (have_some_heterogenous_qoi_bc)
521 0 : _sys.init_context(_femcontext);
522 :
523 13461113 : for (const auto & elem : range)
524 : {
525 13418604 : _femcontext.pre_fe_reinit(_sys, elem);
526 :
527 : // We might have some heterogenous dofs here; let's see for
528 : // certain
529 1216614 : bool elem_has_some_heterogenous_qoi_bc = false;
530 :
531 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
532 : const unsigned int n_dofs =
533 2433228 : cast_int<unsigned int>(_femcontext.get_dof_indices().size());
534 :
535 15851832 : std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.n_qois(), false);
536 13418604 : if (have_some_heterogenous_qoi_bc)
537 : {
538 0 : for (auto q : make_range(_sys.n_qois()))
539 : {
540 0 : if (have_heterogenous_qoi_bc[q])
541 : {
542 0 : for (auto d : make_range(n_dofs))
543 0 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
544 0 : (q, _femcontext.get_dof_indices()[d]) != Number(0))
545 : {
546 0 : elem_has_some_heterogenous_qoi_bc = true;
547 0 : elem_has_heterogenous_qoi_bc[q] = true;
548 0 : break;
549 : }
550 : }
551 : }
552 : }
553 : #endif
554 :
555 13418604 : if (_diff_qoi.assemble_qoi_elements ||
556 : elem_has_some_heterogenous_qoi_bc)
557 13418604 : _femcontext.elem_fe_reinit();
558 :
559 13418604 : if (_diff_qoi.assemble_qoi_elements)
560 13418604 : _diff_qoi.element_qoi(_femcontext, _qoi_indices);
561 :
562 : // If we have some heterogenous dofs here, those are
563 : // themselves part of a regularized flux QoI which the library
564 : // handles integrating
565 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
566 13418604 : if (elem_has_some_heterogenous_qoi_bc)
567 : {
568 0 : _sys.time_solver->element_residual(false, _femcontext);
569 :
570 0 : for (auto q : make_range(_sys.n_qois()))
571 : {
572 0 : if (elem_has_heterogenous_qoi_bc[q])
573 : {
574 0 : for (auto d : make_range(n_dofs))
575 0 : this->qoi[q] -= _femcontext.get_elem_residual()(d) *
576 0 : _sys.get_dof_map().has_heterogenous_adjoint_constraint(q, _femcontext.get_dof_indices()[d]);
577 :
578 : }
579 : }
580 : }
581 : #endif
582 :
583 13418604 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
584 67093020 : for (_femcontext.side = 0; _femcontext.side != n_sides;
585 53674416 : ++_femcontext.side)
586 : {
587 : // Don't compute on non-boundary sides unless requested
588 53675400 : if (!_diff_qoi.assemble_qoi_sides ||
589 11184 : (!_diff_qoi.assemble_qoi_internal_sides &&
590 11184 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
591 48805234 : continue;
592 :
593 2999 : _femcontext.side_fe_reinit();
594 :
595 2999 : _diff_qoi.side_qoi(_femcontext, _qoi_indices);
596 : }
597 : }
598 :
599 44745 : this->_diff_qoi.thread_join( this->qoi, _femcontext.get_qois(), _qoi_indices );
600 42509 : }
601 :
602 2236 : void join (const QoIContributions & other)
603 : {
604 1118 : libmesh_assert_equal_to (this->qoi.size(), other.qoi.size());
605 3354 : this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices );
606 2236 : }
607 :
608 : std::vector<Number> qoi;
609 :
610 : private:
611 :
612 : FEMSystem & _sys;
613 : DifferentiableQoI & _diff_qoi;
614 :
615 : const QoISet _qoi_indices;
616 : };
617 :
618 : class QoIDerivativeContributions
619 : {
620 : public:
621 : /**
622 : * constructor to set context
623 : */
624 464 : QoIDerivativeContributions(FEMSystem & sys,
625 : const QoISet & qoi_indices,
626 : DifferentiableQoI & qoi,
627 : bool include_liftfunc,
628 14654 : bool apply_constraints) :
629 13726 : _sys(sys),
630 13726 : _qoi_indices(qoi_indices),
631 13726 : _qoi(qoi),
632 13726 : _include_liftfunc(include_liftfunc),
633 14654 : _apply_constraints(apply_constraints) {}
634 :
635 : /**
636 : * operator() for use with Threads::parallel_for().
637 : */
638 16016 : void operator()(const ConstElemRange & range) const
639 : {
640 16944 : std::unique_ptr<DiffContext> con = _sys.build_context();
641 928 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
642 16016 : _qoi.init_context(_femcontext);
643 :
644 928 : bool have_some_heterogenous_qoi_bc = false;
645 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
646 17872 : std::vector<bool> have_heterogenous_qoi_bc(_sys.n_qois(), false);
647 16016 : if (_include_liftfunc || _apply_constraints)
648 43173 : for (auto q : make_range(_sys.n_qois()))
649 28757 : if (_qoi_indices.has_index(q) &&
650 3200 : _sys.get_dof_map().has_heterogenous_adjoint_constraints(q))
651 : {
652 48 : have_heterogenous_qoi_bc[q] = true;
653 24 : have_some_heterogenous_qoi_bc = true;
654 : }
655 : #endif
656 :
657 16016 : if (have_some_heterogenous_qoi_bc)
658 212 : _sys.init_context(_femcontext);
659 :
660 2432151 : for (const auto & elem : range)
661 : {
662 2416135 : _femcontext.pre_fe_reinit(_sys, elem);
663 :
664 : // We might have some heterogenous dofs here; let's see for
665 : // certain
666 219086 : bool elem_has_some_heterogenous_qoi_bc = false;
667 :
668 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
669 : const unsigned int n_dofs =
670 438172 : cast_int<unsigned int>(_femcontext.get_dof_indices().size());
671 :
672 2854307 : std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.n_qois(), false);
673 2416135 : if (have_some_heterogenous_qoi_bc)
674 : {
675 85476 : for (auto q : make_range(_sys.n_qois()))
676 : {
677 48450 : if (have_heterogenous_qoi_bc[q])
678 : {
679 412332 : for (auto d : make_range(n_dofs))
680 320980 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
681 421568 : (q, _femcontext.get_dof_indices()[d]) != Number(0))
682 : {
683 140 : elem_has_some_heterogenous_qoi_bc = true;
684 280 : elem_has_heterogenous_qoi_bc[q] = true;
685 1680 : break;
686 : }
687 : }
688 : }
689 : }
690 : #endif
691 :
692 : // If we're going to call a user integral, then we need FE
693 : // information to call element_qoi.
694 : // If we're going to evaluate lift-function-based components
695 : // of a QoI, then we need FE information to assemble the
696 : // element residual.
697 2416135 : if (_qoi.assemble_qoi_elements ||
698 0 : ((_include_liftfunc || _apply_constraints) &&
699 : elem_has_some_heterogenous_qoi_bc))
700 2416135 : _femcontext.elem_fe_reinit();
701 :
702 2416135 : if (_qoi.assemble_qoi_elements)
703 2416135 : _qoi.element_qoi_derivative(_femcontext, _qoi_indices);
704 :
705 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
706 : // If we need to use heterogenous dofs here, we need the
707 : // Jacobian either for the regularized flux QoI integration
708 : // and/or for constraint application.
709 2416135 : if ((_include_liftfunc || _apply_constraints) &&
710 : elem_has_some_heterogenous_qoi_bc)
711 : {
712 1820 : bool jacobian_computed = _sys.time_solver->element_residual(true, _femcontext);
713 :
714 : // If we're using numerical jacobians, above wont compute them
715 1680 : if (!jacobian_computed)
716 : {
717 : // Make sure we didn't compute a jacobian and lie about it
718 140 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
719 : // Logging of numerical jacobians is done separately
720 1680 : _sys.numerical_elem_jacobian(_femcontext);
721 : }
722 : }
723 :
724 : // If we have some heterogenous dofs here, those are
725 : // themselves part of a regularized flux QoI which the library
726 : // may handle integrating
727 2416135 : if (_include_liftfunc && elem_has_some_heterogenous_qoi_bc)
728 : {
729 0 : for (auto q : make_range(_sys.n_qois()))
730 : {
731 0 : if (elem_has_heterogenous_qoi_bc[q])
732 : {
733 0 : for (auto i : make_range(n_dofs))
734 : {
735 : Number liftfunc_val =
736 0 : _sys.get_dof_map().has_heterogenous_adjoint_constraint(q, _femcontext.get_dof_indices()[i]);
737 :
738 0 : if (liftfunc_val != Number(0))
739 : {
740 0 : for (auto j : make_range(n_dofs))
741 0 : _femcontext.get_qoi_derivatives()[q](j) -=
742 0 : _femcontext.get_elem_jacobian()(i,j) *
743 : liftfunc_val;
744 : }
745 : }
746 : }
747 : }
748 : }
749 : #endif
750 :
751 :
752 2416135 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
753 12080675 : for (_femcontext.side = 0; _femcontext.side != n_sides;
754 9664540 : ++_femcontext.side)
755 : {
756 : // Don't compute on non-boundary sides unless requested
757 9723924 : if (!_qoi.assemble_qoi_sides ||
758 649372 : (!_qoi.assemble_qoi_internal_sides &&
759 649372 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
760 8760594 : continue;
761 :
762 30960 : _femcontext.side_fe_reinit();
763 :
764 30960 : _qoi.side_qoi_derivative(_femcontext, _qoi_indices);
765 : }
766 :
767 : // We need some unmodified indices to use for constraining
768 : // multiple vector
769 : // FIXME - there should be a DofMap::constrain_element_vectors
770 : // to do this more efficiently
771 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
772 2635221 : std::vector<dof_id_type> original_dofs = _femcontext.get_dof_indices();
773 : #endif
774 :
775 : { // A lock is necessary around access to the global system
776 438172 : femsystem_mutex::scoped_lock lock(assembly_mutex);
777 :
778 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
779 : // We'll need to see if any heterogenous constraints apply
780 : // to the QoI dofs on this element *or* to any of the dofs
781 : // they depend on, so let's get those dependencies
782 2416135 : if (_apply_constraints)
783 2595285 : _sys.get_dof_map().constrain_nothing(_femcontext.get_dof_indices());
784 : #endif
785 :
786 4929428 : for (auto i : make_range(_sys.n_qois()))
787 2513293 : if (_qoi_indices.has_index(i))
788 : {
789 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
790 2513293 : if (_apply_constraints)
791 : {
792 : #ifndef NDEBUG
793 225458 : bool has_heterogenous_constraint = false;
794 1269226 : for (auto d : make_range(n_dofs))
795 1043908 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
796 1043908 : (i, _femcontext.get_dof_indices()[d]) != Number(0))
797 : {
798 140 : has_heterogenous_constraint = true;
799 140 : libmesh_assert(elem_has_heterogenous_qoi_bc[i]);
800 140 : libmesh_assert(elem_has_some_heterogenous_qoi_bc);
801 140 : break;
802 : }
803 : #else
804 : bool has_heterogenous_constraint =
805 225458 : elem_has_heterogenous_qoi_bc[i];
806 : #endif
807 :
808 2476429 : _femcontext.get_dof_indices() = original_dofs;
809 :
810 2476429 : if (has_heterogenous_constraint)
811 : {
812 : // Q_u gets used for *adjoint* solves, so we
813 : // need K^T here.
814 1960 : DenseMatrix<Number> elem_jacobian_transpose;
815 140 : _femcontext.get_elem_jacobian().get_transpose
816 1680 : (elem_jacobian_transpose);
817 :
818 1680 : _sys.get_dof_map().heterogenously_constrain_element_vector
819 1820 : (elem_jacobian_transpose,
820 280 : _femcontext.get_qoi_derivatives()[i],
821 : _femcontext.get_dof_indices(), false, i);
822 1400 : }
823 : else
824 : {
825 2474749 : _sys.get_dof_map().constrain_element_vector
826 2700067 : (_femcontext.get_qoi_derivatives()[i],
827 : _femcontext.get_dof_indices(), false);
828 : }
829 : }
830 : #endif
831 :
832 2513293 : _sys.get_adjoint_rhs(i).add_vector
833 2513293 : (_femcontext.get_qoi_derivatives()[i], _femcontext.get_dof_indices());
834 : }
835 : }
836 : }
837 16016 : }
838 :
839 : private:
840 :
841 : FEMSystem & _sys;
842 : const QoISet & _qoi_indices;
843 : DifferentiableQoI & _qoi;
844 : bool _include_liftfunc, _apply_constraints;
845 : };
846 :
847 :
848 : }
849 :
850 :
851 : namespace libMesh
852 : {
853 :
854 :
855 :
856 :
857 :
858 5044 : FEMSystem::FEMSystem (EquationSystems & es,
859 : const std::string & name_in,
860 0 : const unsigned int number_in)
861 : : Parent(es, name_in, number_in),
862 4744 : fe_reinit_during_postprocess(true),
863 4744 : numerical_jacobian_h(TOLERANCE),
864 5044 : verify_analytic_jacobians(0.0)
865 : {
866 5044 : }
867 :
868 :
869 4744 : FEMSystem::~FEMSystem () = default;
870 :
871 :
872 :
873 5044 : void FEMSystem::init_data ()
874 : {
875 : // First initialize LinearImplicitSystem data
876 5044 : Parent::init_data();
877 5044 : }
878 :
879 :
880 204784 : void FEMSystem::assembly (bool get_residual, bool get_jacobian,
881 : bool apply_heterogeneous_constraints,
882 : bool apply_no_constraints)
883 : {
884 5900 : libmesh_assert(get_residual || get_jacobian);
885 :
886 : // Log residual and jacobian and combined performance separately
887 : #ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
888 : const char * log_name;
889 11800 : if (get_residual && get_jacobian)
890 1218 : log_name = "assembly()";
891 9364 : else if (get_residual)
892 4030 : log_name = "assembly(get_residual)";
893 : else
894 652 : log_name = "assembly(get_jacobian)";
895 :
896 11800 : LOG_SCOPE(log_name, "FEMSystem");
897 : #endif
898 :
899 11800 : const MeshBase & mesh = this->get_mesh();
900 :
901 204784 : if (print_solution_norms)
902 : {
903 0 : this->solution->close();
904 :
905 0 : std::streamsize old_precision = libMesh::out.precision();
906 0 : libMesh::out.precision(16);
907 0 : libMesh::out << "|U| = "
908 0 : << this->solution->l1_norm()
909 0 : << std::endl;
910 0 : libMesh::out.precision(old_precision);
911 : }
912 204784 : if (print_solutions)
913 : {
914 0 : std::streamsize old_precision = libMesh::out.precision();
915 0 : libMesh::out.precision(16);
916 0 : libMesh::out << "U = [" << *(this->solution)
917 0 : << "];" << std::endl;
918 0 : libMesh::out.precision(old_precision);
919 : }
920 :
921 : // Is this definitely necessary? [RHS]
922 : // Yes. [RHS 2012]
923 204784 : if (get_jacobian)
924 63332 : matrix->zero();
925 204784 : if (get_residual)
926 183460 : rhs->zero();
927 :
928 : // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
929 204784 : if (verify_analytic_jacobians > 0.5)
930 : {
931 0 : libMesh::err << "WARNING! verify_analytic_jacobians was set "
932 0 : << "to absurdly large value of "
933 0 : << verify_analytic_jacobians << std::endl;
934 0 : libMesh::err << "Resetting to 1e-6!" << std::endl;
935 0 : verify_analytic_jacobians = 1e-6;
936 : }
937 :
938 : // In time-dependent problems, the nonlinear function we're trying
939 : // to solve at each timestep may depend on the particular solver
940 : // we're using
941 5900 : libmesh_assert(time_solver.get());
942 :
943 : // Build the residual and jacobian contributions on every active
944 : // mesh element on this processor
945 : Threads::parallel_for
946 415468 : (elem_range.reset(mesh.active_local_elements_begin(),
947 409568 : mesh.active_local_elements_end()),
948 403668 : AssemblyContributions(*this, get_residual, get_jacobian,
949 : apply_heterogeneous_constraints,
950 : apply_no_constraints));
951 :
952 : // Check and see if we have SCALAR variables
953 5900 : bool have_scalar = false;
954 453098 : for (auto i : make_range(this->n_variable_groups()))
955 : {
956 248314 : if (this->variable_group(i).type().family == SCALAR)
957 : {
958 0 : have_scalar = true;
959 0 : break;
960 : }
961 : }
962 :
963 : // SCALAR dofs are stored on the last processor, so we'll evaluate
964 : // their equation terms there and only if we have a SCALAR variable
965 210684 : if (this->processor_id() == (this->n_processors()-1) && have_scalar)
966 : {
967 0 : std::unique_ptr<DiffContext> con = this->build_context();
968 0 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
969 0 : this->init_context(_femcontext);
970 0 : _femcontext.pre_fe_reinit(*this, nullptr);
971 :
972 : bool jacobian_computed =
973 0 : this->time_solver->nonlocal_residual(get_jacobian, _femcontext);
974 :
975 : // Nonlocal residuals are likely to be length 0, in which case we
976 : // don't need to do any more. And we shouldn't try to do any
977 : // more; lots of DenseVector/DenseMatrix code assumes rank>0.
978 0 : if (_femcontext.get_elem_residual().size())
979 : {
980 : // Compute a numeric jacobian if we have to
981 0 : if (get_jacobian && !jacobian_computed)
982 : {
983 : // Make sure we didn't compute a jacobian and lie about it
984 0 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
985 : // Logging of numerical jacobians is done separately
986 0 : this->numerical_nonlocal_jacobian(_femcontext);
987 : }
988 :
989 : // Compute a numeric jacobian if we're asked to verify the
990 : // analytic jacobian we got
991 0 : if (get_jacobian && jacobian_computed &&
992 0 : this->verify_analytic_jacobians != 0.0)
993 : {
994 0 : DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
995 :
996 0 : _femcontext.get_elem_jacobian().zero();
997 : // Logging of numerical jacobians is done separately
998 0 : this->numerical_nonlocal_jacobian(_femcontext);
999 :
1000 0 : Real analytic_norm = analytic_jacobian.l1_norm();
1001 0 : Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
1002 :
1003 : // If we can continue, we'll probably prefer the analytic jacobian
1004 0 : analytic_jacobian.swap(_femcontext.get_elem_jacobian());
1005 :
1006 : // The matrix "analytic_jacobian" will now hold the error matrix
1007 0 : analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
1008 0 : Real error_norm = analytic_jacobian.l1_norm();
1009 :
1010 0 : Real relative_error = error_norm /
1011 0 : std::max(analytic_norm, numerical_norm);
1012 :
1013 0 : if (relative_error > this->verify_analytic_jacobians)
1014 : {
1015 0 : libMesh::err << "Relative error " << relative_error
1016 0 : << " detected in analytic jacobian on nonlocal dofs!"
1017 0 : << std::endl;
1018 :
1019 0 : std::streamsize old_precision = libMesh::out.precision();
1020 0 : libMesh::out.precision(16);
1021 0 : libMesh::out << "J_analytic nonlocal = "
1022 0 : << _femcontext.get_elem_jacobian() << std::endl;
1023 0 : analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
1024 0 : libMesh::out << "J_numeric nonlocal = "
1025 0 : << analytic_jacobian << std::endl;
1026 :
1027 0 : libMesh::out.precision(old_precision);
1028 :
1029 0 : libmesh_error_msg("Relative error too large, exiting!");
1030 : }
1031 0 : }
1032 :
1033 : add_element_system
1034 0 : (*this, get_residual, get_jacobian,
1035 : apply_heterogeneous_constraints, apply_no_constraints, _femcontext);
1036 : }
1037 0 : }
1038 :
1039 204784 : if (get_residual && (print_residual_norms || print_residuals))
1040 0 : this->rhs->close();
1041 184764 : if (get_residual && print_residual_norms)
1042 : {
1043 0 : std::streamsize old_precision = libMesh::out.precision();
1044 0 : libMesh::out.precision(16);
1045 0 : libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
1046 0 : libMesh::out.precision(old_precision);
1047 : }
1048 184764 : if (get_residual && print_residuals)
1049 : {
1050 0 : std::streamsize old_precision = libMesh::out.precision();
1051 0 : libMesh::out.precision(16);
1052 0 : libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
1053 0 : libMesh::out.precision(old_precision);
1054 : }
1055 :
1056 204784 : if (get_jacobian && (print_jacobian_norms || print_jacobians))
1057 0 : this->matrix->close();
1058 71392 : if (get_jacobian && print_jacobian_norms)
1059 : {
1060 0 : std::streamsize old_precision = libMesh::out.precision();
1061 0 : libMesh::out.precision(16);
1062 0 : libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
1063 0 : libMesh::out.precision(old_precision);
1064 : }
1065 71392 : if (get_jacobian && print_jacobians)
1066 : {
1067 0 : std::streamsize old_precision = libMesh::out.precision();
1068 0 : libMesh::out.precision(16);
1069 0 : libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
1070 0 : libMesh::out.precision(old_precision);
1071 : }
1072 204784 : }
1073 :
1074 :
1075 :
1076 25037 : void FEMSystem::solve()
1077 : {
1078 : // We are solving the primal problem
1079 25037 : Parent::solve();
1080 :
1081 : // On a moving mesh we want the mesh to reflect the new solution
1082 25037 : this->mesh_position_set();
1083 25037 : }
1084 :
1085 :
1086 :
1087 27312 : void FEMSystem::mesh_position_set()
1088 : {
1089 : // If we don't need to move the mesh, we're done
1090 27312 : if (_mesh_sys != this)
1091 24712 : return;
1092 :
1093 0 : MeshBase & mesh = this->get_mesh();
1094 :
1095 2600 : std::unique_ptr<DiffContext> con = this->build_context();
1096 0 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
1097 2600 : this->init_context(_femcontext);
1098 :
1099 : // Move every mesh element we can
1100 51280 : for (const auto & elem : mesh.active_local_element_ptr_range())
1101 : {
1102 : // We need the algebraic data
1103 23040 : _femcontext.pre_fe_reinit(*this, elem);
1104 : // And when asserts are on, we also need the FE so
1105 : // we can assert that the mesh data is of the right type.
1106 : #ifndef NDEBUG
1107 0 : _femcontext.elem_fe_reinit();
1108 : #endif
1109 :
1110 : // This code won't handle moving subactive elements
1111 0 : libmesh_assert(!_femcontext.get_elem().has_children());
1112 :
1113 0 : _femcontext.elem_position_set(1.);
1114 2600 : }
1115 :
1116 : // We've now got positions set on all local nodes (and some
1117 : // semilocal nodes); let's request positions for non-local nodes
1118 : // from their processors.
1119 :
1120 2600 : SyncNodalPositions sync_object(mesh);
1121 : Parallel::sync_dofobject_data_by_id
1122 5200 : (this->comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
1123 2600 : }
1124 :
1125 :
1126 :
1127 1976 : void FEMSystem::postprocess ()
1128 : {
1129 102 : LOG_SCOPE("postprocess()", "FEMSystem");
1130 :
1131 204 : const MeshBase & mesh = this->get_mesh();
1132 :
1133 1976 : this->update();
1134 :
1135 : // Get the time solver object associated with the system, and tell it that
1136 : // we are not solving the adjoint problem
1137 102 : this->get_time_solver().set_is_adjoint(false);
1138 :
1139 : // Loop over every active mesh element on this processor
1140 4054 : Threads::parallel_for (elem_range.reset(mesh.active_local_elements_begin(),
1141 2282 : mesh.active_local_elements_end()),
1142 2078 : PostprocessContributions(*this));
1143 1976 : }
1144 :
1145 :
1146 :
1147 39155 : void FEMSystem::assemble_qoi (const QoISet & qoi_indices)
1148 : {
1149 2236 : LOG_SCOPE("assemble_qoi()", "FEMSystem");
1150 :
1151 2236 : const MeshBase & mesh = this->get_mesh();
1152 :
1153 39155 : this->update();
1154 :
1155 1118 : const unsigned int Nq = this->n_qois();
1156 :
1157 : // the quantity of interest is assumed to be a sum of element and
1158 : // side terms
1159 95110 : for (unsigned int i=0; i != Nq; ++i)
1160 55955 : if (qoi_indices.has_index(i))
1161 55955 : this->set_qoi(i, 0);
1162 :
1163 : // Create a non-temporary qoi_contributions object, so we can query
1164 : // its results after the reduction
1165 79428 : QoIContributions qoi_contributions(*this, *(this->get_qoi()), qoi_indices);
1166 :
1167 : // Loop over every active mesh element on this processor
1168 79428 : Threads::parallel_reduce(elem_range.reset(mesh.active_local_elements_begin(),
1169 42509 : mesh.active_local_elements_end()),
1170 : qoi_contributions);
1171 :
1172 39155 : std::vector<Number> global_qoi = this->get_qoi_values();
1173 39155 : this->get_qoi()->parallel_op( this->comm(), global_qoi, qoi_contributions.qoi, qoi_indices );
1174 77192 : this->set_qoi(std::move(global_qoi));
1175 39155 : }
1176 :
1177 :
1178 :
1179 14654 : void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices,
1180 : bool include_liftfunc,
1181 : bool apply_constraints)
1182 : {
1183 928 : LOG_SCOPE("assemble_qoi_derivative()", "FEMSystem");
1184 :
1185 928 : const MeshBase & mesh = this->get_mesh();
1186 :
1187 14654 : this->update();
1188 :
1189 : // The quantity of interest derivative assembly accumulates on
1190 : // initially zero vectors
1191 39471 : for (auto i : make_range(this->n_qois()))
1192 24817 : if (qoi_indices.has_index(i))
1193 24817 : this->add_adjoint_rhs(i).zero();
1194 :
1195 : // Loop over every active mesh element on this processor
1196 29772 : Threads::parallel_for (elem_range.reset(mesh.active_local_elements_begin(),
1197 16046 : mesh.active_local_elements_end()),
1198 14190 : QoIDerivativeContributions(*this, qoi_indices,
1199 14654 : *(this->get_qoi()),
1200 : include_liftfunc,
1201 : apply_constraints));
1202 :
1203 39471 : for (auto i : make_range(this->n_qois()))
1204 24817 : if (qoi_indices.has_index(i))
1205 24817 : this->get_qoi()->finalize_derivative(this->get_adjoint_rhs(i),i);
1206 14654 : }
1207 :
1208 :
1209 :
1210 278138 : void FEMSystem::numerical_jacobian (TimeSolverResPtr res,
1211 : FEMContext & context) const
1212 : {
1213 : // Logging is done by numerical_elem_jacobian
1214 : // or numerical_side_jacobian
1215 :
1216 54632 : DenseVector<Number> original_residual(context.get_elem_residual());
1217 54632 : DenseVector<Number> backwards_residual(context.get_elem_residual());
1218 332770 : DenseMatrix<Number> numeric_jacobian(context.get_elem_jacobian());
1219 : #ifdef DEBUG
1220 54632 : DenseMatrix<Number> old_jacobian(context.get_elem_jacobian());
1221 : #endif
1222 :
1223 27316 : Real numerical_point_h = 0.;
1224 278138 : if (_mesh_sys == this)
1225 0 : numerical_point_h = numerical_jacobian_h * context.get_elem().hmin();
1226 :
1227 : const unsigned int n_dofs =
1228 54632 : cast_int<unsigned int>(context.get_dof_indices().size());
1229 :
1230 608548 : for (auto v : make_range(context.n_vars()))
1231 : {
1232 298738 : const Real my_h = this->numerical_jacobian_h_for_var(v);
1233 :
1234 31672 : unsigned int j_offset = libMesh::invalid_uint;
1235 :
1236 330410 : if (!context.get_dof_indices(v).empty())
1237 : {
1238 4320260 : for (auto i : make_range(n_dofs))
1239 4729306 : if (context.get_dof_indices()[i] ==
1240 369728 : context.get_dof_indices(v)[0])
1241 31672 : j_offset = i;
1242 :
1243 31672 : libmesh_assert_not_equal_to(j_offset, libMesh::invalid_uint);
1244 : }
1245 :
1246 3174308 : for (auto j : make_range(context.get_dof_indices(v).size()))
1247 : {
1248 2843898 : const unsigned int total_j = j + j_offset;
1249 :
1250 : // Take the "minus" side of a central differenced first derivative
1251 2843898 : Number original_solution = context.get_elem_solution(v)(j);
1252 2597836 : context.get_elem_solution(v)(j) -= my_h;
1253 :
1254 : // Make sure to catch any moving mesh terms
1255 274232 : Real * coord = nullptr;
1256 2843898 : if (_mesh_sys == this)
1257 : {
1258 0 : if (_mesh_x_var == v)
1259 0 : coord = &(context.get_elem().point(j)(0));
1260 0 : else if (_mesh_y_var == v)
1261 0 : coord = &(context.get_elem().point(j)(1));
1262 0 : else if (_mesh_z_var == v)
1263 0 : coord = &(context.get_elem().point(j)(2));
1264 : }
1265 274232 : if (coord)
1266 : {
1267 : // We have enough information to scale the perturbations
1268 : // here appropriately
1269 0 : context.get_elem_solution(v)(j) = original_solution - numerical_point_h;
1270 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1271 : }
1272 :
1273 274232 : context.get_elem_residual().zero();
1274 2843898 : ((*time_solver).*(res))(false, context);
1275 : #ifdef DEBUG
1276 274232 : libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1277 : #endif
1278 274232 : backwards_residual = context.get_elem_residual();
1279 :
1280 : // Take the "plus" side of a central differenced first derivative
1281 2843898 : context.get_elem_solution(v)(j) = original_solution + my_h;
1282 2843898 : if (coord)
1283 : {
1284 0 : context.get_elem_solution()(j) = original_solution + numerical_point_h;
1285 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1286 : }
1287 274232 : context.get_elem_residual().zero();
1288 2843898 : ((*time_solver).*(res))(false, context);
1289 : #ifdef DEBUG
1290 274232 : libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1291 : #endif
1292 :
1293 2843898 : context.get_elem_solution(v)(j) = original_solution;
1294 2843898 : if (coord)
1295 : {
1296 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1297 0 : for (auto i : make_range(n_dofs))
1298 : {
1299 0 : numeric_jacobian(i,total_j) =
1300 0 : (context.get_elem_residual()(i) - backwards_residual(i)) /
1301 0 : 2. / numerical_point_h;
1302 : }
1303 : }
1304 : else
1305 : {
1306 37023780 : for (auto i : make_range(n_dofs))
1307 : {
1308 34179882 : numeric_jacobian(i,total_j) =
1309 31249924 : (context.get_elem_residual()(i) - backwards_residual(i)) /
1310 31249924 : 2. / my_h;
1311 : }
1312 : }
1313 : }
1314 : }
1315 :
1316 27316 : context.get_elem_residual() = original_residual;
1317 278138 : context.get_elem_jacobian() = numeric_jacobian;
1318 501644 : }
1319 :
1320 :
1321 :
1322 228972 : void FEMSystem::numerical_elem_jacobian (FEMContext & context) const
1323 : {
1324 43392 : LOG_SCOPE("numerical_elem_jacobian()", "FEMSystem");
1325 228972 : this->numerical_jacobian(&TimeSolver::element_residual, context);
1326 228972 : }
1327 :
1328 :
1329 :
1330 49166 : void FEMSystem::numerical_side_jacobian (FEMContext & context) const
1331 : {
1332 11240 : LOG_SCOPE("numerical_side_jacobian()", "FEMSystem");
1333 49166 : this->numerical_jacobian(&TimeSolver::side_residual, context);
1334 49166 : }
1335 :
1336 :
1337 :
1338 0 : void FEMSystem::numerical_nonlocal_jacobian (FEMContext & context) const
1339 : {
1340 0 : LOG_SCOPE("numerical_nonlocal_jacobian()", "FEMSystem");
1341 0 : this->numerical_jacobian(&TimeSolver::nonlocal_residual, context);
1342 0 : }
1343 :
1344 :
1345 :
1346 284324 : std::unique_ptr<DiffContext> FEMSystem::build_context ()
1347 : {
1348 298728 : auto fc = std::make_unique<FEMContext>(*this);
1349 :
1350 284324 : DifferentiablePhysics * phys = this->get_physics();
1351 :
1352 14404 : libmesh_assert (phys);
1353 :
1354 : // If we are solving a moving mesh problem, tell that to the Context
1355 284324 : fc->set_mesh_system(phys->get_mesh_system());
1356 28808 : fc->set_mesh_x_var(phys->get_mesh_x_var());
1357 28808 : fc->set_mesh_y_var(phys->get_mesh_y_var());
1358 28808 : fc->set_mesh_z_var(phys->get_mesh_z_var());
1359 :
1360 284324 : fc->set_deltat_pointer( &deltat );
1361 :
1362 : // If we are solving the adjoint problem, tell that to the Context
1363 298728 : fc->is_adjoint() = this->get_time_solver().is_adjoint();
1364 :
1365 298728 : return fc;
1366 255516 : }
1367 :
1368 :
1369 :
1370 208445 : void FEMSystem::init_context(DiffContext & c)
1371 : {
1372 : // Parent::init_context(c); // may be a good idea in derived classes
1373 :
1374 : // Although we do this in DiffSystem::build_context() and
1375 : // FEMSystem::build_context() as well, we do it here just to be
1376 : // extra sure that the deltat pointer gets set. Since the
1377 : // intended behavior is for classes derived from FEMSystem to
1378 : // call Parent::init_context() in their own init_context()
1379 : // overloads, we can ensure that those classes get the correct
1380 : // deltat pointers even if they have different build_context()
1381 : // overloads.
1382 208445 : c.set_deltat_pointer ( &deltat );
1383 :
1384 10492 : FEMContext & context = cast_ref<FEMContext &>(c);
1385 :
1386 : // Make sure we're prepared to do mass integration
1387 440211 : for (auto var : make_range(this->n_vars()))
1388 451918 : if (this->get_physics()->is_time_evolving(var))
1389 : {
1390 : // Request shape functions based on FEType
1391 192107 : switch( FEInterface::field_type( this->variable_type(var) ) )
1392 : {
1393 192107 : case( TYPE_SCALAR ):
1394 : {
1395 9596 : FEBase * elem_fe = nullptr;
1396 385369 : for (auto dim : context.elem_dimensions())
1397 : {
1398 9656 : context.get_element_fe(var, elem_fe, dim);
1399 19638 : elem_fe->get_JxW();
1400 9656 : elem_fe->get_phi();
1401 0 : }
1402 : }
1403 9596 : break;
1404 0 : case( TYPE_VECTOR ):
1405 : {
1406 0 : FEGenericBase<RealGradient> * elem_fe = nullptr;
1407 0 : for (auto dim : context.elem_dimensions())
1408 : {
1409 0 : context.get_element_fe(var, elem_fe, dim);
1410 0 : elem_fe->get_JxW();
1411 0 : elem_fe->get_phi();
1412 0 : }
1413 : }
1414 0 : break;
1415 0 : default:
1416 0 : libmesh_error_msg("Unrecognized field type!");
1417 : }
1418 : }
1419 208445 : }
1420 :
1421 :
1422 :
1423 65 : void FEMSystem::mesh_position_get()
1424 : {
1425 : // This function makes no sense unless we've already picked out some
1426 : // variable(s) to reflect mesh position coordinates
1427 65 : libmesh_error_msg_if(!_mesh_sys, "_mesh_sys was nullptr!");
1428 :
1429 : // We currently assume mesh variables are in our own system
1430 65 : if (_mesh_sys != this)
1431 0 : libmesh_not_implemented();
1432 :
1433 : // Loop over every active mesh element on this processor
1434 0 : const MeshBase & mesh = this->get_mesh();
1435 :
1436 65 : std::unique_ptr<DiffContext> con = this->build_context();
1437 0 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
1438 65 : this->init_context(_femcontext);
1439 :
1440 : // Get the solution's mesh variables from every element
1441 1282 : for (const auto & elem : mesh.active_local_element_ptr_range())
1442 : {
1443 576 : _femcontext.pre_fe_reinit(*this, elem);
1444 :
1445 576 : _femcontext.elem_position_get();
1446 :
1447 576 : if (_mesh_x_var != libMesh::invalid_uint)
1448 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_x_var),
1449 0 : _femcontext.get_dof_indices(_mesh_x_var) );
1450 576 : if (_mesh_y_var != libMesh::invalid_uint)
1451 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_y_var),
1452 0 : _femcontext.get_dof_indices(_mesh_y_var));
1453 576 : if (_mesh_z_var != libMesh::invalid_uint)
1454 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_z_var),
1455 0 : _femcontext.get_dof_indices(_mesh_z_var));
1456 65 : }
1457 :
1458 65 : this->solution->close();
1459 :
1460 : // And make sure the current_local_solution is up to date too
1461 65 : this->System::update();
1462 65 : }
1463 :
1464 : } // namespace libMesh
|
