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                 :            :   auto a = edge.first();
      67                 :            :   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                 :            :   const tk::real small = std::numeric_limits< tk::real >::epsilon();
      98                 :            : 
      99                 :            :   const auto& x = coord[0];
     100                 :            :   const auto& y = coord[1];
     101                 :            :   const auto& z = coord[2];
     102                 :            :   auto a = edge.first();
     103                 :            :   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                 :            :   auto dua = tk::dot( ga, h );
     114                 :            :   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                 :            : }