Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/AMR/Error.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Class for computing error estimates for mesh refinement
9 : : \details Class for computing error estimates for mesh refinement.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <cmath>
14 : : #include <limits>
15 : :
16 : : #include "Exception.hpp"
17 : : #include "Error.hpp"
18 : : #include "Vector.hpp"
19 : : #include "Gradients.hpp"
20 : : #include "Inciter/Options/AMRError.hpp"
21 : :
22 : : using AMR::Error;
23 : :
24 : : tk::real
25 : 890258 : Error::scalar( const tk::Fields& u,
26 : : const edge_t& edge,
27 : : ncomp_t c,
28 : : const std::array< std::vector< tk::real >, 3 >& coord,
29 : : const std::vector< std::size_t >& inpoel,
30 : : const std::pair< std::vector< std::size_t >,
31 : : std::vector< std::size_t > >& esup,
32 : : inciter::ctr::AMRErrorType err ) const
33 : : // *****************************************************************************
34 : : // Estimate error for scalar quantity
35 : : //! \param[in] u Solution vector
36 : : //! \param[in] edge Edge defined by its two end-point IDs
37 : : //! \param[in] c Scalar component to compute error of
38 : : //! \param[in] coord Mesh node coordinates
39 : : //! \param[in] inpoel Mesh element connectivity
40 : : //! \param[in] esup Linked lists storing elements surrounding points, see
41 : : //! tk::genEsup()
42 : : //! \param[in] err AMR Error indicator type
43 : : //! \return Error indicator: a real number between [0...1] inclusive
44 : : // *****************************************************************************
45 : : {
46 [ + + ]: 890258 : if (err == inciter::ctr::AMRErrorType::JUMP)
47 : 789382 : return error_jump( u, edge, c );
48 [ + - ]: 100876 : else if (err == inciter::ctr::AMRErrorType::HESSIAN)
49 : 100876 : return error_hessian( u, edge, c, coord, inpoel, esup );
50 : : else
51 [ - - ][ - - ]: 0 : Throw( "No such AMR error indicator type" );
[ - - ]
52 : : }
53 : :
54 : : tk::real
55 : 789382 : Error::error_jump( const tk::Fields& u,
56 : : const edge_t& edge,
57 : : ncomp_t c ) const
58 : : // *****************************************************************************
59 : : // Estimate error for scalar quantity on edge based on jump in solution
60 : : //! \param[in] u Solution vector
61 : : //! \param[in] edge Edge defined by its two end-point IDs
62 : : //! \param[in] c Scalar component to compute error of
63 : : //! \return Error indicator: a real number between [0...1] inclusive
64 : : // *****************************************************************************
65 : : {
66 : 789382 : auto a = edge.first();
67 : 789382 : auto b = edge.second();
68 : :
69 : : // Normalization factor with a noise filter
70 : 789382 : auto norm = std::abs(u(a,c)) + std::abs(u(b,c)) + 1e-3;
71 : :
72 : : // Normalized error
73 : 789382 : return std::abs( u(a,c) - u(b,c) ) / norm;
74 : : }
75 : :
76 : : tk::real
77 : 100876 : Error::error_hessian( const tk::Fields& u,
78 : : const edge_t& edge,
79 : : ncomp_t c,
80 : : const std::array< std::vector< tk::real >, 3 >& coord,
81 : : const std::vector< std::size_t >& inpoel,
82 : : const std::pair< std::vector< std::size_t >,
83 : : std::vector< std::size_t > >& esup )
84 : : const
85 : : // *****************************************************************************
86 : : // Estimate error for scalar quantity on edge based on Hessian of solution
87 : : //! \param[in] u Solution vector
88 : : //! \param[in] edge Edge defined by its two end-point IDs
89 : : //! \param[in] c Scalar component to compute error of
90 : : //! \param[in] coord Mesh node coordinates
91 : : //! \param[in] inpoel Mesh element connectivity
92 : : //! \param[in] esup Linked lists storing elements surrounding points, see
93 : : //! tk::genEsup()
94 : : //! \return Error indicator: a real number between [0...1] inclusive
95 : : // *****************************************************************************
96 : : {
97 : 100876 : const tk::real small = std::numeric_limits< tk::real >::epsilon();
98 : :
99 : 100876 : const auto& x = coord[0];
100 : 100876 : const auto& y = coord[1];
101 : 100876 : const auto& z = coord[2];
102 : 100876 : auto a = edge.first();
103 : 100876 : auto b = edge.second();
104 : :
105 : : // Compute edge vector
106 : 100876 : std::array< tk::real, 3 > h {{ x[a]-x[b], y[a]-y[b], z[a]-z[b] }};
107 : :
108 : : // Compute gradients at edge-end points
109 [ + - ]: 100876 : auto ga = nodegrad( a, coord, inpoel, esup, u, c );
110 [ + - ]: 100876 : auto gb = nodegrad( b, coord, inpoel, esup, u, c );
111 : :
112 : : // Compute dot products of gradients and edge vectors
113 : 100876 : auto dua = tk::dot( ga, h );
114 : 100876 : auto dub = tk::dot( gb, h );
115 : :
116 : : // If the normalization factor is zero, return zero error
117 : 100876 : auto norm = std::abs(dua) + std::abs(dub);
118 [ + + ]: 100876 : if (norm < small) return 0.0;
119 : :
120 : 50753 : return std::abs(dub-dua) / norm;
121 : : }
|
