adjoint_residual_error_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 
19 // C++ includes
20 #include <iostream>
21 #include <iomanip>
22 #include <sstream>
23 
24 // Local Includes
26 #include "libmesh/error_vector.h"
29 #include "libmesh/numeric_vector.h"
30 #include "libmesh/system.h"
31 #include "libmesh/system_norm.h"
32 #include "libmesh/qoi_set.h"
33 
34 
35 namespace libMesh
36 {
37 
38 //-----------------------------------------------------------------
39 // AdjointResidualErrorEstimator implementations
42  error_plot_suffix(),
43  _primal_error_estimator(new PatchRecoveryErrorEstimator()),
44  _dual_error_estimator(new PatchRecoveryErrorEstimator()),
45  _qoi_set(QoISet())
46 {
47 }
48 
49 
50 
52  ErrorVector& error_per_cell,
53  const NumericVector<Number>* solution_vector,
54  bool estimate_parent_error)
55 {
56  START_LOG("estimate_error()", "AdjointResidualErrorEstimator");
57 
58  // The current mesh
59  const MeshBase& mesh = _system.get_mesh();
60 
61  // Resize the error_per_cell vector to be
62  // the number of elements, initialize it to 0.
63  error_per_cell.resize (mesh.max_elem_id());
64  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
65 
66  // Get the number of variables in the system
67  unsigned int n_vars = _system.n_vars();
68 
69  // We need to make a map of the pointer to the solution vector
70  std::map<const System*, const NumericVector<Number>*>solutionvecs;
71  solutionvecs[&_system] = _system.solution.get();
72 
73  // Solve the dual problem if we have to
74  if (!_system.is_adjoint_already_solved())
75  {
76  // FIXME - we'll need to change a lot of APIs to make this trick
77  // work with a const System...
78  System& system = const_cast<System&>(_system);
79  system.adjoint_solve(_qoi_set);
80  }
81 
82  // Flag to check whether we have not been asked to weight the variable error contributions in any specific manner
83  bool error_norm_is_identity = error_norm.is_identity();
84 
85  // Create an ErrorMap/ErrorVector to store the primal, dual and total_dual variable errors
86  ErrorMap primal_errors_per_cell;
87  ErrorMap dual_errors_per_cell;
88  ErrorMap total_dual_errors_per_cell;
89  // Allocate ErrorVectors to this map if we're going to use it
90  if (!error_norm_is_identity)
91  for(unsigned int v = 0; v < n_vars; v++)
92  {
93  primal_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
94  dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
95  total_dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
96  }
97  ErrorVector primal_error_per_cell;
98  ErrorVector dual_error_per_cell;
99  ErrorVector total_dual_error_per_cell;
100 
101  // Have we been asked to weight the variable error contributions in any specific manner
102  if(!error_norm_is_identity) // If we do
103  {
104  // Estimate the primal problem error for each variable
105  _primal_error_estimator->estimate_errors
106  (_system.get_equation_systems(), primal_errors_per_cell, &solutionvecs, estimate_parent_error);
107  }
108  else // If not
109  {
110  // Just get the combined error estimate
111  _primal_error_estimator->estimate_error
112  (_system, primal_error_per_cell, solution_vector, estimate_parent_error);
113  }
114 
115  // Sum and weight the dual error estimate based on our QoISet
116  for (unsigned int i = 0; i != _system.qoi.size(); ++i)
117  {
118  if (_qoi_set.has_index(i))
119  {
120  // Get the weight for the current QoI
121  Real error_weight = _qoi_set.weight(i);
122 
123  // We need to make a map of the pointer to the adjoint solution vector
124  std::map<const System*, const NumericVector<Number>*>adjointsolutionvecs;
125  adjointsolutionvecs[&_system] = &_system.get_adjoint_solution(i);
126 
127  // Have we been asked to weight the variable error contributions in any specific manner
128  if(!error_norm_is_identity) // If we have
129  {
130  _dual_error_estimator->estimate_errors
131  (_system.get_equation_systems(), dual_errors_per_cell, &adjointsolutionvecs,
132  estimate_parent_error);
133  }
134  else // If not
135  {
136  // Just get the combined error estimate
137  _dual_error_estimator->estimate_error
138  (_system, dual_error_per_cell, &(_system.get_adjoint_solution(i)), estimate_parent_error);
139  }
140 
141  std::size_t error_size;
142 
143  // Get the size of the first ErrorMap vector; this will give us the number of elements
144  if(!error_norm_is_identity) // If in non default weights case
145  {
146  error_size = dual_errors_per_cell[std::make_pair(&_system, 0)]->size();
147  }
148  else // If in the standard default weights case
149  {
150  error_size = dual_error_per_cell.size();
151  }
152 
153  // Resize the ErrorVector(s)
154  if(!error_norm_is_identity)
155  {
156  // Loop over variables
157  for(unsigned int v = 0; v < n_vars; v++)
158  {
159  libmesh_assert(!total_dual_errors_per_cell[std::make_pair(&_system, v)]->size() ||
160  total_dual_errors_per_cell[std::make_pair(&_system, v)]->size() == error_size) ;
161  total_dual_errors_per_cell[std::make_pair(&_system, v)]->resize(error_size);
162  }
163  }
164  else
165  {
166  libmesh_assert(!total_dual_error_per_cell.size() ||
167  total_dual_error_per_cell.size() == error_size);
168  total_dual_error_per_cell.resize(error_size);
169  }
170 
171  for (std::size_t e = 0; e != error_size; ++e)
172  {
173  // Have we been asked to weight the variable error contributions in any specific manner
174  if(!error_norm_is_identity) // If we have
175  {
176  // Loop over variables
177  for(unsigned int v = 0; v < n_vars; v++)
178  {
179  // Now fill in total_dual_error ErrorMap with the weight
180  (*total_dual_errors_per_cell[std::make_pair(&_system, v)])[e] +=
181  static_cast<ErrorVectorReal>
182  (error_weight *
183  (*dual_errors_per_cell[std::make_pair(&_system, v)])[e]);
184  }
185  }
186  else // If not
187  {
188  total_dual_error_per_cell[e] +=
189  static_cast<ErrorVectorReal>(error_weight * dual_error_per_cell[e]);
190  }
191  }
192  }
193  }
194 
195  // Do some debugging plots if requested
196  if (!error_plot_suffix.empty())
197  {
198  if(!error_norm_is_identity) // If we have
199  {
200  // Loop over variables
201  for(unsigned int v = 0; v < n_vars; v++)
202  {
203  std::ostringstream primal_out;
204  std::ostringstream dual_out;
205  primal_out << "primal_" << error_plot_suffix << ".";
206  dual_out << "dual_" << error_plot_suffix << ".";
207 
208  primal_out << std::setw(1)
209  << std::setprecision(0)
210  << std::setfill('0')
211  << std::right
212  << v;
213 
214  dual_out << std::setw(1)
215  << std::setprecision(0)
216  << std::setfill('0')
217  << std::right
218  << v;
219 
220  (*primal_errors_per_cell[std::make_pair(&_system, v)]).plot_error(primal_out.str(), _system.get_mesh());
221  (*total_dual_errors_per_cell[std::make_pair(&_system, v)]).plot_error(dual_out.str(), _system.get_mesh());
222 
223  primal_out.clear();
224  dual_out.clear();
225  }
226  }
227  else // If not
228  {
229  std::ostringstream primal_out;
230  std::ostringstream dual_out;
231  primal_out << "primal_" << error_plot_suffix ;
232  dual_out << "dual_" << error_plot_suffix ;
233 
234  primal_error_per_cell.plot_error(primal_out.str(), _system.get_mesh());
235  total_dual_error_per_cell.plot_error(dual_out.str(), _system.get_mesh());
236 
237  primal_out.clear();
238  dual_out.clear();
239  }
240  }
241 
242  // Weight the primal error by the dual error using the system norm object
243  // FIXME: we ought to thread this
244  for (unsigned int i=0; i != error_per_cell.size(); ++i)
245  {
246  // Have we been asked to weight the variable error contributions in any specific manner
247  if(!error_norm_is_identity) // If we do
248  {
249  // Create Error Vectors to pass to calculate_norm
250  std::vector<Real> cell_primal_error;
251  std::vector<Real> cell_dual_error;
252 
253  for(unsigned int v = 0; v < n_vars; v++)
254  {
255  cell_primal_error.push_back((*primal_errors_per_cell[std::make_pair(&_system, v)])[i]);
256  cell_dual_error.push_back((*total_dual_errors_per_cell[std::make_pair(&_system, v)])[i]);
257  }
258 
259  error_per_cell[i] =
260  static_cast<ErrorVectorReal>
261  (error_norm.calculate_norm(cell_primal_error, cell_dual_error));
262  }
263  else // If not
264  {
265  error_per_cell[i] = primal_error_per_cell[i]*total_dual_error_per_cell[i];
266  }
267  }
268 
269  // Deallocate the ErrorMap contents if we allocated them earlier
270  if (!error_norm_is_identity)
271  for(unsigned int v = 0; v < n_vars; v++)
272  {
273  delete primal_errors_per_cell[std::make_pair(&_system, v)];
274  delete dual_errors_per_cell[std::make_pair(&_system, v)];
275  delete total_dual_errors_per_cell[std::make_pair(&_system, v)];
276  }
277 
278  STOP_LOG("estimate_error()", "AdjointResidualErrorEstimator");
279 }
280 
281 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo