Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // C++ includes
19 : #include <algorithm> // for std::fill
20 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
21 : #include <cmath> // for sqrt
22 : #include <set>
23 : #include <sstream> // for ostringstream
24 : #include <unordered_map>
25 : #include <unordered_set>
26 :
27 : // Local Includes
28 : #include "libmesh/dof_map.h"
29 : #include "libmesh/elem.h"
30 : #include "libmesh/equation_systems.h"
31 : #include "libmesh/error_vector.h"
32 : #include "libmesh/fe.h"
33 : #include "libmesh/fe_interface.h"
34 : #include "libmesh/libmesh_common.h"
35 : #include "libmesh/libmesh_logging.h"
36 : #include "libmesh/mesh_base.h"
37 : #include "libmesh/mesh_refinement.h"
38 : #include "libmesh/numeric_vector.h"
39 : #include "libmesh/quadrature.h"
40 : #include "libmesh/system.h"
41 : #include "libmesh/diff_physics.h"
42 : #include "libmesh/fem_system.h"
43 : #include "libmesh/implicit_system.h"
44 : #include "libmesh/partitioner.h"
45 : #include "libmesh/adjoint_refinement_estimator.h"
46 : #include "libmesh/enum_error_estimator_type.h"
47 : #include "libmesh/enum_norm_type.h"
48 : #include "libmesh/utility.h"
49 :
50 : #ifdef LIBMESH_ENABLE_AMR
51 :
52 : namespace libMesh
53 : {
54 :
55 : //-----------------------------------------------------------------
56 : // AdjointRefinementErrorEstimator
57 :
58 : // As of 10/2/2012, this function implements a 'brute-force' adjoint based QoI
59 : // error estimator, using the following relationship:
60 : // Q(u) - Q(u_h) \approx - R( (u_h)_(h/2), z_(h/2) ) .
61 : // where: Q(u) is the true QoI
62 : // u_h is the approximate primal solution on the current FE space
63 : // Q(u_h) is the approximate QoI
64 : // z_(h/2) is the adjoint corresponding to Q, on a richer FE space
65 : // (u_h)_(h/2) is a projection of the primal solution on the richer FE space
66 : // By richer FE space, we mean a grid that has been refined once and a polynomial order
67 : // that has been increased once, i.e. one h and one p refinement
68 :
69 : // Both a global QoI error estimate and element wise error indicators are included
70 : // Note that the element wise error indicators slightly over estimate the error in
71 : // each element
72 :
73 986 : AdjointRefinementEstimator::AdjointRefinementEstimator() :
74 : ErrorEstimator(),
75 930 : number_h_refinements(1),
76 930 : number_p_refinements(0),
77 930 : _residual_evaluation_physics(nullptr),
78 930 : _adjoint_evaluation_physics(nullptr),
79 986 : _qoi_set(QoISet())
80 : {
81 : // We're not actually going to use error_norm; our norms are
82 : // absolute values of QoI error.
83 986 : error_norm = INVALID_NORM;
84 986 : }
85 :
86 0 : ErrorEstimatorType AdjointRefinementEstimator::type() const
87 : {
88 0 : return ADJOINT_REFINEMENT;
89 : }
90 :
91 16666 : void AdjointRefinementEstimator::estimate_error (const System & _system,
92 : ErrorVector & error_per_cell,
93 : const NumericVector<Number> * solution_vector,
94 : bool /*estimate_parent_error*/)
95 : {
96 : // We have to break the rules here, because we can't refine a const System
97 476 : System & mutable_system = const_cast<System &>(_system);
98 :
99 : // We really have to break the rules, because we can't do an
100 : // adjoint_solve without a matrix.
101 16666 : ImplicitSystem & system = dynamic_cast<ImplicitSystem &>(mutable_system);
102 :
103 : // An EquationSystems reference will be convenient.
104 952 : EquationSystems & es = system.get_equation_systems();
105 :
106 : // The current mesh
107 952 : MeshBase & mesh = es.get_mesh();
108 :
109 : // Get coarse grid adjoint solutions. This should be a relatively
110 : // quick (especially with preconditioner reuse) way to get a good
111 : // initial guess for the fine grid adjoint solutions. More
112 : // importantly, subtracting off a coarse adjoint approximation gives
113 : // us better local error indication, and subtracting off *some* lift
114 : // function is necessary for correctness if we have heterogeneous
115 : // adjoint Dirichlet conditions.
116 :
117 : // Solve the adjoint problem(s) on the coarse FE space
118 : // Only if the user didn't already solve it for us
119 : // If _adjoint_evaluation_physics pointer is not null, swap
120 : // the current physics with the one held by _adjoint_evaluation_physics
121 : // before assembling the adjoint problem
122 16666 : if (!system.is_adjoint_already_solved())
123 : {
124 : // Swap in different physics if needed
125 0 : if (_adjoint_evaluation_physics)
126 0 : dynamic_cast<DifferentiableSystem &>(system).push_physics(*_adjoint_evaluation_physics);
127 :
128 : // Solve the adjoint problem, remember physics swap also resets the cache, so
129 : // we will assemble again, otherwise we just take the transpose
130 0 : system.adjoint_solve(_qoi_set);
131 :
132 0 : if (_adjoint_evaluation_physics)
133 0 : dynamic_cast<DifferentiableSystem &>(system).pop_physics();
134 : }
135 :
136 : // Loop over all the adjoint problems and, if any have heterogenous
137 : // Dirichlet conditions, get the corresponding coarse lift
138 : // function(s)
139 48833 : for (unsigned int j=0,
140 49785 : n_qois = system.n_qois(); j != n_qois; j++)
141 : {
142 : // Skip this QoI if it is not in the QoI Set or if there are no
143 : // heterogeneous Dirichlet boundaries for it
144 66238 : if (_qoi_set.has_index(j) &&
145 33119 : system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
146 : {
147 : // Next, we are going to build up the residual for evaluating the flux QoI
148 934 : NumericVector<Number> * coarse_residual = nullptr;
149 :
150 : // The definition of a flux QoI is R(u^h, L) where R is the residual as defined
151 : // by a conservation law. Therefore, if we are using stabilization, the
152 : // R should be specified by the user via the residual_evaluation_physics
153 :
154 : // If the residual physics pointer is not null, use it when
155 : // evaluating here.
156 : {
157 32693 : if (_residual_evaluation_physics)
158 32480 : dynamic_cast<DifferentiableSystem &>(system).push_physics(*_residual_evaluation_physics);
159 :
160 : // Assemble without applying constraints, to capture the solution values on the boundary
161 32693 : system.assembly(true, false, false, true);
162 :
163 : // Get the residual vector (no constraints applied on boundary, so we wont blow away the lift)
164 32693 : coarse_residual = &system.get_vector("RHS Vector");
165 32693 : coarse_residual->close();
166 :
167 : // Now build the lift function and add it to the system vectors
168 34561 : std::ostringstream liftfunc_name;
169 31759 : liftfunc_name << "adjoint_lift_function" << j;
170 :
171 64452 : system.add_vector(liftfunc_name.str());
172 :
173 : // Initialize lift with coarse adjoint solve associate with this flux QoI to begin with
174 64452 : system.get_vector(liftfunc_name.str()).init(system.get_adjoint_solution(j), false);
175 :
176 : // Build the actual lift using adjoint dirichlet conditions
177 934 : system.get_dof_map().enforce_adjoint_constraints_exactly
178 33627 : (system.get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
179 :
180 : // Compute the flux R(u^h, L)
181 63518 : std::cout<<"The flux QoI "<<static_cast<unsigned int>(j)<<" is: "<<coarse_residual->dot(system.get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
182 :
183 : // Restore the original physics if needed
184 32693 : if (_residual_evaluation_physics)
185 32480 : dynamic_cast<DifferentiableSystem &>(system).pop_physics();
186 30825 : }
187 : } // End if QoI in set and flux/dirichlet boundary QoI
188 : } // End loop over QoIs
189 :
190 : // We'll want to back up all coarse grid vectors
191 952 : std::map<std::string, std::unique_ptr<NumericVector<Number>>> coarse_vectors;
192 180983 : for (const auto & [var_name, vec] : as_range(system.vectors_begin(), system.vectors_end()))
193 323940 : coarse_vectors[var_name] = vec->clone();
194 :
195 : // Back up the coarse solution and coarse local solution
196 17142 : std::unique_ptr<NumericVector<Number>> coarse_solution = system.solution->clone();
197 17142 : std::unique_ptr<NumericVector<Number>> coarse_local_solution = system.current_local_solution->clone();
198 :
199 : // We need to make sure that the coarse adjoint vectors used in the
200 : // calculations below are preserved during reinit, regardless of how
201 : // the user is treating them in their code
202 : // The adjoint lift function we have defined above is set to be preserved
203 : // by default
204 17142 : std::vector<bool> old_adjoints_projection_settings(system.n_qois());
205 49785 : for (auto j : make_range(system.n_qois()))
206 : {
207 33119 : if (_qoi_set.has_index(j))
208 : {
209 : // Get the vector preservation setting for this adjoint vector
210 34065 : auto adjoint_vector_name = system.vector_name(system.get_adjoint_solution(j));
211 33119 : auto old_adjoint_vector_projection_setting = system.vector_preservation(adjoint_vector_name);
212 :
213 : // Save for restoration later on
214 33119 : old_adjoints_projection_settings[j] = old_adjoint_vector_projection_setting;
215 :
216 : // Set the preservation to true for the upcoming reinits
217 33119 : system.set_vector_preservation(adjoint_vector_name, true);
218 : }
219 : }
220 :
221 : // And we'll need to temporarily change solution projection settings
222 : bool old_projection_setting;
223 16666 : old_projection_setting = system.project_solution_on_reinit();
224 :
225 : // Make sure the solution is projected when we refine the mesh
226 16666 : system.project_solution_on_reinit() = true;
227 :
228 : // And we need to make sure we dont reapply constraints after refining the mesh
229 : bool old_project_with_constraints_setting;
230 952 : old_project_with_constraints_setting = system.get_project_with_constraints();
231 :
232 476 : system.set_project_with_constraints(false);
233 :
234 : // And it'll be best to avoid any repartitioning
235 17142 : std::unique_ptr<Partitioner> old_partitioner(mesh.partitioner().release());
236 :
237 : // And we can't allow any renumbering
238 952 : const bool old_renumbering_setting = mesh.allow_renumbering();
239 476 : mesh.allow_renumbering(false);
240 :
241 : // Use a non-standard solution vector if necessary
242 16666 : if (solution_vector && solution_vector != system.solution.get())
243 : {
244 0 : NumericVector<Number> * newsol =
245 : const_cast<NumericVector<Number> *> (solution_vector);
246 0 : newsol->swap(*system.solution);
247 0 : system.update();
248 : }
249 :
250 : // Resize the error_per_cell vector to be
251 : // the number of elements, initialized to 0.
252 16666 : error_per_cell.clear();
253 16666 : error_per_cell.resize (mesh.max_elem_id(), 0.);
254 :
255 : #ifndef NDEBUG
256 : // These variables are only used in assertions later so
257 : // avoid declaring them unless asserts are active.
258 476 : const dof_id_type n_coarse_elem = mesh.n_active_elem();
259 :
260 476 : dof_id_type local_dof_bearing_nodes = 0;
261 476 : const unsigned int sysnum = system.number();
262 5146 : for (const auto * node : mesh.local_node_ptr_range())
263 4670 : for (unsigned int v=0, nvars = node->n_vars(sysnum);
264 4670 : v != nvars; ++v)
265 4670 : if (node->n_comp(sysnum, v))
266 : {
267 4670 : local_dof_bearing_nodes++;
268 4670 : break;
269 : }
270 : #endif // NDEBUG
271 :
272 : // Uniformly refine the mesh
273 17618 : MeshRefinement mesh_refinement(mesh);
274 :
275 : // We only need to worry about Galerkin orthogonality if we
276 : // are estimating discretization error in a single model setting
277 : {
278 476 : const bool swapping_adjoint_physics = _adjoint_evaluation_physics;
279 476 : if(!swapping_adjoint_physics)
280 12 : libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
281 : }
282 :
283 : // FIXME: this may break if there is more than one System
284 : // on this mesh but estimate_error was still called instead of
285 : // estimate_errors
286 17518 : for (unsigned int i = 0; i != number_h_refinements; ++i)
287 : {
288 852 : mesh_refinement.uniformly_refine(1);
289 852 : es.reinit();
290 : }
291 :
292 16666 : for (unsigned int i = 0; i != number_p_refinements; ++i)
293 : {
294 0 : mesh_refinement.uniformly_p_refine(1);
295 0 : es.reinit();
296 : }
297 :
298 : // Copy the projected coarse grid solutions, which will be
299 : // overwritten by solve()
300 1428 : std::vector<std::unique_ptr<NumericVector<Number>>> coarse_adjoints;
301 49785 : for (auto j : make_range(system.n_qois()))
302 : {
303 33119 : if (_qoi_set.has_index(j))
304 : {
305 34065 : auto coarse_adjoint = NumericVector<Number>::build(mesh.comm());
306 :
307 : // Can do "fast" init since we're overwriting this in a sec
308 34065 : coarse_adjoint->init(system.get_adjoint_solution(j),
309 1892 : /* fast = */ true);
310 :
311 33119 : *coarse_adjoint = system.get_adjoint_solution(j);
312 :
313 33119 : coarse_adjoints.emplace_back(std::move(coarse_adjoint));
314 31227 : }
315 : else
316 0 : coarse_adjoints.emplace_back(nullptr);
317 : }
318 :
319 : // Next, we are going to build up the residual for evaluating the
320 : // error estimate.
321 :
322 : // If the residual physics pointer is not null, use it when
323 : // evaluating here.
324 : {
325 16666 : if (_residual_evaluation_physics)
326 16240 : dynamic_cast<DifferentiableSystem &>(system).push_physics(*_residual_evaluation_physics);
327 :
328 : // Rebuild the rhs with the projected primal solution, do not apply constraints
329 16666 : system.assembly(true, false, false, true);
330 :
331 : // Restore the original physics if needed
332 16666 : if (_residual_evaluation_physics)
333 16240 : dynamic_cast<DifferentiableSystem &>(system).pop_physics();
334 : }
335 :
336 16666 : NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem &>(system)).get_vector("RHS Vector");
337 16666 : projected_residual.close();
338 :
339 16666 : if (_adjoint_evaluation_physics)
340 16240 : dynamic_cast<DifferentiableSystem &>(system).push_physics(*_adjoint_evaluation_physics);
341 :
342 : // Solve the adjoint problem(s) on the refined FE space
343 : // The matrix will be reassembled again because we have refined the mesh
344 : // If we have no h or p refinements, no need to solve for a fine adjoint
345 16666 : if(number_h_refinements > 0 || number_p_refinements > 0)
346 426 : system.adjoint_solve(_qoi_set);
347 :
348 : // Swap back if needed, recall that _adjoint_evaluation_physics now holds the pointer
349 : // to the pre-swap physics
350 16666 : if (_adjoint_evaluation_physics)
351 16240 : dynamic_cast<DifferentiableSystem &>(system).pop_physics();
352 :
353 : // Now that we have the refined adjoint solution and the projected primal solution,
354 : // we first compute the global QoI error estimate
355 :
356 : // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
357 17142 : computed_global_QoI_errors.resize(system.n_qois());
358 :
359 : // Loop over all the adjoint solutions and get the QoI error
360 : // contributions from all of them. While we're looping anyway we'll
361 : // pull off the coarse adjoints
362 49785 : for (auto j : make_range(system.n_qois()))
363 : {
364 : // Skip this QoI if not in the QoI Set
365 33119 : if (_qoi_set.has_index(j))
366 : {
367 :
368 : // If the adjoint solution has heterogeneous dirichlet
369 : // values, then to get a proper error estimate here we need
370 : // to subtract off a coarse grid lift function.
371 : // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [L]+)| + HOT
372 33119 : if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
373 : {
374 : // Need to create a string with current loop index to retrieve
375 : // the correct vector from the liftvectors map
376 32693 : std::ostringstream liftfunc_name;
377 31759 : liftfunc_name << "adjoint_lift_function" << j;
378 :
379 : // Subtract off the corresponding lift vector from the adjoint solution
380 32693 : system.get_adjoint_solution(j) -= system.get_vector(liftfunc_name.str());
381 :
382 : // Now evaluate R(u^h, z^h+ - lift)
383 32693 : computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
384 :
385 : // Add the lift back to get back the adjoint
386 33627 : system.get_adjoint_solution(j) += system.get_vector(liftfunc_name.str());
387 :
388 30825 : }
389 : else
390 : {
391 : // Usual dual weighted residual error estimate
392 : // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+)| + HOT
393 426 : computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
394 : }
395 :
396 : }
397 : }
398 :
399 : // We are all done with Dirichlet lift vectors and they should be removed, lest we run into I/O issues later
400 49785 : for (auto j : make_range(system.n_qois()))
401 : {
402 : // Skip this QoI if not in the QoI Set
403 33119 : if (_qoi_set.has_index(j))
404 : {
405 : // Lifts are created only for adjoint dirichlet QoIs
406 33119 : if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j))
407 : {
408 : // Need to create a string with current loop index to retrieve
409 : // the correct vector from the liftvectors map
410 32693 : std::ostringstream liftfunc_name;
411 31759 : liftfunc_name << "adjoint_lift_function" << j;
412 :
413 : // Remove the lift vector from system since we did not write it to file and it cannot be retrieved
414 33627 : system.remove_vector(liftfunc_name.str());
415 30825 : }
416 : }
417 : }
418 :
419 : // Done with the global error estimates, now construct the element wise error indicators
420 :
421 : // To get a better element wise breakdown of the error estimate,
422 : // we subtract off a coarse representation of the adjoint solution.
423 : // |Q(u) - Q(u^h)| = |R([u^h]+, z^h+ - [z^h]+)|
424 : // This remains valid for all combinations of heterogenous adjoint bcs and
425 : // stabilized/non-stabilized formulations, except for the case where we not using a
426 : // heterogenous adjoint bc and have a stabilized formulation.
427 : // Then, R(u^h_s, z^h_s) != 0 (no Galerkin orthogonality w.r.t the non-stabilized residual)
428 49785 : for (auto j : make_range(system.n_qois()))
429 : {
430 : // Skip this QoI if not in the QoI Set
431 33119 : if (_qoi_set.has_index(j))
432 : {
433 : // If we have a nullptr residual evaluation physics pointer, we
434 : // assume the user's formulation is consistent from mesh to
435 : // mesh, so we have Galerkin orthogonality and we can get
436 : // better indicator results by subtracting a coarse adjoint.
437 :
438 : // If we have a residual evaluation physics pointer, but we
439 : // also have heterogeneous adjoint dirichlet boundaries,
440 : // then we have to subtract off *some* lift function for
441 : // consistency, and we choose the coarse adjoint in lieu of
442 : // having any better ideas.
443 :
444 : // If we have a residual evaluation physics pointer and we
445 : // have homogeneous adjoint dirichlet boundaries, then we
446 : // don't have to subtract off anything, and with stabilized
447 : // formulations we get the best results if we don't.
448 33119 : if(system.get_dof_map().has_adjoint_dirichlet_boundaries(j)
449 33119 : || !_residual_evaluation_physics)
450 : {
451 : // z^h+ -> z^h+ - [z^h]+
452 33119 : system.get_adjoint_solution(j) -= *coarse_adjoints[j];
453 : }
454 : }
455 : }
456 :
457 : // We ought to account for 'spill-over' effects while computing the
458 : // element error indicators This happens because the same dof is
459 : // shared by multiple elements, one way of mitigating this is to
460 : // scale the contribution from each dof by the number of elements it
461 : // belongs to We first obtain the number of elements each node
462 : // belongs to
463 :
464 : // A map that relates a node id to an int that will tell us how many elements it is a node of
465 952 : std::unordered_map<dof_id_type, unsigned int> shared_element_count;
466 :
467 : // To fill this map, we will loop over elements, and then in each element, we will loop
468 : // over the nodes each element contains, and then query it for the number of coarse
469 : // grid elements it was a node of
470 :
471 : // Keep track of which nodes we have already dealt with
472 952 : std::unordered_set<dof_id_type> processed_node_ids;
473 :
474 : // We will be iterating over all the active elements in the fine mesh that live on
475 : // this processor.
476 258392 : for (const auto & elem : mesh.active_local_element_ptr_range())
477 1190336 : for (const Node & node : elem->node_ref_range())
478 : {
479 : // Get the id of this node
480 1056896 : dof_id_type node_id = node.id();
481 :
482 : // If we haven't already processed this node, do so now
483 1056896 : if (processed_node_ids.find(node_id) == processed_node_ids.end())
484 : {
485 : // Declare a neighbor_set to be filled by the find_point_neighbors
486 83464 : std::set<const Elem *> fine_grid_neighbor_set;
487 :
488 : // Call find_point_neighbors to fill the neighbor_set
489 514112 : elem->find_point_neighbors(node, fine_grid_neighbor_set);
490 :
491 : // A vector to hold the coarse grid parents neighbors
492 83464 : std::vector<dof_id_type> coarse_grid_neighbors;
493 :
494 : // Loop over all the fine neighbors of this node
495 1654482 : for (const auto & fine_elem : fine_grid_neighbor_set)
496 : {
497 : // Find the element id for the corresponding coarse grid element
498 1140370 : const Elem * coarse_elem = fine_elem;
499 3266134 : for (unsigned int j = 0; j != number_h_refinements; ++j)
500 : {
501 171798 : libmesh_assert (coarse_elem->parent());
502 :
503 343596 : coarse_elem = coarse_elem->parent();
504 : }
505 :
506 : // Loop over the existing coarse neighbors and check if this one is
507 : // already in there
508 1140370 : const dof_id_type coarse_id = coarse_elem->id();
509 91467 : std::size_t j = 0;
510 :
511 : // If the set already contains this element break out of the loop
512 1401377 : for (std::size_t cgns = coarse_grid_neighbors.size(); j != cgns; j++)
513 717061 : if (coarse_grid_neighbors[j] == coarse_id)
514 37350 : break;
515 :
516 : // If we didn't leave the loop even at the last element,
517 : // this is a new neighbour, put in the coarse_grid_neighbor_set
518 1140370 : if (j == coarse_grid_neighbors.size())
519 684316 : coarse_grid_neighbors.push_back(coarse_id);
520 : } // End loop over fine neighbors
521 :
522 : // Set the shared_neighbour index for this node to the
523 : // size of the coarse grid neighbor set
524 555844 : shared_element_count[node_id] =
525 83464 : cast_int<unsigned int>(coarse_grid_neighbors.size());
526 :
527 : // Add this node to processed_node_ids vector
528 41732 : processed_node_ids.insert(node_id);
529 : } // End if not processed node
530 15714 : } // End loop over nodes
531 :
532 : // Get a DoF map, will be used to get the nodal dof_indices for each element
533 476 : DofMap & dof_map = system.get_dof_map();
534 :
535 : // The global DOF indices, we will use these later on when we compute the element wise indicators
536 952 : std::vector<dof_id_type> dof_indices;
537 :
538 : // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processors) onto a
539 : // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
540 : // an element it owns
541 17142 : std::unique_ptr<NumericVector<Number>> localized_projected_residual = NumericVector<Number>::build(system.comm());
542 17142 : localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
543 17142 : projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
544 :
545 : // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
546 17142 : std::unique_ptr<NumericVector<Number>> localized_adjoint_solution = NumericVector<Number>::build(system.comm());
547 17142 : localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
548 :
549 : // We will loop over each adjoint solution, localize that adjoint
550 : // solution and then loop over local elements
551 49785 : for (auto i : make_range(system.n_qois()))
552 : {
553 : // Skip this QoI if not in the QoI Set
554 33119 : if (_qoi_set.has_index(i))
555 : {
556 : // Get the weight for the current QoI
557 946 : Real error_weight = _qoi_set.weight(i);
558 :
559 33119 : (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
560 :
561 : // Loop over elements
562 398092 : for (const auto & elem : mesh.active_local_element_ptr_range())
563 : {
564 : // Go up number_h_refinements levels up to find the coarse parent
565 181696 : const Elem * coarse = elem;
566 :
567 504256 : for (unsigned int j = 0; j != number_h_refinements; ++j)
568 : {
569 26880 : libmesh_assert (coarse->parent());
570 :
571 53760 : coarse = coarse->parent();
572 : }
573 :
574 30592 : const dof_id_type e_id = coarse->id();
575 :
576 : // Get the local to global degree of freedom maps for this element
577 181696 : dof_map.dof_indices (elem, dof_indices);
578 :
579 : // We will have to manually do the dot products.
580 28736 : Number local_contribution = 0.;
581 :
582 : // Sum the contribution to the error indicator for each element from the current QoI
583 1714880 : for (const auto & dof : dof_indices)
584 1533184 : local_contribution += (*localized_projected_residual)(dof) * (*localized_adjoint_solution)(dof);
585 :
586 : // Multiply by the error weight for this QoI
587 168256 : local_contribution *= error_weight;
588 :
589 : // FIXME: we're throwing away information in the
590 : // --enable-complex case
591 196992 : error_per_cell[e_id] += static_cast<ErrorVectorReal>
592 15296 : (std::abs(local_contribution));
593 :
594 31227 : } // End loop over elements
595 : } // End if belong to QoI set
596 : } // End loop over QoIs
597 :
598 : // Don't bother projecting the solution; we'll restore from backup
599 : // after coarsening
600 16666 : system.project_solution_on_reinit() = false;
601 :
602 : // Uniformly coarsen the mesh, without projecting the solution
603 : // Only need to do this if we are estimating discretization error
604 : // with a single physics residual
605 476 : libmesh_assert (_adjoint_evaluation_physics ||
606 : number_h_refinements > 0 || number_p_refinements > 0);
607 :
608 17518 : for (unsigned int i = 0; i != number_h_refinements; ++i)
609 : {
610 852 : mesh_refinement.uniformly_coarsen(1);
611 : // FIXME - should the reinits here be necessary? - RHS
612 852 : es.reinit();
613 : }
614 :
615 16666 : for (unsigned int i = 0; i != number_p_refinements; ++i)
616 : {
617 0 : mesh_refinement.uniformly_p_coarsen(1);
618 0 : es.reinit();
619 : }
620 :
621 : // We should have the same number of active elements as when we started,
622 : // but not necessarily the same number of elements since reinit() doesn't
623 : // always call contract()
624 476 : libmesh_assert_equal_to (n_coarse_elem, mesh.n_active_elem());
625 :
626 : // We should have the same number of dof-bearing nodes as when we
627 : // started
628 : #ifndef NDEBUG
629 476 : dof_id_type final_local_dof_bearing_nodes = 0;
630 5146 : for (const auto * node : mesh.local_node_ptr_range())
631 4670 : for (unsigned int v=0, nvars = node->n_vars(sysnum);
632 4670 : v != nvars; ++v)
633 4670 : if (node->n_comp(sysnum, v))
634 : {
635 4670 : final_local_dof_bearing_nodes++;
636 4670 : break;
637 : }
638 476 : libmesh_assert_equal_to (local_dof_bearing_nodes,
639 : final_local_dof_bearing_nodes);
640 : #endif // NDEBUG
641 :
642 : // Restore old solutions and clean up the heap
643 16666 : system.project_solution_on_reinit() = old_projection_setting;
644 476 : system.set_project_with_constraints(old_project_with_constraints_setting);
645 :
646 : // Restore the adjoint vector preservation settings
647 49785 : for (auto j : make_range(system.n_qois()))
648 : {
649 33119 : if (_qoi_set.has_index(j))
650 : {
651 33119 : auto adjoint_vector_name = system.vector_name(system.get_adjoint_solution(j));
652 34065 : system.set_vector_preservation(adjoint_vector_name, old_adjoints_projection_settings[j]);
653 : }
654 : }
655 :
656 : // Restore the coarse solution vectors and delete their copies
657 17142 : *system.solution = *coarse_solution;
658 17142 : *system.current_local_solution = *coarse_local_solution;
659 :
660 148290 : for (const auto & pr : as_range(system.vectors_begin(), system.vectors_end()))
661 : {
662 : // The (string) name of this vector
663 131624 : const std::string & var_name = pr.first;
664 :
665 : // If it's a vector we already had (and not a newly created
666 : // vector like an adjoint rhs), we need to restore it.
667 131624 : if (auto it = coarse_vectors.find(var_name);
668 3760 : it != coarse_vectors.end())
669 : {
670 3760 : NumericVector<Number> * coarsevec = it->second.get();
671 131624 : system.get_vector(var_name) = *coarsevec;
672 :
673 131624 : coarsevec->clear();
674 : }
675 : }
676 :
677 : // Restore old partitioner and renumbering settings
678 16666 : mesh.partitioner().reset(old_partitioner.release());
679 476 : mesh.allow_renumbering(old_renumbering_setting);
680 :
681 : // Finally sum the vector of estimated error values.
682 16666 : this->reduce_error(error_per_cell, system.comm());
683 :
684 : // We don't take a square root here; this is a goal-oriented
685 : // estimate not a Hilbert norm estimate.
686 63808 : } // end estimate_error function
687 :
688 : }// namespace libMesh
689 :
690 : #endif // #ifdef LIBMESH_ENABLE_AMR
|
