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 : // 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/mesh_tools.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/parallel_algebra.h"
31 : #include "libmesh/parallel_ghost_sync.h"
32 : #include "libmesh/quadrature.h"
33 : #include "libmesh/sparse_matrix.h"
34 : #include "libmesh/time_solver.h"
35 : #include "libmesh/unsteady_solver.h" // For eulerian_residual
36 : #include "libmesh/fe_interface.h"
37 :
38 : namespace {
39 : using namespace libMesh;
40 :
41 : typedef Threads::spin_mutex femsystem_mutex;
42 : femsystem_mutex assembly_mutex;
43 :
44 21265831 : void assemble_unconstrained_element_system(const FEMSystem & _sys,
45 : const bool _get_jacobian,
46 : const bool _constrain_heterogeneously,
47 : FEMContext & _femcontext)
48 : {
49 21265831 : if (_sys.print_element_solutions)
50 : {
51 0 : std::streamsize old_precision = libMesh::out.precision();
52 0 : libMesh::out.precision(16);
53 0 : if (_femcontext.has_elem())
54 0 : libMesh::out << "U_elem " << _femcontext.get_elem().id();
55 : else
56 0 : libMesh::out << "U_scalar ";
57 0 : libMesh::out << " = " << _femcontext.get_elem_solution() << std::endl;
58 :
59 0 : if (_sys.use_fixed_solution)
60 : {
61 0 : if (_femcontext.has_elem())
62 0 : libMesh::out << "Ufixed_elem " << _femcontext.get_elem().id();
63 : else
64 0 : libMesh::out << "Ufixed_scalar ";
65 0 : libMesh::out << " = " << _femcontext.get_elem_fixed_solution() << std::endl;
66 0 : libMesh::out.precision(old_precision);
67 : }
68 : }
69 :
70 : // We need jacobians to do heterogeneous residual constraints
71 21265831 : const bool need_jacobian =
72 1898390 : (_get_jacobian || _constrain_heterogeneously);
73 :
74 : bool jacobian_computed =
75 21265831 : _sys.time_solver->element_residual(need_jacobian, _femcontext);
76 :
77 : // Compute a numeric jacobian if we have to
78 21265831 : if (need_jacobian && !jacobian_computed)
79 : {
80 : // Make sure we didn't compute a jacobian and lie about it
81 21472 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
82 : // Logging of numerical jacobians is done separately
83 227064 : _sys.numerical_elem_jacobian(_femcontext);
84 : }
85 :
86 : // Compute a numeric jacobian if we're asked to verify the
87 : // analytic jacobian we got
88 21265831 : if (need_jacobian && jacobian_computed &&
89 6252220 : _sys.verify_analytic_jacobians != 0.0)
90 : {
91 0 : DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
92 :
93 0 : _femcontext.get_elem_jacobian().zero();
94 : // Logging of numerical jacobians is done separately
95 0 : _sys.numerical_elem_jacobian(_femcontext);
96 :
97 0 : Real analytic_norm = analytic_jacobian.l1_norm();
98 0 : Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
99 :
100 : // If we can continue, we'll probably prefer the analytic jacobian
101 0 : analytic_jacobian.swap(_femcontext.get_elem_jacobian());
102 :
103 : // The matrix "analytic_jacobian" will now hold the error matrix
104 0 : analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
105 0 : Real error_norm = analytic_jacobian.l1_norm();
106 :
107 0 : Real relative_error = error_norm /
108 0 : std::max(analytic_norm, numerical_norm);
109 :
110 0 : if (relative_error > _sys.verify_analytic_jacobians)
111 : {
112 0 : libMesh::err << "Relative error " << relative_error
113 0 : << " detected in analytic jacobian on element "
114 0 : << _femcontext.get_elem().id() << '!' << std::endl;
115 :
116 0 : std::streamsize old_precision = libMesh::out.precision();
117 0 : libMesh::out.precision(16);
118 0 : libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
119 0 : << _femcontext.get_elem_jacobian() << std::endl;
120 0 : analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
121 0 : libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
122 0 : << analytic_jacobian << std::endl;
123 :
124 0 : libMesh::out.precision(old_precision);
125 :
126 0 : libmesh_error_msg("Relative error too large, exiting!");
127 : }
128 0 : }
129 :
130 21265831 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
131 106813139 : for (_femcontext.side = 0; _femcontext.side != n_sides;
132 85547308 : ++_femcontext.side)
133 : {
134 : // Don't compute on non-boundary sides unless requested
135 171094616 : if (!_sys.get_physics()->compute_internal_sides &&
136 85547308 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr)
137 82998150 : continue;
138 :
139 : // Any mesh movement has already been done (and restored,
140 : // if the TimeSolver isn't broken), but
141 : // reinitializing the side FE objects is still necessary
142 2549158 : _femcontext.side_fe_reinit();
143 :
144 2993478 : DenseMatrix<Number> old_jacobian;
145 : // If we're in DEBUG mode, we should always verify that the
146 : // user's side_residual function doesn't alter our existing
147 : // jacobian and then lie about it
148 : #ifndef DEBUG
149 : // Even if we're not in DEBUG mode, when we're verifying
150 : // analytic jacobians we'll want to verify each side's
151 : // jacobian contribution separately.
152 2326998 : if (_sys.verify_analytic_jacobians != 0.0 && need_jacobian)
153 : #endif // ifndef DEBUG
154 : {
155 222160 : old_jacobian = _femcontext.get_elem_jacobian();
156 222160 : _femcontext.get_elem_jacobian().zero();
157 : }
158 :
159 : jacobian_computed =
160 2549158 : _sys.time_solver->side_residual(need_jacobian, _femcontext);
161 :
162 : // Compute a numeric jacobian if we have to
163 2549158 : if (need_jacobian && !jacobian_computed)
164 : {
165 : // If we have already backed up old_jacobian,
166 : // we can make sure side_residual didn't compute a
167 : // jacobian and lie about it.
168 : //
169 : // If we haven't, then we need to, to let
170 : // numerical_side_jacobian work.
171 49182 : if (old_jacobian.m())
172 5624 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
173 : else
174 : {
175 43558 : old_jacobian = _femcontext.get_elem_jacobian();
176 0 : _femcontext.get_elem_jacobian().zero();
177 : }
178 :
179 : // Logging of numerical jacobians is done separately
180 49182 : _sys.numerical_side_jacobian(_femcontext);
181 :
182 : // Add back in element interior numerical Jacobian
183 43558 : _femcontext.get_elem_jacobian() += old_jacobian;
184 : }
185 :
186 : // Compute a numeric jacobian if we're asked to verify the
187 : // analytic jacobian we got
188 2499976 : else if (need_jacobian && jacobian_computed &&
189 827981 : _sys.verify_analytic_jacobians != 0.0)
190 : {
191 0 : DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
192 :
193 0 : _femcontext.get_elem_jacobian().zero();
194 : // Logging of numerical jacobians is done separately
195 0 : _sys.numerical_side_jacobian(_femcontext);
196 :
197 0 : Real analytic_norm = analytic_jacobian.l1_norm();
198 0 : Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
199 :
200 : // If we can continue, we'll probably prefer the analytic jacobian
201 0 : analytic_jacobian.swap(_femcontext.get_elem_jacobian());
202 :
203 : // The matrix "analytic_jacobian" will now hold the error matrix
204 0 : analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
205 0 : Real error_norm = analytic_jacobian.l1_norm();
206 :
207 0 : Real relative_error = error_norm /
208 0 : std::max(analytic_norm, numerical_norm);
209 :
210 0 : if (relative_error > _sys.verify_analytic_jacobians)
211 : {
212 0 : libMesh::err << "Relative error " << relative_error
213 0 : << " detected in analytic jacobian on element "
214 0 : << _femcontext.get_elem().id()
215 0 : << ", side "
216 0 : << static_cast<unsigned int>(_femcontext.side) << '!' << std::endl;
217 :
218 0 : std::streamsize old_precision = libMesh::out.precision();
219 0 : libMesh::out.precision(16);
220 0 : libMesh::out << "J_analytic " << _femcontext.get_elem().id() << " = "
221 0 : << _femcontext.get_elem_jacobian() << std::endl;
222 0 : analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
223 0 : libMesh::out << "J_numeric " << _femcontext.get_elem().id() << " = "
224 0 : << analytic_jacobian << std::endl;
225 0 : libMesh::out.precision(old_precision);
226 :
227 0 : libmesh_error_msg("Relative error too large, exiting!");
228 : }
229 : // Once we've verified a side, we'll want to add back the
230 : // rest of the accumulated jacobian
231 0 : _femcontext.get_elem_jacobian() += old_jacobian;
232 0 : }
233 :
234 : // In DEBUG mode, we've set elem_jacobian == 0, and we
235 : // may have yet to add the old jacobian back
236 : #ifdef DEBUG
237 : else
238 : {
239 216536 : _femcontext.get_elem_jacobian() += old_jacobian;
240 : }
241 : #endif // ifdef DEBUG
242 2104838 : }
243 21265831 : }
244 :
245 21265831 : void add_element_system(FEMSystem & _sys,
246 : const bool _get_residual,
247 : const bool _get_jacobian,
248 : const bool _constrain_heterogeneously,
249 : const bool _no_constraints,
250 : FEMContext & _femcontext)
251 : {
252 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
253 21265831 : if (_get_residual && _sys.print_element_residuals)
254 : {
255 0 : std::streamsize old_precision = libMesh::out.precision();
256 0 : libMesh::out.precision(16);
257 0 : if (_femcontext.has_elem())
258 0 : libMesh::out << "Rraw_elem " << _femcontext.get_elem().id();
259 : else
260 0 : libMesh::out << "Rraw_scalar ";
261 0 : libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
262 0 : libMesh::out.precision(old_precision);
263 : }
264 :
265 21265831 : if (_get_jacobian && _sys.print_element_jacobians)
266 : {
267 0 : std::streamsize old_precision = libMesh::out.precision();
268 0 : libMesh::out.precision(16);
269 0 : if (_femcontext.has_elem())
270 0 : libMesh::out << "Jraw_elem " << _femcontext.get_elem().id();
271 : else
272 0 : libMesh::out << "Jraw_scalar ";
273 0 : libMesh::out << " = " << _femcontext.get_elem_jacobian() << std::endl;
274 0 : libMesh::out.precision(old_precision);
275 : }
276 :
277 : // We turn off the asymmetric constraint application iff we expect
278 : // enforce_constraints_exactly() to be called in the solver
279 21265831 : const bool constrain_in_solver = _sys.get_constrain_in_solver();
280 :
281 21265831 : if (_get_residual && _get_jacobian)
282 : {
283 4307493 : if (constrain_in_solver)
284 : {
285 3954375 : if (_constrain_heterogeneously)
286 0 : _sys.get_dof_map().heterogenously_constrain_element_matrix_and_vector
287 0 : (_femcontext.get_elem_jacobian(),
288 : _femcontext.get_elem_residual(),
289 : _femcontext.get_dof_indices(), false);
290 3954375 : else if (!_no_constraints)
291 351198 : _sys.get_dof_map().constrain_element_matrix_and_vector
292 3954375 : (_femcontext.get_elem_jacobian(),
293 : _femcontext.get_elem_residual(),
294 : _femcontext.get_dof_indices(), false);
295 : }
296 1760 : else if (!_no_constraints)
297 160 : _sys.get_dof_map().heterogeneously_constrain_element_jacobian_and_residual
298 1760 : (_femcontext.get_elem_jacobian(),
299 : _femcontext.get_elem_residual(),
300 : _femcontext.get_dof_indices(),
301 160 : *_sys.current_local_solution);
302 : // Do nothing if (_no_constraints)
303 : }
304 17309696 : else if (_get_residual)
305 : {
306 14788307 : if (constrain_in_solver)
307 : {
308 14778891 : if (_constrain_heterogeneously)
309 0 : _sys.get_dof_map().heterogenously_constrain_element_vector
310 0 : (_femcontext.get_elem_jacobian(),
311 : _femcontext.get_elem_residual(),
312 : _femcontext.get_dof_indices(), false);
313 14778891 : else if (!_no_constraints)
314 1305852 : _sys.get_dof_map().constrain_element_vector
315 14638995 : (_femcontext.get_elem_residual(),
316 : _femcontext.get_dof_indices(), false);
317 : }
318 9416 : else if (!_no_constraints)
319 160 : _sys.get_dof_map().heterogeneously_constrain_element_residual
320 1760 : (_femcontext.get_elem_residual(),
321 : _femcontext.get_dof_indices(),
322 160 : *_sys.current_local_solution);
323 : // Do nothing if (_no_constraints)
324 : }
325 2521389 : else if (_get_jacobian)
326 : {
327 : // Heterogeneous and homogeneous constraints are the same on the
328 : // matrix
329 : // Only get these contribs if we are applying some constraints
330 2521389 : if (!_no_constraints)
331 2749881 : _sys.get_dof_map().constrain_element_matrix (_femcontext.get_elem_jacobian(),
332 : _femcontext.get_dof_indices(),
333 2521389 : !constrain_in_solver);
334 : }
335 : #else
336 : libmesh_ignore(_constrain_heterogeneously, _no_constraints);
337 : #endif // #ifdef LIBMESH_ENABLE_CONSTRAINTS
338 :
339 21265831 : if (_get_residual && _sys.print_element_residuals)
340 : {
341 0 : std::streamsize old_precision = libMesh::out.precision();
342 0 : libMesh::out.precision(16);
343 0 : if (_femcontext.has_elem())
344 0 : libMesh::out << "R_elem " << _femcontext.get_elem().id();
345 : else
346 0 : libMesh::out << "R_scalar ";
347 0 : libMesh::out << " = " << _femcontext.get_elem_residual() << std::endl;
348 0 : libMesh::out.precision(old_precision);
349 : }
350 :
351 21265831 : if (_get_jacobian && _sys.print_element_jacobians)
352 : {
353 0 : std::streamsize old_precision = libMesh::out.precision();
354 0 : libMesh::out.precision(16);
355 0 : if (_femcontext.has_elem())
356 0 : libMesh::out << "J_elem " << _femcontext.get_elem().id();
357 : else
358 0 : libMesh::out << "J_scalar ";
359 0 : libMesh::out << " = " << _femcontext.get_elem_jacobian() << std::endl;
360 0 : libMesh::out.precision(old_precision);
361 : }
362 :
363 : { // A lock is necessary around access to the global system
364 3796780 : femsystem_mutex::scoped_lock lock(assembly_mutex);
365 :
366 21265831 : if (_get_jacobian)
367 7057374 : _sys.get_system_matrix().add_matrix (_femcontext.get_elem_jacobian(),
368 1159700 : _femcontext.get_dof_indices());
369 21265831 : if (_get_residual)
370 20414340 : _sys.rhs->add_vector (_femcontext.get_elem_residual(),
371 1669898 : _femcontext.get_dof_indices());
372 : } // Scope for assembly mutex
373 21265831 : }
374 :
375 :
376 :
377 : class AssemblyContributions
378 : {
379 : public:
380 : /**
381 : * constructor to set context
382 : */
383 6268 : AssemblyContributions(FEMSystem & sys,
384 : bool get_residual,
385 : bool get_jacobian,
386 : bool constrain_heterogeneously,
387 218854 : bool no_constraints) :
388 206318 : _sys(sys),
389 206318 : _get_residual(get_residual),
390 206318 : _get_jacobian(get_jacobian),
391 206318 : _constrain_heterogeneously(constrain_heterogeneously),
392 218854 : _no_constraints(no_constraints) {}
393 :
394 : /**
395 : * operator() for use with Threads::parallel_for().
396 : */
397 221050 : void operator()(const ConstElemRange & range) const
398 : {
399 228050 : std::unique_ptr<DiffContext> con = _sys.build_context();
400 7000 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
401 221050 : _sys.init_context(_femcontext);
402 :
403 21486881 : for (const auto & elem : range)
404 : {
405 21265831 : _femcontext.pre_fe_reinit(_sys, elem);
406 21265831 : _femcontext.elem_fe_reinit();
407 :
408 : assemble_unconstrained_element_system
409 21265831 : (_sys, _get_jacobian, _constrain_heterogeneously, _femcontext);
410 :
411 : add_element_system
412 25062611 : (_sys, _get_residual, _get_jacobian,
413 21265831 : _constrain_heterogeneously, _no_constraints, _femcontext);
414 : }
415 221050 : }
416 :
417 : private:
418 :
419 : FEMSystem & _sys;
420 :
421 : const bool _get_residual, _get_jacobian, _constrain_heterogeneously, _no_constraints;
422 : };
423 :
424 : class PostprocessContributions
425 : {
426 : public:
427 : /**
428 : * constructor to set context
429 : */
430 : explicit
431 1976 : PostprocessContributions(FEMSystem & sys) : _sys(sys) {}
432 :
433 : /**
434 : * operator() for use with Threads::parallel_for().
435 : */
436 1874 : void operator()(const ConstElemRange & range) const
437 : {
438 1976 : std::unique_ptr<DiffContext> con = _sys.build_context();
439 102 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
440 1874 : _sys.init_context(_femcontext);
441 :
442 64542 : for (const auto & elem : range)
443 : {
444 62668 : _femcontext.pre_fe_reinit(_sys, elem);
445 :
446 : // Optionally initialize all the interior FE objects on elem.
447 62668 : if (_sys.fe_reinit_during_postprocess)
448 62668 : _femcontext.elem_fe_reinit();
449 :
450 62668 : _sys.element_postprocess(_femcontext);
451 :
452 62668 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
453 310524 : for (_femcontext.side = 0; _femcontext.side != n_sides;
454 247856 : ++_femcontext.side)
455 : {
456 : // Don't compute on non-boundary sides unless requested
457 268760 : if (!_sys.postprocess_sides ||
458 316584 : (!_sys.get_physics()->compute_internal_sides &&
459 168744 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
460 232237 : continue;
461 :
462 : // Optionally initialize all the FE objects on this side.
463 15619 : if (_sys.fe_reinit_during_postprocess)
464 15619 : _femcontext.side_fe_reinit();
465 :
466 15619 : _sys.side_postprocess(_femcontext);
467 : }
468 : }
469 1874 : }
470 :
471 : private:
472 :
473 : FEMSystem & _sys;
474 : };
475 :
476 : class QoIContributions
477 : {
478 : public:
479 : /**
480 : * constructor to set context
481 : */
482 : explicit
483 39155 : QoIContributions(FEMSystem & sys,
484 : DifferentiableQoI & diff_qoi,
485 39155 : const QoISet & qoi_indices) :
486 40273 : qoi(sys.n_qois(), 0.), _sys(sys), _diff_qoi(diff_qoi),_qoi_indices(qoi_indices) {}
487 :
488 : /**
489 : * splitting constructor
490 : */
491 1152 : QoIContributions(const QoIContributions & other,
492 1152 : Threads::split) :
493 1152 : qoi(other._sys.n_qois(), 0.), _sys(other._sys), _diff_qoi(other._diff_qoi) {}
494 :
495 : /**
496 : * operator() for use with Threads::parallel_reduce().
497 : */
498 40883 : void operator()(const ConstElemRange & range)
499 : {
500 42577 : std::unique_ptr<DiffContext> con = _sys.build_context();
501 1694 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
502 40883 : _diff_qoi.init_context(_femcontext);
503 :
504 1694 : bool have_some_heterogenous_qoi_bc = false;
505 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
506 44271 : std::vector<bool> have_heterogenous_qoi_bc(_sys.n_qois(), false);
507 98566 : for (auto q : make_range(_sys.n_qois()))
508 59857 : if (_qoi_indices.has_index(q) &&
509 4348 : _sys.get_dof_map().has_heterogenous_adjoint_constraints(q))
510 : {
511 0 : have_heterogenous_qoi_bc[q] = true;
512 0 : have_some_heterogenous_qoi_bc = true;
513 : }
514 : #endif
515 :
516 40883 : if (have_some_heterogenous_qoi_bc)
517 0 : _sys.init_context(_femcontext);
518 :
519 13459487 : for (const auto & elem : range)
520 : {
521 13418604 : _femcontext.pre_fe_reinit(_sys, elem);
522 :
523 : // We might have some heterogenous dofs here; let's see for
524 : // certain
525 1216614 : bool elem_has_some_heterogenous_qoi_bc = false;
526 :
527 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
528 : const unsigned int n_dofs =
529 2433228 : cast_int<unsigned int>(_femcontext.get_dof_indices().size());
530 :
531 15851832 : std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.n_qois(), false);
532 13418604 : if (have_some_heterogenous_qoi_bc)
533 : {
534 0 : for (auto q : make_range(_sys.n_qois()))
535 : {
536 0 : if (have_heterogenous_qoi_bc[q])
537 : {
538 0 : for (auto d : make_range(n_dofs))
539 0 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
540 0 : (q, _femcontext.get_dof_indices()[d]) != Number(0))
541 : {
542 0 : elem_has_some_heterogenous_qoi_bc = true;
543 0 : elem_has_heterogenous_qoi_bc[q] = true;
544 0 : break;
545 : }
546 : }
547 : }
548 : }
549 : #endif
550 :
551 13418604 : if (_diff_qoi.assemble_qoi_elements ||
552 : elem_has_some_heterogenous_qoi_bc)
553 13418604 : _femcontext.elem_fe_reinit();
554 :
555 13418604 : if (_diff_qoi.assemble_qoi_elements)
556 13418604 : _diff_qoi.element_qoi(_femcontext, _qoi_indices);
557 :
558 : // If we have some heterogenous dofs here, those are
559 : // themselves part of a regularized flux QoI which the library
560 : // handles integrating
561 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
562 13418604 : if (elem_has_some_heterogenous_qoi_bc)
563 : {
564 0 : _sys.time_solver->element_residual(false, _femcontext);
565 :
566 0 : for (auto q : make_range(_sys.n_qois()))
567 : {
568 0 : if (elem_has_heterogenous_qoi_bc[q])
569 : {
570 0 : for (auto d : make_range(n_dofs))
571 0 : this->qoi[q] -= _femcontext.get_elem_residual()(d) *
572 0 : _sys.get_dof_map().has_heterogenous_adjoint_constraint(q, _femcontext.get_dof_indices()[d]);
573 :
574 : }
575 : }
576 : }
577 : #endif
578 :
579 13418604 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
580 67093020 : for (_femcontext.side = 0; _femcontext.side != n_sides;
581 53674416 : ++_femcontext.side)
582 : {
583 : // Don't compute on non-boundary sides unless requested
584 53675400 : if (!_diff_qoi.assemble_qoi_sides ||
585 11184 : (!_diff_qoi.assemble_qoi_internal_sides &&
586 11184 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
587 48805234 : continue;
588 :
589 2999 : _femcontext.side_fe_reinit();
590 :
591 2999 : _diff_qoi.side_qoi(_femcontext, _qoi_indices);
592 : }
593 : }
594 :
595 42577 : this->_diff_qoi.thread_join( this->qoi, _femcontext.get_qois(), _qoi_indices );
596 40883 : }
597 :
598 1152 : void join (const QoIContributions & other)
599 : {
600 576 : libmesh_assert_equal_to (this->qoi.size(), other.qoi.size());
601 1728 : this->_diff_qoi.thread_join( this->qoi, other.qoi, _qoi_indices );
602 1152 : }
603 :
604 : std::vector<Number> qoi;
605 :
606 : private:
607 :
608 : FEMSystem & _sys;
609 : DifferentiableQoI & _diff_qoi;
610 :
611 : const QoISet _qoi_indices;
612 : };
613 :
614 : class QoIDerivativeContributions
615 : {
616 : public:
617 : /**
618 : * constructor to set context
619 : */
620 464 : QoIDerivativeContributions(FEMSystem & sys,
621 : const QoISet & qoi_indices,
622 : DifferentiableQoI & qoi,
623 : bool include_liftfunc,
624 14654 : bool apply_constraints) :
625 13726 : _sys(sys),
626 13726 : _qoi_indices(qoi_indices),
627 13726 : _qoi(qoi),
628 13726 : _include_liftfunc(include_liftfunc),
629 14654 : _apply_constraints(apply_constraints) {}
630 :
631 : /**
632 : * operator() for use with Threads::parallel_for().
633 : */
634 14948 : void operator()(const ConstElemRange & range) const
635 : {
636 15510 : std::unique_ptr<DiffContext> con = _sys.build_context();
637 562 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
638 14948 : _qoi.init_context(_femcontext);
639 :
640 562 : bool have_some_heterogenous_qoi_bc = false;
641 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
642 16072 : std::vector<bool> have_heterogenous_qoi_bc(_sys.n_qois(), false);
643 14948 : if (_include_liftfunc || _apply_constraints)
644 40059 : for (auto q : make_range(_sys.n_qois()))
645 26009 : if (_qoi_indices.has_index(q) &&
646 1796 : _sys.get_dof_map().has_heterogenous_adjoint_constraints(q))
647 : {
648 28 : have_heterogenous_qoi_bc[q] = true;
649 14 : have_some_heterogenous_qoi_bc = true;
650 : }
651 : #endif
652 :
653 14948 : if (have_some_heterogenous_qoi_bc)
654 182 : _sys.init_context(_femcontext);
655 :
656 2430969 : for (const auto & elem : range)
657 : {
658 2416021 : _femcontext.pre_fe_reinit(_sys, elem);
659 :
660 : // We might have some heterogenous dofs here; let's see for
661 : // certain
662 219044 : bool elem_has_some_heterogenous_qoi_bc = false;
663 :
664 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
665 : const unsigned int n_dofs =
666 438088 : cast_int<unsigned int>(_femcontext.get_dof_indices().size());
667 :
668 2854109 : std::vector<bool> elem_has_heterogenous_qoi_bc(_sys.n_qois(), false);
669 2416021 : if (have_some_heterogenous_qoi_bc)
670 : {
671 85476 : for (auto q : make_range(_sys.n_qois()))
672 : {
673 48450 : if (have_heterogenous_qoi_bc[q])
674 : {
675 412332 : for (auto d : make_range(n_dofs))
676 320980 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
677 421568 : (q, _femcontext.get_dof_indices()[d]) != Number(0))
678 : {
679 140 : elem_has_some_heterogenous_qoi_bc = true;
680 280 : elem_has_heterogenous_qoi_bc[q] = true;
681 1680 : break;
682 : }
683 : }
684 : }
685 : }
686 : #endif
687 :
688 : // If we're going to call a user integral, then we need FE
689 : // information to call element_qoi.
690 : // If we're going to evaluate lift-function-based components
691 : // of a QoI, then we need FE information to assemble the
692 : // element residual.
693 2416021 : if (_qoi.assemble_qoi_elements ||
694 0 : ((_include_liftfunc || _apply_constraints) &&
695 : elem_has_some_heterogenous_qoi_bc))
696 2416021 : _femcontext.elem_fe_reinit();
697 :
698 2416021 : if (_qoi.assemble_qoi_elements)
699 2416021 : _qoi.element_qoi_derivative(_femcontext, _qoi_indices);
700 :
701 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
702 : // If we need to use heterogenous dofs here, we need the
703 : // Jacobian either for the regularized flux QoI integration
704 : // and/or for constraint application.
705 2416021 : if ((_include_liftfunc || _apply_constraints) &&
706 : elem_has_some_heterogenous_qoi_bc)
707 : {
708 1820 : bool jacobian_computed = _sys.time_solver->element_residual(true, _femcontext);
709 :
710 : // If we're using numerical jacobians, above wont compute them
711 1680 : if (!jacobian_computed)
712 : {
713 : // Make sure we didn't compute a jacobian and lie about it
714 140 : libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
715 : // Logging of numerical jacobians is done separately
716 1680 : _sys.numerical_elem_jacobian(_femcontext);
717 : }
718 : }
719 :
720 : // If we have some heterogenous dofs here, those are
721 : // themselves part of a regularized flux QoI which the library
722 : // may handle integrating
723 2416021 : if (_include_liftfunc && elem_has_some_heterogenous_qoi_bc)
724 : {
725 0 : for (auto q : make_range(_sys.n_qois()))
726 : {
727 0 : if (elem_has_heterogenous_qoi_bc[q])
728 : {
729 0 : for (auto i : make_range(n_dofs))
730 : {
731 : Number liftfunc_val =
732 0 : _sys.get_dof_map().has_heterogenous_adjoint_constraint(q, _femcontext.get_dof_indices()[i]);
733 :
734 0 : if (liftfunc_val != Number(0))
735 : {
736 0 : for (auto j : make_range(n_dofs))
737 0 : _femcontext.get_qoi_derivatives()[q](j) -=
738 0 : _femcontext.get_elem_jacobian()(i,j) *
739 : liftfunc_val;
740 : }
741 : }
742 : }
743 : }
744 : }
745 : #endif
746 :
747 :
748 2416021 : const unsigned char n_sides = _femcontext.get_elem().n_sides();
749 12080105 : for (_femcontext.side = 0; _femcontext.side != n_sides;
750 9664084 : ++_femcontext.side)
751 : {
752 : // Don't compute on non-boundary sides unless requested
753 9723300 : if (!_qoi.assemble_qoi_sides ||
754 648916 : (!_qoi.assemble_qoi_internal_sides &&
755 648916 : _femcontext.get_elem().neighbor_ptr(_femcontext.side) != nullptr))
756 8760300 : continue;
757 :
758 30968 : _femcontext.side_fe_reinit();
759 :
760 30968 : _qoi.side_qoi_derivative(_femcontext, _qoi_indices);
761 : }
762 :
763 : // We need some unmodified indices to use for constraining
764 : // multiple vector
765 : // FIXME - there should be a DofMap::constrain_element_vectors
766 : // to do this more efficiently
767 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
768 2635065 : std::vector<dof_id_type> original_dofs = _femcontext.get_dof_indices();
769 : #endif
770 :
771 : { // A lock is necessary around access to the global system
772 438088 : femsystem_mutex::scoped_lock lock(assembly_mutex);
773 :
774 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
775 : // We'll need to see if any heterogenous constraints apply
776 : // to the QoI dofs on this element *or* to any of the dofs
777 : // they depend on, so let's get those dependencies
778 2416021 : if (_apply_constraints)
779 2595129 : _sys.get_dof_map().constrain_nothing(_femcontext.get_dof_indices());
780 : #endif
781 :
782 4929086 : for (auto i : make_range(_sys.n_qois()))
783 2513065 : if (_qoi_indices.has_index(i))
784 : {
785 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
786 2513065 : if (_apply_constraints)
787 : {
788 : #ifndef NDEBUG
789 225374 : bool has_heterogenous_constraint = false;
790 1268386 : for (auto d : make_range(n_dofs))
791 1043152 : if (_sys.get_dof_map().has_heterogenous_adjoint_constraint
792 1043152 : (i, _femcontext.get_dof_indices()[d]) != Number(0))
793 : {
794 140 : has_heterogenous_constraint = true;
795 140 : libmesh_assert(elem_has_heterogenous_qoi_bc[i]);
796 140 : libmesh_assert(elem_has_some_heterogenous_qoi_bc);
797 140 : break;
798 : }
799 : #else
800 : bool has_heterogenous_constraint =
801 225374 : elem_has_heterogenous_qoi_bc[i];
802 : #endif
803 :
804 2476201 : _femcontext.get_dof_indices() = original_dofs;
805 :
806 2476201 : if (has_heterogenous_constraint)
807 : {
808 : // Q_u gets used for *adjoint* solves, so we
809 : // need K^T here.
810 1960 : DenseMatrix<Number> elem_jacobian_transpose;
811 140 : _femcontext.get_elem_jacobian().get_transpose
812 1680 : (elem_jacobian_transpose);
813 :
814 1680 : _sys.get_dof_map().heterogenously_constrain_element_vector
815 1820 : (elem_jacobian_transpose,
816 280 : _femcontext.get_qoi_derivatives()[i],
817 : _femcontext.get_dof_indices(), false, i);
818 1400 : }
819 : else
820 : {
821 2474521 : _sys.get_dof_map().constrain_element_vector
822 2699755 : (_femcontext.get_qoi_derivatives()[i],
823 : _femcontext.get_dof_indices(), false);
824 : }
825 : }
826 : #endif
827 :
828 2513065 : _sys.get_adjoint_rhs(i).add_vector
829 2513065 : (_femcontext.get_qoi_derivatives()[i], _femcontext.get_dof_indices());
830 : }
831 : }
832 : }
833 14948 : }
834 :
835 : private:
836 :
837 : FEMSystem & _sys;
838 : const QoISet & _qoi_indices;
839 : DifferentiableQoI & _qoi;
840 : bool _include_liftfunc, _apply_constraints;
841 : };
842 :
843 :
844 : }
845 :
846 :
847 : namespace libMesh
848 : {
849 :
850 :
851 :
852 :
853 :
854 6647 : FEMSystem::FEMSystem (EquationSystems & es,
855 : const std::string & name_in,
856 0 : const unsigned int number_in)
857 : : Parent(es, name_in, number_in),
858 6279 : fe_reinit_during_postprocess(true),
859 6279 : numerical_jacobian_h(TOLERANCE),
860 6647 : verify_analytic_jacobians(0.0)
861 : {
862 6647 : }
863 :
864 :
865 6279 : FEMSystem::~FEMSystem () = default;
866 :
867 :
868 :
869 6647 : void FEMSystem::init_data ()
870 : {
871 : // First initialize LinearImplicitSystem data
872 6647 : Parent::init_data();
873 6647 : }
874 :
875 :
876 218854 : void FEMSystem::assembly (bool get_residual, bool get_jacobian,
877 : bool apply_heterogeneous_constraints,
878 : bool apply_no_constraints)
879 : {
880 6268 : libmesh_assert(get_residual || get_jacobian);
881 :
882 : // Log residual and jacobian and combined performance separately
883 : #ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
884 : const char * log_name;
885 12536 : if (get_residual && get_jacobian)
886 1374 : log_name = "assembly()";
887 9788 : else if (get_residual)
888 4242 : log_name = "assembly(get_residual)";
889 : else
890 652 : log_name = "assembly(get_jacobian)";
891 :
892 12536 : LOG_SCOPE(log_name, "FEMSystem");
893 : #endif
894 :
895 12536 : const MeshBase & mesh = this->get_mesh();
896 :
897 6268 : libmesh_assert(mesh.is_prepared());
898 : #if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED)
899 : MeshTools::libmesh_assert_valid_is_prepared(mesh);
900 : #endif
901 :
902 218854 : if (print_solution_norms)
903 : {
904 0 : this->solution->close();
905 :
906 0 : std::streamsize old_precision = libMesh::out.precision();
907 0 : libMesh::out.precision(16);
908 0 : libMesh::out << "|U| = "
909 0 : << this->solution->l1_norm()
910 0 : << std::endl;
911 0 : libMesh::out.precision(old_precision);
912 : }
913 218854 : if (print_solutions)
914 : {
915 0 : std::streamsize old_precision = libMesh::out.precision();
916 0 : libMesh::out.precision(16);
917 0 : libMesh::out << "U = [" << *(this->solution)
918 0 : << "];" << std::endl;
919 0 : libMesh::out.precision(old_precision);
920 : }
921 :
922 : // Is this definitely necessary? [RHS]
923 : // Yes. [RHS 2012]
924 218854 : if (get_jacobian)
925 69277 : matrix->zero();
926 218854 : if (get_residual)
927 197530 : rhs->zero();
928 :
929 : // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
930 218854 : if (verify_analytic_jacobians > 0.5)
931 : {
932 0 : libMesh::err << "WARNING! verify_analytic_jacobians was set "
933 0 : << "to absurdly large value of "
934 0 : << verify_analytic_jacobians << std::endl;
935 0 : libMesh::err << "Resetting to 1e-6!" << std::endl;
936 0 : verify_analytic_jacobians = 1e-6;
937 : }
938 :
939 : // In time-dependent problems, the nonlinear function we're trying
940 : // to solve at each timestep may depend on the particular solver
941 : // we're using
942 6268 : libmesh_assert(time_solver.get());
943 :
944 : // Build the residual and jacobian contributions on every active
945 : // mesh element on this processor
946 : Threads::parallel_for
947 225122 : (mesh.active_local_element_stored_range(),
948 218854 : 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 6268 : bool have_scalar = false;
954 481238 : for (auto i : make_range(this->n_variable_groups()))
955 : {
956 262384 : 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 225122 : 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 218854 : if (get_residual && (print_residual_norms || print_residuals))
1040 0 : this->rhs->close();
1041 198834 : 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 198834 : 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 218854 : if (get_jacobian && (print_jacobian_norms || print_jacobians))
1057 0 : this->matrix->close();
1058 77761 : 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 77761 : 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 218854 : }
1073 :
1074 :
1075 :
1076 26782 : void FEMSystem::solve()
1077 : {
1078 : // We are solving the primal problem
1079 26782 : Parent::solve();
1080 :
1081 : // On a moving mesh we want the mesh to reflect the new solution
1082 26782 : this->mesh_position_set();
1083 26782 : }
1084 :
1085 :
1086 :
1087 29057 : void FEMSystem::mesh_position_set()
1088 : {
1089 : // If we don't need to move the mesh, we're done
1090 29057 : if (_mesh_sys != this)
1091 26457 : 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 2078 : Threads::parallel_for (mesh.active_local_element_stored_range(),
1141 2078 : PostprocessContributions(*this));
1142 1976 : }
1143 :
1144 :
1145 :
1146 39155 : void FEMSystem::assemble_qoi (const QoISet & qoi_indices)
1147 : {
1148 2236 : LOG_SCOPE("assemble_qoi()", "FEMSystem");
1149 :
1150 2236 : const MeshBase & mesh = this->get_mesh();
1151 :
1152 1118 : libmesh_assert(mesh.is_prepared());
1153 : #if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED)
1154 : MeshTools::libmesh_assert_valid_is_prepared(mesh);
1155 : #endif
1156 :
1157 39155 : this->update();
1158 :
1159 1118 : const unsigned int Nq = this->n_qois();
1160 :
1161 : // the quantity of interest is assumed to be a sum of element and
1162 : // side terms
1163 95110 : for (unsigned int i=0; i != Nq; ++i)
1164 55955 : if (qoi_indices.has_index(i))
1165 55955 : this->set_qoi(i, 0);
1166 :
1167 : // Create a non-temporary qoi_contributions object, so we can query
1168 : // its results after the reduction
1169 79428 : QoIContributions qoi_contributions(*this, *(this->get_qoi()), qoi_indices);
1170 :
1171 : // Loop over every active mesh element on this processor
1172 39155 : Threads::parallel_reduce(mesh.active_local_element_stored_range(),
1173 : qoi_contributions);
1174 :
1175 39155 : std::vector<Number> global_qoi = this->get_qoi_values();
1176 39155 : this->get_qoi()->parallel_op( this->comm(), global_qoi, qoi_contributions.qoi, qoi_indices );
1177 77192 : this->set_qoi(std::move(global_qoi));
1178 39155 : }
1179 :
1180 :
1181 :
1182 14654 : void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices,
1183 : bool include_liftfunc,
1184 : bool apply_constraints)
1185 : {
1186 928 : LOG_SCOPE("assemble_qoi_derivative()", "FEMSystem");
1187 :
1188 928 : const MeshBase & mesh = this->get_mesh();
1189 :
1190 464 : libmesh_assert(mesh.is_prepared());
1191 : #if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED)
1192 : MeshTools::libmesh_assert_valid_is_prepared(mesh);
1193 : #endif
1194 :
1195 14654 : this->update();
1196 :
1197 : // The quantity of interest derivative assembly accumulates on
1198 : // initially zero vectors
1199 39471 : for (auto i : make_range(this->n_qois()))
1200 24817 : if (qoi_indices.has_index(i))
1201 24817 : this->add_adjoint_rhs(i).zero();
1202 :
1203 : // Loop over every active mesh element on this processor
1204 15118 : Threads::parallel_for (mesh.active_local_element_stored_range(),
1205 14190 : QoIDerivativeContributions(*this, qoi_indices,
1206 14654 : *(this->get_qoi()),
1207 : include_liftfunc,
1208 : apply_constraints));
1209 :
1210 39471 : for (auto i : make_range(this->n_qois()))
1211 24817 : if (qoi_indices.has_index(i))
1212 24817 : this->get_qoi()->finalize_derivative(this->get_adjoint_rhs(i),i);
1213 14654 : }
1214 :
1215 :
1216 :
1217 277926 : void FEMSystem::numerical_jacobian (TimeSolverResPtr res,
1218 : FEMContext & context) const
1219 : {
1220 : // Logging is done by numerical_elem_jacobian
1221 : // or numerical_side_jacobian
1222 :
1223 54472 : DenseVector<Number> original_residual(context.get_elem_residual());
1224 54472 : DenseVector<Number> backwards_residual(context.get_elem_residual());
1225 332398 : DenseMatrix<Number> numeric_jacobian(context.get_elem_jacobian());
1226 : #ifdef DEBUG
1227 54472 : DenseMatrix<Number> old_jacobian(context.get_elem_jacobian());
1228 : #endif
1229 :
1230 27236 : Real numerical_point_h = 0.;
1231 277926 : if (_mesh_sys == this)
1232 0 : numerical_point_h = numerical_jacobian_h * context.get_elem().hmin();
1233 :
1234 : const unsigned int n_dofs =
1235 54472 : cast_int<unsigned int>(context.get_dof_indices().size());
1236 :
1237 608124 : for (auto v : make_range(context.n_vars()))
1238 : {
1239 298606 : const Real my_h = this->numerical_jacobian_h_for_var(v);
1240 :
1241 31592 : unsigned int j_offset = libMesh::invalid_uint;
1242 :
1243 330198 : if (!context.get_dof_indices(v).empty())
1244 : {
1245 4318140 : for (auto i : make_range(n_dofs))
1246 4725958 : if (context.get_dof_indices()[i] ==
1247 369008 : context.get_dof_indices(v)[0])
1248 31592 : j_offset = i;
1249 :
1250 31592 : libmesh_assert_not_equal_to(j_offset, libMesh::invalid_uint);
1251 : }
1252 :
1253 3172188 : for (auto j : make_range(context.get_dof_indices(v).size()))
1254 : {
1255 2841990 : const unsigned int total_j = j + j_offset;
1256 :
1257 : // Take the "minus" side of a central differenced first derivative
1258 2841990 : Number original_solution = context.get_elem_solution(v)(j);
1259 2595928 : context.get_elem_solution(v)(j) -= my_h;
1260 :
1261 : // Make sure to catch any moving mesh terms
1262 273512 : Real * coord = nullptr;
1263 2841990 : if (_mesh_sys == this)
1264 : {
1265 0 : if (_mesh_x_var == v)
1266 0 : coord = &(context.get_elem().point(j)(0));
1267 0 : else if (_mesh_y_var == v)
1268 0 : coord = &(context.get_elem().point(j)(1));
1269 0 : else if (_mesh_z_var == v)
1270 0 : coord = &(context.get_elem().point(j)(2));
1271 : }
1272 273512 : if (coord)
1273 : {
1274 : // We have enough information to scale the perturbations
1275 : // here appropriately
1276 0 : context.get_elem_solution(v)(j) = original_solution - numerical_point_h;
1277 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1278 : }
1279 :
1280 273512 : context.get_elem_residual().zero();
1281 2841990 : ((*time_solver).*(res))(false, context);
1282 : #ifdef DEBUG
1283 273512 : libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1284 : #endif
1285 273512 : backwards_residual = context.get_elem_residual();
1286 :
1287 : // Take the "plus" side of a central differenced first derivative
1288 2841990 : context.get_elem_solution(v)(j) = original_solution + my_h;
1289 2841990 : if (coord)
1290 : {
1291 0 : context.get_elem_solution()(j) = original_solution + numerical_point_h;
1292 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1293 : }
1294 273512 : context.get_elem_residual().zero();
1295 2841990 : ((*time_solver).*(res))(false, context);
1296 : #ifdef DEBUG
1297 273512 : libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
1298 : #endif
1299 :
1300 2841990 : context.get_elem_solution(v)(j) = original_solution;
1301 2841990 : if (coord)
1302 : {
1303 0 : *coord = libmesh_real(context.get_elem_solution(v)(j));
1304 0 : for (auto i : make_range(n_dofs))
1305 : {
1306 0 : numeric_jacobian(i,total_j) =
1307 0 : (context.get_elem_residual()(i) - backwards_residual(i)) /
1308 0 : 2. / numerical_point_h;
1309 : }
1310 : }
1311 : else
1312 : {
1313 37004700 : for (auto i : make_range(n_dofs))
1314 : {
1315 34162710 : numeric_jacobian(i,total_j) =
1316 31232752 : (context.get_elem_residual()(i) - backwards_residual(i)) /
1317 31232752 : 2. / my_h;
1318 : }
1319 : }
1320 : }
1321 : }
1322 :
1323 27236 : context.get_elem_residual() = original_residual;
1324 277926 : context.get_elem_jacobian() = numeric_jacobian;
1325 501380 : }
1326 :
1327 :
1328 :
1329 228744 : void FEMSystem::numerical_elem_jacobian (FEMContext & context) const
1330 : {
1331 43224 : LOG_SCOPE("numerical_elem_jacobian()", "FEMSystem");
1332 228744 : this->numerical_jacobian(&TimeSolver::element_residual, context);
1333 228744 : }
1334 :
1335 :
1336 :
1337 49182 : void FEMSystem::numerical_side_jacobian (FEMContext & context) const
1338 : {
1339 11248 : LOG_SCOPE("numerical_side_jacobian()", "FEMSystem");
1340 49182 : this->numerical_jacobian(&TimeSolver::side_residual, context);
1341 49182 : }
1342 :
1343 :
1344 :
1345 0 : void FEMSystem::numerical_nonlocal_jacobian (FEMContext & context) const
1346 : {
1347 0 : LOG_SCOPE("numerical_nonlocal_jacobian()", "FEMSystem");
1348 0 : this->numerical_jacobian(&TimeSolver::nonlocal_residual, context);
1349 0 : }
1350 :
1351 :
1352 :
1353 317376 : std::unique_ptr<DiffContext> FEMSystem::build_context ()
1354 : {
1355 327738 : auto fc = std::make_unique<FEMContext>(*this);
1356 :
1357 317376 : DifferentiablePhysics * phys = this->get_physics();
1358 :
1359 10362 : libmesh_assert (phys);
1360 :
1361 : // If we are solving a moving mesh problem, tell that to the Context
1362 317376 : fc->set_mesh_system(phys->get_mesh_system());
1363 20724 : fc->set_mesh_x_var(phys->get_mesh_x_var());
1364 20724 : fc->set_mesh_y_var(phys->get_mesh_y_var());
1365 20724 : fc->set_mesh_z_var(phys->get_mesh_z_var());
1366 :
1367 317376 : fc->set_deltat_pointer( &deltat );
1368 :
1369 : // If we are solving the adjoint problem, tell that to the Context
1370 327738 : fc->is_adjoint() = this->get_time_solver().is_adjoint();
1371 :
1372 327738 : return fc;
1373 296652 : }
1374 :
1375 :
1376 :
1377 246841 : void FEMSystem::init_context(DiffContext & c)
1378 : {
1379 : // Parent::init_context(c); // may be a good idea in derived classes
1380 :
1381 : // Although we do this in DiffSystem::build_context() and
1382 : // FEMSystem::build_context() as well, we do it here just to be
1383 : // extra sure that the deltat pointer gets set. Since the
1384 : // intended behavior is for classes derived from FEMSystem to
1385 : // call Parent::init_context() in their own init_context()
1386 : // overloads, we can ensure that those classes get the correct
1387 : // deltat pointers even if they have different build_context()
1388 : // overloads.
1389 246841 : c.set_deltat_pointer ( &deltat );
1390 :
1391 8286 : FEMContext & context = cast_ref<FEMContext &>(c);
1392 :
1393 : // Make sure we're prepared to do mass integration
1394 594642 : for (auto var : make_range(this->n_vars()))
1395 684484 : if (this->get_physics()->is_time_evolving(var))
1396 : {
1397 : // Request shape functions based on FEType
1398 182801 : switch( FEInterface::field_type( this->variable_type(var) ) )
1399 : {
1400 182801 : case( TYPE_SCALAR ):
1401 : {
1402 6494 : FEBase * elem_fe = nullptr;
1403 366667 : for (auto dim : context.elem_dimensions())
1404 : {
1405 6524 : context.get_element_fe(var, elem_fe, dim);
1406 13374 : elem_fe->get_JxW();
1407 6524 : elem_fe->get_phi();
1408 0 : }
1409 : }
1410 6494 : break;
1411 0 : case( TYPE_VECTOR ):
1412 : {
1413 0 : FEGenericBase<RealGradient> * elem_fe = nullptr;
1414 0 : for (auto dim : context.elem_dimensions())
1415 : {
1416 0 : context.get_element_fe(var, elem_fe, dim);
1417 0 : elem_fe->get_JxW();
1418 0 : elem_fe->get_phi();
1419 0 : }
1420 : }
1421 0 : break;
1422 0 : default:
1423 0 : libmesh_error_msg("Unrecognized field type!");
1424 : }
1425 : }
1426 246841 : }
1427 :
1428 :
1429 :
1430 65 : void FEMSystem::mesh_position_get()
1431 : {
1432 : // This function makes no sense unless we've already picked out some
1433 : // variable(s) to reflect mesh position coordinates
1434 65 : libmesh_error_msg_if(!_mesh_sys, "_mesh_sys was nullptr!");
1435 :
1436 : // We currently assume mesh variables are in our own system
1437 65 : if (_mesh_sys != this)
1438 0 : libmesh_not_implemented();
1439 :
1440 : // Loop over every active mesh element on this processor
1441 0 : const MeshBase & mesh = this->get_mesh();
1442 :
1443 65 : std::unique_ptr<DiffContext> con = this->build_context();
1444 0 : FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
1445 65 : this->init_context(_femcontext);
1446 :
1447 : // Get the solution's mesh variables from every element
1448 1282 : for (const auto & elem : mesh.active_local_element_ptr_range())
1449 : {
1450 576 : _femcontext.pre_fe_reinit(*this, elem);
1451 :
1452 576 : _femcontext.elem_position_get();
1453 :
1454 576 : if (_mesh_x_var != libMesh::invalid_uint)
1455 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_x_var),
1456 0 : _femcontext.get_dof_indices(_mesh_x_var) );
1457 576 : if (_mesh_y_var != libMesh::invalid_uint)
1458 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_y_var),
1459 0 : _femcontext.get_dof_indices(_mesh_y_var));
1460 576 : if (_mesh_z_var != libMesh::invalid_uint)
1461 576 : this->solution->insert(_femcontext.get_elem_solution(_mesh_z_var),
1462 0 : _femcontext.get_dof_indices(_mesh_z_var));
1463 65 : }
1464 :
1465 65 : this->solution->close();
1466 :
1467 : // And make sure the current_local_solution is up to date too
1468 65 : this->System::update();
1469 65 : }
1470 :
1471 : } // namespace libMesh
|
