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