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 :
19 :
20 : // Local Includes
21 : #include "libmesh/libmesh_common.h"
22 : #include "libmesh/exact_error_estimator.h"
23 : #include "libmesh/dof_map.h"
24 : #include "libmesh/equation_systems.h"
25 : #include "libmesh/error_vector.h"
26 : #include "libmesh/fe_base.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/elem.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/mesh_function.h"
31 : #include "libmesh/numeric_vector.h"
32 : #include "libmesh/quadrature.h"
33 : #include "libmesh/system.h"
34 : #include "libmesh/tensor_tools.h"
35 : #include "libmesh/enum_error_estimator_type.h"
36 : #include "libmesh/enum_norm_type.h"
37 :
38 : // C++ includes
39 : #include <algorithm> // for std::fill
40 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
41 : #include <cmath> // for sqrt
42 : #include <memory>
43 :
44 : namespace libMesh
45 : {
46 :
47 : //-----------------------------------------------------------------
48 : // ErrorEstimator implementations
49 0 : ExactErrorEstimator::ExactErrorEstimator() :
50 : ErrorEstimator(),
51 0 : _exact_value(nullptr),
52 0 : _exact_deriv(nullptr),
53 0 : _exact_hessian(nullptr),
54 0 : _equation_systems_fine(nullptr),
55 0 : _extra_order(1)
56 : {
57 0 : error_norm = H1;
58 0 : }
59 :
60 :
61 0 : ErrorEstimatorType ExactErrorEstimator::type() const
62 : {
63 0 : return EXACT;
64 : }
65 :
66 :
67 0 : void ExactErrorEstimator::attach_exact_value (ValueFunctionPointer fptr)
68 : {
69 0 : libmesh_assert(fptr);
70 0 : _exact_value = fptr;
71 :
72 : // We're not using a fine grid solution
73 0 : _equation_systems_fine = nullptr;
74 :
75 : // We're not using user-provided functors
76 0 : this->clear_functors();
77 0 : }
78 :
79 :
80 0 : void ExactErrorEstimator::attach_exact_values (std::vector<FunctionBase<Number> *> f)
81 : {
82 : // Automatically delete any previous _exact_values entries, then add a new
83 : // entry for each system.
84 0 : _exact_values.clear();
85 :
86 0 : for (auto ptr : f)
87 0 : _exact_values.emplace_back(ptr ? ptr->clone() : nullptr);
88 0 : }
89 :
90 :
91 0 : void ExactErrorEstimator::attach_exact_value (unsigned int sys_num,
92 : FunctionBase<Number> * f)
93 : {
94 0 : if (_exact_values.size() <= sys_num)
95 0 : _exact_values.resize(sys_num+1);
96 :
97 0 : if (f)
98 0 : _exact_values[sys_num] = f->clone();
99 0 : }
100 :
101 :
102 0 : void ExactErrorEstimator::attach_exact_deriv (GradientFunctionPointer gptr)
103 : {
104 0 : libmesh_assert(gptr);
105 0 : _exact_deriv = gptr;
106 :
107 : // We're not using a fine grid solution
108 0 : _equation_systems_fine = nullptr;
109 :
110 : // We're not using user-provided functors
111 0 : this->clear_functors();
112 0 : }
113 :
114 :
115 0 : void ExactErrorEstimator::attach_exact_derivs (std::vector<FunctionBase<Gradient> *> g)
116 : {
117 : // Automatically delete any previous _exact_derivs entries, then add a new
118 : // entry for each system.
119 0 : _exact_derivs.clear();
120 :
121 0 : for (auto ptr : g)
122 0 : _exact_derivs.emplace_back(ptr ? ptr->clone() : nullptr);
123 0 : }
124 :
125 :
126 0 : void ExactErrorEstimator::attach_exact_deriv (unsigned int sys_num,
127 : FunctionBase<Gradient> * g)
128 : {
129 0 : if (_exact_derivs.size() <= sys_num)
130 0 : _exact_derivs.resize(sys_num+1);
131 :
132 0 : if (g)
133 0 : _exact_derivs[sys_num] = g->clone();
134 0 : }
135 :
136 :
137 :
138 :
139 0 : void ExactErrorEstimator::attach_exact_hessian (HessianFunctionPointer hptr)
140 : {
141 0 : libmesh_assert(hptr);
142 0 : _exact_hessian = hptr;
143 :
144 : // We're not using a fine grid solution
145 0 : _equation_systems_fine = nullptr;
146 :
147 : // We're not using user-provided functors
148 0 : this->clear_functors();
149 0 : }
150 :
151 :
152 0 : void ExactErrorEstimator::attach_exact_hessians (std::vector<FunctionBase<Tensor> *> h)
153 : {
154 : // Automatically delete any previous _exact_hessians entries, then add a new
155 : // entry for each system.
156 0 : _exact_hessians.clear();
157 :
158 0 : for (auto ptr : h)
159 0 : _exact_hessians.emplace_back(ptr ? ptr->clone() : nullptr);
160 0 : }
161 :
162 :
163 0 : void ExactErrorEstimator::attach_exact_hessian (unsigned int sys_num,
164 : FunctionBase<Tensor> * h)
165 : {
166 0 : if (_exact_hessians.size() <= sys_num)
167 0 : _exact_hessians.resize(sys_num+1);
168 :
169 0 : if (h)
170 0 : _exact_hessians[sys_num] = h->clone();
171 0 : }
172 :
173 :
174 0 : void ExactErrorEstimator::attach_reference_solution (EquationSystems * es_fine)
175 : {
176 0 : libmesh_assert(es_fine);
177 0 : _equation_systems_fine = es_fine;
178 :
179 : // If we're using a fine grid solution, we're not using exact value
180 : // function pointers or functors.
181 0 : _exact_value = nullptr;
182 0 : _exact_deriv = nullptr;
183 0 : _exact_hessian = nullptr;
184 :
185 0 : this->clear_functors();
186 0 : }
187 :
188 0 : void ExactErrorEstimator::estimate_error (const System & system,
189 : ErrorVector & error_per_cell,
190 : const NumericVector<Number> * solution_vector,
191 : bool estimate_parent_error)
192 : {
193 : // Ignore the fact that this variable is unused when !LIBMESH_ENABLE_AMR
194 0 : libmesh_ignore(estimate_parent_error);
195 :
196 : // The current mesh
197 0 : const MeshBase & mesh = system.get_mesh();
198 :
199 : // The dimensionality of the mesh
200 0 : const unsigned int dim = mesh.mesh_dimension();
201 :
202 : // The number of variables in the system
203 0 : const unsigned int n_vars = system.n_vars();
204 :
205 : // The DofMap for this system
206 0 : const DofMap & dof_map = system.get_dof_map();
207 :
208 : // Resize the error_per_cell vector to be
209 : // the number of elements, initialize it to 0.
210 0 : error_per_cell.resize (mesh.max_elem_id());
211 0 : std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
212 :
213 : // Prepare current_local_solution to localize a non-standard
214 : // solution vector if necessary
215 0 : if (solution_vector && solution_vector != system.solution.get())
216 : {
217 0 : NumericVector<Number> * newsol =
218 : const_cast<NumericVector<Number> *>(solution_vector);
219 0 : System & sys = const_cast<System &>(system);
220 0 : newsol->swap(*sys.solution);
221 0 : sys.update();
222 : }
223 :
224 : // Loop over all the variables in the system
225 0 : for (unsigned int var=0; var<n_vars; var++)
226 : {
227 : // Possibly skip this variable
228 0 : if (error_norm.weight(var) == 0.0) continue;
229 :
230 : // The (string) name of this variable
231 0 : const std::string & var_name = system.variable_name(var);
232 :
233 : // The type of finite element to use for this variable
234 0 : const FEType & fe_type = dof_map.variable_type (var);
235 :
236 0 : std::unique_ptr<FEBase> fe (FEBase::build (dim, fe_type));
237 :
238 : // Build an appropriate Gaussian quadrature rule
239 : std::unique_ptr<QBase> qrule =
240 : fe_type.default_quadrature_rule (dim,
241 0 : _extra_order);
242 :
243 0 : fe->attach_quadrature_rule (qrule.get());
244 :
245 : // Prepare a global solution and a MeshFunction of the fine system if we need one
246 0 : std::unique_ptr<MeshFunction> fine_values;
247 0 : std::unique_ptr<NumericVector<Number>> fine_soln = NumericVector<Number>::build(system.comm());
248 0 : if (_equation_systems_fine)
249 : {
250 0 : const System & fine_system = _equation_systems_fine->get_system(system.name());
251 :
252 0 : std::vector<Number> global_soln;
253 : // FIXME - we're assuming that the fine system solution gets
254 : // used even when a different vector is used for the coarse
255 : // system
256 0 : fine_system.update_global_solution(global_soln);
257 0 : fine_soln->init
258 0 : (cast_int<numeric_index_type>(global_soln.size()), true,
259 0 : SERIAL);
260 0 : (*fine_soln) = global_soln;
261 :
262 : fine_values = std::make_unique<MeshFunction>
263 0 : (*_equation_systems_fine,
264 0 : *fine_soln,
265 : fine_system.get_dof_map(),
266 0 : fine_system.variable_number(var_name));
267 0 : fine_values->init();
268 : } else {
269 : // Initialize functors if we're using them
270 0 : for (auto & ev : _exact_values)
271 0 : if (ev)
272 0 : ev->init();
273 :
274 0 : for (auto & ed : _exact_derivs)
275 0 : if (ed)
276 0 : ed->init();
277 :
278 0 : for (auto & eh : _exact_hessians)
279 0 : if (eh)
280 0 : eh->init();
281 : }
282 :
283 : // Request the data we'll need to compute with
284 0 : fe->get_JxW();
285 0 : fe->get_phi();
286 0 : fe->get_dphi();
287 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
288 0 : fe->get_d2phi();
289 : #endif
290 0 : fe->get_xyz();
291 :
292 : #ifdef LIBMESH_ENABLE_AMR
293 : // If we compute on parent elements, we'll want to do so only
294 : // once on each, so we need to keep track of which we've done.
295 0 : std::vector<bool> computed_var_on_parent;
296 :
297 0 : if (estimate_parent_error)
298 0 : computed_var_on_parent.resize(error_per_cell.size(), false);
299 : #endif
300 :
301 : // TODO: this ought to be threaded (and using subordinate
302 : // MeshFunction objects in each thread rather than a single
303 : // master)
304 :
305 : // Iterate over all the active elements in the mesh
306 : // that live on this processor.
307 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
308 : {
309 0 : const dof_id_type e_id = elem->id();
310 :
311 : #ifdef LIBMESH_ENABLE_AMR
312 : // See if the parent of element e has been examined yet;
313 : // if not, we may want to compute the estimator on it
314 0 : const Elem * parent = elem->parent();
315 :
316 : // We only can compute and only need to compute on
317 : // parents with all active children
318 0 : bool compute_on_parent = true;
319 0 : if (!parent || !estimate_parent_error)
320 0 : compute_on_parent = false;
321 : else
322 0 : for (auto & child : parent->child_ref_range())
323 0 : if (!child.active())
324 0 : compute_on_parent = false;
325 :
326 0 : if (compute_on_parent &&
327 0 : !computed_var_on_parent[parent->id()])
328 : {
329 0 : computed_var_on_parent[parent->id()] = true;
330 :
331 : // Compute a projection onto the parent
332 0 : DenseVector<Number> Uparent;
333 0 : FEBase::coarsened_dof_values(*(system.current_local_solution),
334 : dof_map, parent, Uparent,
335 : var, false);
336 :
337 0 : error_per_cell[parent->id()] +=
338 0 : static_cast<ErrorVectorReal>
339 0 : (find_squared_element_error(system, var_name,
340 : parent, Uparent,
341 : fe.get(),
342 : fine_values.get()));
343 : }
344 : #endif
345 :
346 : // Get the local to global degree of freedom maps
347 0 : std::vector<dof_id_type> dof_indices;
348 0 : dof_map.dof_indices (elem, dof_indices, var);
349 : const unsigned int n_dofs =
350 0 : cast_int<unsigned int>(dof_indices.size());
351 0 : DenseVector<Number> Uelem(n_dofs);
352 0 : for (unsigned int i=0; i != n_dofs; ++i)
353 0 : Uelem(i) = system.current_solution(dof_indices[i]);
354 :
355 0 : error_per_cell[e_id] +=
356 0 : static_cast<ErrorVectorReal>
357 0 : (find_squared_element_error(system, var_name, elem,
358 : Uelem, fe.get(),
359 : fine_values.get()));
360 :
361 0 : } // End loop over active local elements
362 0 : } // End loop over variables
363 :
364 :
365 :
366 : // Each processor has now computed the error contributions
367 : // for its local elements. We need to sum the vector
368 : // and then take the square-root of each component. Note
369 : // that we only need to sum if we are running on multiple
370 : // processors, and we only need to take the square-root
371 : // if the value is nonzero. There will in general be many
372 : // zeros for the inactive elements.
373 :
374 : // First sum the vector of estimated error values
375 0 : this->reduce_error(error_per_cell, system.comm());
376 :
377 : // Compute the square-root of each component.
378 : {
379 0 : LOG_SCOPE("std::sqrt()", "ExactErrorEstimator");
380 0 : for (auto & val : error_per_cell)
381 0 : if (val != 0.)
382 : {
383 0 : libmesh_assert_greater (val, 0.);
384 0 : val = std::sqrt(val);
385 : }
386 : }
387 :
388 : // If we used a non-standard solution before, now is the time to fix
389 : // the current_local_solution
390 0 : if (solution_vector && solution_vector != system.solution.get())
391 : {
392 0 : NumericVector<Number> * newsol =
393 : const_cast<NumericVector<Number> *>(solution_vector);
394 0 : System & sys = const_cast<System &>(system);
395 0 : newsol->swap(*sys.solution);
396 0 : sys.update();
397 : }
398 0 : }
399 :
400 :
401 :
402 0 : Real ExactErrorEstimator::find_squared_element_error(const System & system,
403 : const std::string & var_name,
404 : const Elem * elem,
405 : const DenseVector<Number> & Uelem,
406 : FEBase * fe,
407 : MeshFunction * fine_values) const
408 : {
409 : // The (string) name of this system
410 0 : const std::string & sys_name = system.name();
411 0 : const unsigned int sys_num = system.number();
412 :
413 0 : const unsigned int var = system.variable_number(var_name);
414 : const unsigned int var_component =
415 0 : system.variable_scalar_number(var, 0);
416 :
417 0 : const Parameters & parameters = system.get_equation_systems().parameters;
418 :
419 : // reinitialize the element-specific data
420 : // for the current element
421 0 : fe->reinit (elem);
422 :
423 : // Get the data we need to compute with
424 0 : const std::vector<Real> & JxW = fe->get_JxW();
425 0 : const std::vector<std::vector<Real>> & phi_values = fe->get_phi();
426 0 : const std::vector<std::vector<RealGradient>> & dphi_values = fe->get_dphi();
427 0 : const std::vector<Point> & q_point = fe->get_xyz();
428 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
429 0 : const std::vector<std::vector<RealTensor>> & d2phi_values = fe->get_d2phi();
430 : #endif
431 :
432 : // The number of shape functions
433 : const unsigned int n_sf =
434 0 : cast_int<unsigned int>(Uelem.size());
435 :
436 : // The number of quadrature points
437 : const unsigned int n_qp =
438 0 : cast_int<unsigned int>(JxW.size());
439 :
440 0 : Real error_val = 0;
441 :
442 : // Begin the loop over the Quadrature points.
443 : //
444 0 : for (unsigned int qp=0; qp<n_qp; qp++)
445 : {
446 : // Real u_h = 0.;
447 : // RealGradient grad_u_h;
448 :
449 0 : Number u_h = 0.;
450 :
451 0 : Gradient grad_u_h;
452 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
453 0 : Tensor grad2_u_h;
454 : #endif
455 :
456 : // Compute solution values at the current
457 : // quadrature point. This requires a sum
458 : // over all the shape functions evaluated
459 : // at the quadrature point.
460 0 : for (unsigned int i=0; i<n_sf; i++)
461 : {
462 : // Values from current solution.
463 0 : u_h += phi_values[i][qp]*Uelem(i);
464 0 : grad_u_h += dphi_values[i][qp]*Uelem(i);
465 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
466 0 : grad2_u_h += d2phi_values[i][qp]*Uelem(i);
467 : #endif
468 : }
469 :
470 : // Compute the value of the error at this quadrature point
471 0 : if (error_norm.type(var) == L2 ||
472 0 : error_norm.type(var) == H1 ||
473 0 : error_norm.type(var) == H2)
474 : {
475 0 : Number val_error = u_h;
476 0 : if (_exact_value)
477 0 : val_error -= _exact_value(q_point[qp],parameters,sys_name,var_name);
478 0 : else if (_exact_values.size() > sys_num && _exact_values[sys_num])
479 0 : val_error -= _exact_values[sys_num]->
480 0 : component(var_component, q_point[qp], system.time);
481 0 : else if (_equation_systems_fine)
482 0 : val_error -= (*fine_values)(q_point[qp]);
483 :
484 : // Add the squares of the error to each contribution
485 0 : error_val += JxW[qp]*TensorTools::norm_sq(val_error);
486 : }
487 :
488 : // Compute the value of the error in the gradient at this
489 : // quadrature point
490 0 : if (error_norm.type(var) == H1 ||
491 0 : error_norm.type(var) == H1_SEMINORM ||
492 0 : error_norm.type(var) == H2)
493 : {
494 0 : Gradient grad_error = grad_u_h;
495 0 : if (_exact_deriv)
496 0 : grad_error -= _exact_deriv(q_point[qp],parameters,sys_name,var_name);
497 0 : else if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num])
498 0 : grad_error -= _exact_derivs[sys_num]->
499 0 : component(var_component, q_point[qp], system.time);
500 0 : else if (_equation_systems_fine)
501 0 : grad_error -= fine_values->gradient(q_point[qp]);
502 :
503 0 : error_val += JxW[qp]*grad_error.norm_sq();
504 : }
505 :
506 :
507 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
508 : // Compute the value of the error in the hessian at this
509 : // quadrature point
510 0 : if ((error_norm.type(var) == H2_SEMINORM ||
511 0 : error_norm.type(var) == H2))
512 : {
513 0 : Tensor grad2_error = grad2_u_h;
514 0 : if (_exact_hessian)
515 0 : grad2_error -= _exact_hessian(q_point[qp],parameters,sys_name,var_name);
516 0 : else if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num])
517 0 : grad2_error -= _exact_hessians[sys_num]->
518 0 : component(var_component, q_point[qp], system.time);
519 0 : else if (_equation_systems_fine)
520 0 : grad2_error -= fine_values->hessian(q_point[qp]);
521 :
522 0 : error_val += JxW[qp]*grad2_error.norm_sq();
523 : }
524 : #endif
525 :
526 : } // end qp loop
527 :
528 0 : libmesh_assert_greater_equal (error_val, 0.);
529 :
530 0 : return error_val;
531 : }
532 :
533 :
534 :
535 0 : void ExactErrorEstimator::clear_functors()
536 : {
537 0 : _exact_values.clear();
538 0 : _exact_derivs.clear();
539 0 : _exact_hessians.clear();
540 0 : }
541 :
542 :
543 :
544 : } // namespace libMesh
|
