adjoint_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
25 // Local Includes
26 #include "libmesh/dof_map.h"
27 #include "libmesh/elem.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/fe.h"
31 #include "libmesh/fe_interface.h"
32 #include "libmesh/libmesh_common.h"
34 #include "libmesh/mesh_base.h"
36 #include "libmesh/numeric_vector.h"
37 #include "libmesh/quadrature.h"
38 #include "libmesh/system.h"
40 #include "libmesh/partitioner.h"
42 
43 #include LIBMESH_INCLUDE_UNORDERED_MAP
44 #include LIBMESH_INCLUDE_UNORDERED_SET
45 
46 #ifdef LIBMESH_ENABLE_AMR
47 
48 namespace libMesh
49 {
50 
51 //-----------------------------------------------------------------
52 // AdjointRefinementErrorEstimator
53 
54 // As of 10/2/2012, this function implements a 'brute-force' adjoint based QoI
55 // error estimator, using the following relationship:
56 // Q(u) - Q(u_h) \aprrox - R( (u_h)_(h/2), z_(h/2) ) .
57 // where: Q(u) is the true QoI
58 // u_h is the approximate primal solution on the current FE space
59 // Q(u_h) is the approximate QoI
60 // z_(h/2) is the adjoint corresponding to Q, on a richer FE space
61 // (u_h)_(h/2) is a projection of the primal solution on the richer FE space
62 // By richer FE space, we mean a grid that has been refined once and a polynomial order
63 // that has been increased once, i.e. one h and one p refinement
64 
65 // Both a global QoI error estimate and element wise error indicators are included
66 // Note that the element wise error indicators slightly over estimate the error in
67 // each element
68 
70  ErrorVector& error_per_cell,
71  const NumericVector<Number>* solution_vector,
72  bool /*estimate_parent_error*/)
73 {
74  // We have to break the rules here, because we can't refine a const System
75  System& system = const_cast<System&>(_system);
76 
77  // An EquationSystems reference will be convenient.
79 
80  // The current mesh
81  MeshBase& mesh = es.get_mesh();
82 
83  // Resize the error_per_cell vector to be
84  // the number of elements, initialized to 0.
85  error_per_cell.clear();
86  error_per_cell.resize (mesh.max_elem_id(), 0.);
87 
88  // We'll want to back up all coarse grid vectors
89  std::map<std::string, NumericVector<Number> *> coarse_vectors;
90  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
91  system.vectors_end(); ++vec)
92  {
93  // The (string) name of this vector
94  const std::string& var_name = vec->first;
95 
96  coarse_vectors[var_name] = vec->second->clone().release();
97  }
98  // Back up the coarse solution and coarse local solution
99  NumericVector<Number> * coarse_solution =
100  system.solution->clone().release();
101  NumericVector<Number> * coarse_local_solution =
102  system.current_local_solution->clone().release();
103  // And make copies of the projected solution
104  NumericVector<Number> * projected_solution;
105 
106  // And we'll need to temporarily change solution projection settings
107  bool old_projection_setting;
108  old_projection_setting = system.project_solution_on_reinit();
109 
110  // Make sure the solution is projected when we refine the mesh
111  system.project_solution_on_reinit() = true;
112 
113  // And it'll be best to avoid any repartitioning
114  AutoPtr<Partitioner> old_partitioner = mesh.partitioner();
115  mesh.partitioner().reset(NULL);
116 
117  // And we can't allow any renumbering
118  const bool old_renumbering_setting = mesh.allow_renumbering();
119  mesh.allow_renumbering(false);
120 
121  // Use a non-standard solution vector if necessary
122  if (solution_vector && solution_vector != system.solution.get())
123  {
124  NumericVector<Number> *newsol =
125  const_cast<NumericVector<Number>*> (solution_vector);
126  newsol->swap(*system.solution);
127  system.update();
128  }
129 
130  // Loop over all the adjoint problems and, if any have heterogenous
131  // Dirichlet conditions, get the corresponding coarse lift
132  // function(s)
133  for (unsigned int j=0; j != system.qoi.size(); j++)
134  {
135  // Skip this QoI if it is not in the QoI Set or if there are no
136  // heterogeneous Dirichlet boundaries for it
137  if (_qoi_set.has_index(j) &&
139  {
140  std::ostringstream liftfunc_name;
141  liftfunc_name << "adjoint_lift_function" << j;
142  NumericVector<Number> &liftvec =
143  system.add_vector(liftfunc_name.str());
144 
146  (system, &liftvec, true);
147  }
148  }
149 
150 
151 
152 
153 #ifndef NDEBUG
154  // n_coarse_elem is only used in an assertion later so
155  // avoid declaring it unless asserts are active.
156  const dof_id_type n_coarse_elem = mesh.n_elem();
157 #endif
158 
159  // Uniformly refine the mesh
160  MeshRefinement mesh_refinement(mesh);
161 
163 
164  // FIXME: this may break if there is more than one System
165  // on this mesh but estimate_error was still called instead of
166  // estimate_errors
167  for (unsigned int i = 0; i != number_h_refinements; ++i)
168  {
169  mesh_refinement.uniformly_refine(1);
170  es.reinit();
171  }
172 
173  for (unsigned int i = 0; i != number_p_refinements; ++i)
174  {
175  mesh_refinement.uniformly_p_refine(1);
176  es.reinit();
177  }
178 
179  // Copy the projected coarse grid solutions, which will be
180  // overwritten by solve()
181  projected_solution = NumericVector<Number>::build(mesh.comm()).release();
182  projected_solution->init(system.solution->size(), true, SERIAL);
183  system.solution->localize(*projected_solution,
184  system.get_dof_map().get_send_list());
185 
186  // Rebuild the rhs with the projected primal solution
187  (dynamic_cast<ImplicitSystem&>(system)).assembly(true, false);
188  NumericVector<Number> & projected_residual = (dynamic_cast<ExplicitSystem&>(system)).get_vector("RHS Vector");
189  projected_residual.close();
190 
191  // Solve the adjoint problem(s) on the refined FE space
192  system.adjoint_solve(_qoi_set);
193 
194  // Now that we have the refined adjoint solution and the projected primal solution,
195  // we first compute the global QoI error estimate
196 
197  // Resize the computed_global_QoI_errors vector to hold the error estimates for each QoI
198  computed_global_QoI_errors.resize(system.qoi.size());
199 
200  // Loop over all the adjoint solutions and get the QoI error
201  // contributions from all of them. While we're looping anyway we'll
202  // handle heterogenous adjoint BCs
203  for (unsigned int j=0; j != system.qoi.size(); j++)
204  {
205  // Skip this QoI if not in the QoI Set
206  if (_qoi_set.has_index(j))
207  {
208  // If the adjoint solution has heterogeneous dirichlet
209  // values, then to get a proper error estimate here we need
210  // to subtract off a coarse grid lift function. We won't
211  // need the lift function afterward. We won't even need the
212  // fine adjoint solution afterward, so we'll modify it here.
214  {
215  std::ostringstream liftfunc_name;
216  liftfunc_name << "adjoint_lift_function" << j;
217  system.get_adjoint_solution(j) -=
218  system.get_vector(liftfunc_name.str());
219 
220  system.remove_vector(liftfunc_name.str());
221  }
222 
223  computed_global_QoI_errors[j] = projected_residual.dot(system.get_adjoint_solution(j));
224  }
225  }
226 
227  // Done with the global error estimates, now construct the element wise error indicators
228 
229  // We ought to account for 'spill-over' effects while computing the
230  // element error indicators This happens because the same dof is
231  // shared by multiple elements, one way of mitigating this is to
232  // scale the contribution from each dof by the number of elements it
233  // belongs to We first obtain the number of elements each node
234  // belongs to
235 
236  // A map that relates a node id to an int that will tell us how many elements it is a node of
237  LIBMESH_BEST_UNORDERED_MAP<dof_id_type, unsigned int>shared_element_count;
238 
239  // To fill this map, we will loop over elements, and then in each element, we will loop
240  // over the nodes each element contains, and then query it for the number of coarse
241  // grid elements it was a node of
242 
243  // We will be iterating over all the active elements in the fine mesh that live on
244  // this processor
247 
248  // Keep track of which nodes we have already dealt with
249  LIBMESH_BEST_UNORDERED_SET<dof_id_type> processed_node_ids;
250 
251  // Start loop over elems
252  for(; elem_it != elem_end; ++elem_it)
253  {
254  // Pointer to this element
255  const Elem* elem = *elem_it;
256 
257  // Loop over the nodes in the element
258  for(unsigned int n=0; n != elem->n_nodes(); ++n)
259  {
260  // Get a pointer to the current node
261  Node* node = elem->get_node(n);
262 
263  // Get the id of this node
264  dof_id_type node_id = node->id();
265 
266  // If we havent already processed this node, do so now
267  if(processed_node_ids.find(node_id) == processed_node_ids.end())
268  {
269  // Declare a neighbor_set to be filled by the find_point_neighbors
270  std::set<const Elem *> fine_grid_neighbor_set;
271 
272  // Call find_point_neighbors to fill the neighbor_set
273  elem->find_point_neighbors(*node, fine_grid_neighbor_set);
274 
275  // A vector to hold the coarse grid parents neighbors
276  std::vector<dof_id_type> coarse_grid_neighbors;
277 
278  // Iterators over the fine grid neighbors set
279  std::set<const Elem*>::iterator fine_neighbor_it = fine_grid_neighbor_set.begin();
280  const std::set<const Elem*>::iterator fine_neighbor_end = fine_grid_neighbor_set.end();
281 
282  // Loop over all the fine neighbors of this node
283  for(; fine_neighbor_it != fine_neighbor_end ; ++fine_neighbor_it)
284  {
285  // Pointer to the current fine neighbor element
286  const Elem* fine_elem = *fine_neighbor_it;
287 
288  // Find the element id for the corresponding coarse grid element
289  const Elem* coarse_elem = fine_elem;
290  for (unsigned int j = 0; j != number_h_refinements; ++j)
291  {
292  libmesh_assert (coarse_elem->parent());
293 
294  coarse_elem = coarse_elem->parent();
295  }
296 
297  // Loop over the existing coarse neighbors and check if this one is
298  // already in there
299  const dof_id_type coarse_id = coarse_elem->id();
300  std::size_t j = 0;
301  for (; j != coarse_grid_neighbors.size(); j++)
302  {
303  // If the set already contains this element break out of the loop
304  if(coarse_grid_neighbors[j] == coarse_id)
305  {
306  break;
307  }
308  }
309 
310  // If we didn't leave the loop even at the last element,
311  // this is a new neighbour, put in the coarse_grid_neighbor_set
312  if(j == coarse_grid_neighbors.size())
313  {
314  coarse_grid_neighbors.push_back(coarse_id);
315  }
316 
317  } // End loop over fine neighbors
318 
319  // Set the shared_neighbour index for this node to the
320  // size of the coarse grid neighbor set
321  shared_element_count[node_id] =
322  libmesh_cast_int<unsigned int>(coarse_grid_neighbors.size());
323 
324  // Add this node to processed_node_ids vector
325  processed_node_ids.insert(node_id);
326 
327  } // End if not processed node
328 
329  } // End loop over nodes
330 
331  } // End loop over elems
332 
333  // Get a DoF map, will be used to get the nodal dof_indices for each element
334  DofMap &dof_map = system.get_dof_map();
335 
336  // The global DOF indices, we will use these later on when we compute the element wise indicators
337  std::vector<dof_id_type> dof_indices;
338 
339  // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processsors) onto a
340  // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
341  // an element it owns
342  AutoPtr<NumericVector<Number> > localized_projected_residual = NumericVector<Number>::build(system.comm());
343  localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
344  projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
345 
346  // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
347  AutoPtr<NumericVector<Number> > localized_adjoint_solution = NumericVector<Number>::build(system.comm());
348  localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
349 
350  // We will loop over each adjoint solution, localize that adjoint
351  // solution and then loop over local elements
352  for (unsigned int i=0; i != system.qoi.size(); i++)
353  {
354  // Skip this QoI if not in the QoI Set
355  if (_qoi_set.has_index(i))
356  {
357  // Get the weight for the current QoI
358  Real error_weight = _qoi_set.weight(i);
359 
360  (system.get_adjoint_solution(i)).localize(*localized_adjoint_solution, system.get_dof_map().get_send_list());
361 
362  // Loop over elements
365 
366  for(; elem_it != elem_end; ++elem_it)
367  {
368  // Pointer to the element
369  const Elem* elem = *elem_it;
370 
371  // Go up number_h_refinements levels up to find the coarse parent
372  const Elem* coarse = elem;
373 
374  for (unsigned int j = 0; j != number_h_refinements; ++j)
375  {
376  libmesh_assert (coarse->parent());
377 
378  coarse = coarse->parent();
379  }
380 
381  const dof_id_type e_id = coarse->id();
382 
383  // Get the local to global degree of freedom maps for this element
384  dof_map.dof_indices (elem, dof_indices);
385 
386  // We will have to manually do the dot products.
387  Number local_contribution = 0.;
388 
389  for (unsigned int j=0; j != dof_indices.size(); j++)
390  {
391  // The contribution to the error indicator for this element from the current QoI
392  local_contribution += (*localized_projected_residual)(dof_indices[j]) * (*localized_adjoint_solution)(dof_indices[j]);
393  }
394 
395  // Multiply by the error weight for this QoI
396  local_contribution *= error_weight;
397 
398  // FIXME: we're throwing away information in the
399  // --enable-complex case
400  error_per_cell[e_id] += static_cast<ErrorVectorReal>
401  (libmesh_real(local_contribution));
402 
403  } // End loop over elements
404 
405  } // End if belong to QoI set
406 
407  } // End loop over QoIs
408 
409  // Don't bother projecting the solution; we'll restore from backup
410  // after coarsening
411  system.project_solution_on_reinit() = false;
412 
413  // Uniformly coarsen the mesh, without projecting the solution
415 
416  for (unsigned int i = 0; i != number_h_refinements; ++i)
417  {
418  mesh_refinement.uniformly_coarsen(1);
419  // FIXME - should the reinits here be necessary? - RHS
420  es.reinit();
421  }
422 
423  for (unsigned int i = 0; i != number_p_refinements; ++i)
424  {
425  mesh_refinement.uniformly_p_coarsen(1);
426  es.reinit();
427  }
428 
429  // We should be back where we started
430  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
431 
432  // Restore old solutions and clean up the heap
433  system.project_solution_on_reinit() = old_projection_setting;
434 
435  // Restore the coarse solution vectors and delete their copies
436  *system.solution = *coarse_solution;
437  delete coarse_solution;
438  *system.current_local_solution = *coarse_local_solution;
439  delete coarse_local_solution;
440  delete projected_solution;
441 
442  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
443  system.vectors_end(); ++vec)
444  {
445  // The (string) name of this vector
446  const std::string& var_name = vec->first;
447 
448  // If it's a vector we already had (and not a newly created
449  // vector like an adjoint rhs), we need to restore it.
450  std::map<std::string, NumericVector<Number> *>::iterator it =
451  coarse_vectors.find(var_name);
452  if (it != coarse_vectors.end())
453  {
454  NumericVector<Number> *coarsevec = it->second;
455  system.get_vector(var_name) = *coarsevec;
456 
457  coarsevec->clear();
458  delete coarsevec;
459  }
460  }
461 
462  // Restore old partitioner and renumbering settings
463  mesh.partitioner() = old_partitioner;
464  mesh.allow_renumbering(old_renumbering_setting);
465 
466  // Fiinally sum the vector of estimated error values.
467  this->reduce_error(error_per_cell, system.comm());
468 
469  // We don't take a square root here; this is a goal-oriented
470  // estimate not a Hilbert norm estimate.
471 } // end estimate_error function
472 
473 }// namespace libMesh
474 
475 #endif // #ifdef LIBMESH_ENABLE_AMR

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo