Xyst test code coverage report
Current view: top level - Physics - Lohner.cpp (source / functions) Hit Total Coverage
Commit: d790e211db0f2cf155e72db869cf1a5a372c10f5 Lines: 546 548 99.6 %
Date: 2025-02-17 15:57:47 Functions: 14 14 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 378 718 52.6 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Physics/Lohner.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.,
       7                 :            :              2022-2025 J. Bakosi
       8                 :            :              All rights reserved. See the LICENSE file for details.
       9                 :            :   \brief     LohCG: Artificial compressibility solver for incompressible flow
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : 
      13                 :            : #include "Vector.hpp"
      14                 :            : #include "Around.hpp"
      15                 :            : #include "DerivedData.hpp"
      16                 :            : #include "EOS.hpp"
      17                 :            : #include "Lohner.hpp"
      18                 :            : #include "Problems.hpp"
      19                 :            : #include "InciterConfig.hpp"
      20                 :            : 
      21                 :            : namespace inciter {
      22                 :            : 
      23                 :            : extern ctr::Config g_cfg;
      24                 :            : 
      25                 :            : } // ::inciter
      26                 :            : 
      27                 :            : namespace lohner {
      28                 :            : 
      29                 :            : static const tk::real muscl_eps = 1.0e-9;
      30                 :            : static const tk::real muscl_const = 1.0/3.0;
      31                 :            : 
      32                 :            : using inciter::g_cfg;
      33                 :            : 
      34                 :            : void
      35                 :        770 : div( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
      36                 :            :      const std::array< std::vector< tk::real >, 3 >& dsupint,
      37                 :            :      const std::array< std::vector< tk::real >, 3 >& coord,
      38                 :            :      const std::vector< std::size_t >& triinpoel,
      39                 :            :      const tk::Fields& U,
      40                 :            :      std::vector< tk::real >& D,
      41                 :            :      std::size_t pos )
      42                 :            : // *****************************************************************************
      43                 :            : //  Compute divergence of a vector in all points
      44                 :            : //! \param[in] dsupedge Domain superedges
      45                 :            : //! \param[in] dsupint Domain superedge integrals
      46                 :            : //! \param[in] coord Mesh node coordinates
      47                 :            : //! \param[in] triinpoel Boundary face connectivity
      48                 :            : //! \param[in] U Vector whose divergence to compute
      49                 :            : //! \param[in,out] D Vector added to
      50                 :            : //! \param[in] pos Position at which the three vector components are in U
      51                 :            : // *****************************************************************************
      52                 :            : {
      53 [ -  + ][ -  - ]:        770 :   Assert( U.nunk() == D.size(), "Size mismatch" );
         [ -  - ][ -  - ]
      54                 :            : 
      55                 :            :   // Lambda to compute the velocity divergence contribution for edge p-q
      56                 :     197600 :   auto div = [&]( const tk::real d[], std::size_t p, std::size_t q ){
      57                 :     197600 :     return d[0] * (U(p,pos+0) + U(q,pos+0)) +
      58                 :     197600 :            d[1] * (U(p,pos+1) + U(q,pos+1)) +
      59                 :     197600 :            d[2] * (U(p,pos+2) + U(q,pos+2));
      60                 :        770 :   };
      61                 :            : 
      62                 :            :   #if defined(__clang__)
      63                 :            :     #pragma clang diagnostic push
      64                 :            :     #pragma clang diagnostic ignored "-Wvla"
      65                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
      66                 :            :   #elif defined(STRICT_GNUC)
      67                 :            :     #pragma GCC diagnostic push
      68                 :            :     #pragma GCC diagnostic ignored "-Wvla"
      69                 :            :   #endif
      70                 :            : 
      71                 :            :   // domain integral
      72                 :            : 
      73                 :            :   // domain edge contributions: tetrahedron superedges
      74         [ +  + ]:      17772 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
      75                 :      17002 :     const auto N = dsupedge[0].data() + e*4;
      76                 :      17002 :     const auto d = dsupint[0].data();
      77                 :            :     // edge fluxes
      78                 :            :     tk::real f[] = {
      79         [ +  - ]:      17002 :       div( d+(e*6+0)*4, N[0], N[1] ),
      80                 :      34004 :       div( d+(e*6+1)*4, N[1], N[2] ),
      81                 :      34004 :       div( d+(e*6+2)*4, N[2], N[0] ),
      82                 :      34004 :       div( d+(e*6+3)*4, N[0], N[3] ),
      83                 :      34004 :       div( d+(e*6+4)*4, N[1], N[3] ),
      84 [ +  - ][ +  - ]:      17002 :       div( d+(e*6+5)*4, N[2], N[3] ) };
         [ +  - ][ +  - ]
                 [ +  - ]
      85                 :            :     // edge flux contributions
      86                 :      17002 :     D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
      87                 :      17002 :     D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
      88                 :      17002 :     D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
      89                 :      17002 :     D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
      90                 :            :   }
      91                 :            : 
      92                 :            :   // domain edge contributions: triangle superedges
      93         [ +  + ]:      17038 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
      94                 :      16268 :     const auto N = dsupedge[1].data() + e*3;
      95                 :      16268 :     const auto d = dsupint[1].data();
      96                 :            :     // edge fluxes
      97                 :            :     tk::real f[] = {
      98         [ +  - ]:      16268 :       div( d+(e*3+0)*4, N[0], N[1] ),
      99                 :      32536 :       div( d+(e*3+1)*4, N[1], N[2] ),
     100 [ +  - ][ +  - ]:      16268 :       div( d+(e*3+2)*4, N[2], N[0] ) };
     101                 :            :     // edge flux contributions
     102                 :      16268 :     D[N[0]] = D[N[0]] - f[0] + f[2];
     103                 :      16268 :     D[N[1]] = D[N[1]] + f[0] - f[1];
     104                 :      16268 :     D[N[2]] = D[N[2]] + f[1] - f[2];
     105                 :            :   }
     106                 :            : 
     107                 :            :   // domain edge contributions: edges
     108         [ +  + ]:      47554 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     109                 :      46784 :     const auto N = dsupedge[2].data() + e*2;
     110                 :      46784 :     const auto d = dsupint[2].data();
     111                 :            :     // edge flux
     112         [ +  - ]:      46784 :     tk::real f = div( d+e*4, N[0], N[1] );
     113                 :            :     // edge flux contributions
     114                 :      46784 :     D[N[0]] -= f;
     115                 :      46784 :     D[N[1]] += f;
     116                 :            :   }
     117                 :            : 
     118                 :            :   // boundary integral
     119                 :            : 
     120                 :        770 :   const auto& x = coord[0];
     121                 :        770 :   const auto& y = coord[1];
     122                 :        770 :   const auto& z = coord[2];
     123                 :            : 
     124         [ +  + ]:      67074 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
     125                 :      66304 :     const auto N = triinpoel.data() + e*3;
     126                 :            :     tk::real n[3];
     127                 :      66304 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
     128                 :      66304 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
     129                 :            :                   n[0], n[1], n[2] );
     130         [ +  - ]:      66304 :     auto uxA = U(N[0],pos+0);
     131         [ +  - ]:      66304 :     auto uyA = U(N[0],pos+1);
     132         [ +  - ]:      66304 :     auto uzA = U(N[0],pos+2);
     133         [ +  - ]:      66304 :     auto uxB = U(N[1],pos+0);
     134         [ +  - ]:      66304 :     auto uyB = U(N[1],pos+1);
     135         [ +  - ]:      66304 :     auto uzB = U(N[1],pos+2);
     136         [ +  - ]:      66304 :     auto uxC = U(N[2],pos+0);
     137         [ +  - ]:      66304 :     auto uyC = U(N[2],pos+1);
     138         [ +  - ]:      66304 :     auto uzC = U(N[2],pos+2);
     139                 :      66304 :     auto ux = (6.0*uxA + uxB + uxC)/8.0;
     140                 :      66304 :     auto uy = (6.0*uyA + uyB + uyC)/8.0;
     141                 :      66304 :     auto uz = (6.0*uzA + uzB + uzC)/8.0;
     142                 :      66304 :     D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
     143                 :      66304 :     ux = (uxA + 6.0*uxB + uxC)/8.0;
     144                 :      66304 :     uy = (uyA + 6.0*uyB + uyC)/8.0;
     145                 :      66304 :     uz = (uzA + 6.0*uzB + uzC)/8.0;
     146                 :      66304 :     D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
     147                 :      66304 :     ux = (uxA + uxB + 6.0*uxC)/8.0;
     148                 :      66304 :     uy = (uyA + uyB + 6.0*uyC)/8.0;
     149                 :      66304 :     uz = (uzA + uzB + 6.0*uzC)/8.0;
     150                 :      66304 :     D[N[2]] += ux*n[0] + uy*n[1] + uz*n[2];
     151                 :            :   }
     152                 :            : 
     153                 :            :   #if defined(__clang__)
     154                 :            :     #pragma clang diagnostic pop
     155                 :            :   #elif defined(STRICT_GNUC)
     156                 :            :     #pragma GCC diagnostic pop
     157                 :            :   #endif
     158                 :        770 : }
     159                 :            : 
     160                 :            : void
     161                 :        385 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     162                 :            :       const std::array< std::vector< tk::real >, 3 >& dsupint,
     163                 :            :       const std::array< std::vector< tk::real >, 3 >& coord,
     164                 :            :       const std::vector< std::size_t >& triinpoel,
     165                 :            :       const std::vector< tk::real >& U,
     166                 :            :       tk::Fields& G )
     167                 :            : // *****************************************************************************
     168                 :            : //  Compute gradients of scalar in all points
     169                 :            : //! \param[in] dsupedge Domain superedges
     170                 :            : //! \param[in] dsupint Domain superedge integrals
     171                 :            : //! \param[in] coord Mesh node coordinates
     172                 :            : //! \param[in] triinpoel Boundary face connectivity
     173                 :            : //! \param[in] U Scalar whose gradient to compute
     174                 :            : //! \param[in,out] G Nodal gradient of scalar in all points
     175                 :            : // *****************************************************************************
     176                 :            : {
     177         [ -  + ]:        385 :   if (G.nprop() == 0) return;
     178                 :            : 
     179 [ -  + ][ -  - ]:        385 :   Assert( G.nunk() == U.size(), "Size mismatch" );
         [ -  - ][ -  - ]
     180 [ -  + ][ -  - ]:        385 :   Assert( G.nprop() > 2, "Size mismatch" );
         [ -  - ][ -  - ]
     181 [ -  + ][ -  - ]:        385 :   Assert( G.nprop() % 3 == 0, "Size mismatch" );
         [ -  - ][ -  - ]
     182                 :            : 
     183                 :            :   #if defined(__clang__)
     184                 :            :     #pragma clang diagnostic push
     185                 :            :     #pragma clang diagnostic ignored "-Wvla"
     186                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     187                 :            :   #elif defined(STRICT_GNUC)
     188                 :            :     #pragma GCC diagnostic push
     189                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     190                 :            :   #endif
     191                 :            : 
     192                 :            :   // domain integral
     193                 :            : 
     194                 :            :   // domain edge contributions: tetrahedron superedges
     195         [ +  + ]:       8886 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
     196                 :       8501 :     const auto N = dsupedge[0].data() + e*4;
     197                 :       8501 :     const auto d = dsupint[0].data();
     198                 :       8501 :     tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };
     199         [ +  + ]:      34004 :     for (std::size_t j=0; j<3; ++j) {
     200                 :            :       tk::real f[] = {
     201                 :      25503 :         d[(e*6+0)*4+j] * (u[1] + u[0]),
     202                 :      25503 :         d[(e*6+1)*4+j] * (u[2] + u[1]),
     203                 :      25503 :         d[(e*6+2)*4+j] * (u[0] + u[2]),
     204                 :      25503 :         d[(e*6+3)*4+j] * (u[3] + u[0]),
     205                 :      25503 :         d[(e*6+4)*4+j] * (u[3] + u[1]),
     206                 :      25503 :         d[(e*6+5)*4+j] * (u[3] + u[2]) };
     207 [ +  - ][ +  - ]:      25503 :       G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
     208 [ +  - ][ +  - ]:      25503 :       G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
     209 [ +  - ][ +  - ]:      25503 :       G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
     210 [ +  - ][ +  - ]:      25503 :       G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
     211                 :            :     }
     212                 :            :   }
     213                 :            : 
     214                 :            :   // domain edge contributions: triangle superedges
     215         [ +  + ]:       8519 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
     216                 :       8134 :     const auto N = dsupedge[1].data() + e*3;
     217                 :       8134 :     const auto d = dsupint[1].data();
     218                 :       8134 :     tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };
     219         [ +  + ]:      32536 :     for (std::size_t j=0; j<3; ++j) {
     220                 :            :       tk::real f[] = {
     221                 :      24402 :         d[(e*3+0)*4+j] * (u[1] + u[0]),
     222                 :      24402 :         d[(e*3+1)*4+j] * (u[2] + u[1]),
     223                 :      24402 :         d[(e*3+2)*4+j] * (u[0] + u[2]) };
     224 [ +  - ][ +  - ]:      24402 :       G(N[0],j) = G(N[0],j) - f[0] + f[2];
     225 [ +  - ][ +  - ]:      24402 :       G(N[1],j) = G(N[1],j) + f[0] - f[1];
     226 [ +  - ][ +  - ]:      24402 :       G(N[2],j) = G(N[2],j) + f[1] - f[2];
     227                 :            :     }
     228                 :            :   }
     229                 :            : 
     230                 :            :   // domain edge contributions: edges
     231         [ +  + ]:      23777 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     232                 :      23392 :     const auto N = dsupedge[2].data() + e*2;
     233                 :      23392 :     const auto d = dsupint[2].data() + e*4;
     234                 :      23392 :     tk::real u[] = { U[N[0]], U[N[1]] };
     235         [ +  + ]:      93568 :     for (std::size_t j=0; j<3; ++j) {
     236                 :      70176 :       tk::real f = d[j] * (u[1] + u[0]);
     237         [ +  - ]:      70176 :       G(N[0],j) -= f;
     238         [ +  - ]:      70176 :       G(N[1],j) += f;
     239                 :            :     }
     240                 :            :   }
     241                 :            : 
     242                 :            :   // boundary integral
     243                 :            : 
     244                 :        385 :   const auto& x = coord[0];
     245                 :        385 :   const auto& y = coord[1];
     246                 :        385 :   const auto& z = coord[2];
     247                 :            : 
     248         [ +  + ]:      33537 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
     249                 :      33152 :     const auto N = triinpoel.data() + e*3;
     250                 :            :     tk::real n[3];
     251                 :      33152 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
     252                 :      33152 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
     253                 :            :                   n[0], n[1], n[2] );
     254                 :      33152 :     auto uA = U[N[0]];
     255                 :      33152 :     auto uB = U[N[1]];
     256                 :      33152 :     auto uC = U[N[2]];
     257                 :      33152 :     auto f = (6.0*uA + uB + uC)/8.0;
     258         [ +  - ]:      33152 :     G(N[0],0) += f * n[0];
     259         [ +  - ]:      33152 :     G(N[0],1) += f * n[1];
     260         [ +  - ]:      33152 :     G(N[0],2) += f * n[2];
     261                 :      33152 :     f = (uA + 6.0*uB + uC)/8.0;
     262         [ +  - ]:      33152 :     G(N[1],0) += f * n[0];
     263         [ +  - ]:      33152 :     G(N[1],1) += f * n[1];
     264         [ +  - ]:      33152 :     G(N[1],2) += f * n[2];
     265                 :      33152 :     f = (uA + uB + 6.0*uC)/8.0;
     266         [ +  - ]:      33152 :     G(N[2],0) += f * n[0];
     267         [ +  - ]:      33152 :     G(N[2],1) += f * n[1];
     268         [ +  - ]:      33152 :     G(N[2],2) += f * n[2];
     269                 :            :   }
     270                 :            : 
     271                 :            :   #if defined(__clang__)
     272                 :            :     #pragma clang diagnostic pop
     273                 :            :   #elif defined(STRICT_GNUC)
     274                 :            :     #pragma GCC diagnostic pop
     275                 :            :   #endif
     276                 :            : }
     277                 :            : 
     278                 :            : void
     279                 :        385 : vgrad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     280                 :            :        const std::array< std::vector< tk::real >, 3 >& dsupint,
     281                 :            :        const std::array< std::vector< tk::real >, 3 >& coord,
     282                 :            :        const std::vector< std::size_t >& triinpoel,
     283                 :            :        const tk::Fields& U,
     284                 :            :        tk::Fields& G )
     285                 :            : // *****************************************************************************
     286                 :            : //  Compute velocity gradients in all points
     287                 :            : //! \param[in] dsupedge Domain superedges
     288                 :            : //! \param[in] dsupint Domain superedge integrals
     289                 :            : //! \param[in] coord Mesh node coordinates
     290                 :            : //! \param[in] triinpoel Boundary face connectivity
     291                 :            : //! \param[in] U Velocity whose gradient to compute
     292                 :            : //! \param[in,out] G Nodal velocity gradients (9 components) in all points
     293                 :            : // *****************************************************************************
     294                 :            : {
     295 [ -  + ][ -  - ]:        385 :   Assert( G.nunk() == U.nunk(), "Size mismatch" );
         [ -  - ][ -  - ]
     296 [ -  + ][ -  - ]:        385 :   Assert( G.nprop() == 9, "Size mismatch" );
         [ -  - ][ -  - ]
     297                 :            : 
     298                 :            :   #if defined(__clang__)
     299                 :            :     #pragma clang diagnostic push
     300                 :            :     #pragma clang diagnostic ignored "-Wvla"
     301                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     302                 :            :   #elif defined(STRICT_GNUC)
     303                 :            :     #pragma GCC diagnostic push
     304                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     305                 :            :   #endif
     306                 :            : 
     307                 :            :   // domain integral
     308                 :            : 
     309                 :            :   // domain edge contributions: tetrahedron superedges
     310         [ +  + ]:       8886 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
     311                 :       8501 :     const auto N = dsupedge[0].data() + e*4;
     312                 :       8501 :     const auto d = dsupint[0].data();
     313         [ +  + ]:      34004 :     for (std::size_t i=0; i<3; ++i) {
     314 [ +  - ][ +  - ]:      25503 :       tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1), U(N[3],i+1) };
         [ +  - ][ +  - ]
     315                 :      25503 :       auto i3 = i*3;
     316         [ +  + ]:     102012 :       for (std::size_t j=0; j<3; ++j) {
     317                 :      76509 :         tk::real f[] = { d[(e*6+0)*4+j] * (u[1] + u[0]),
     318                 :      76509 :                          d[(e*6+1)*4+j] * (u[2] + u[1]),
     319                 :      76509 :                          d[(e*6+2)*4+j] * (u[0] + u[2]),
     320                 :      76509 :                          d[(e*6+3)*4+j] * (u[3] + u[0]),
     321                 :      76509 :                          d[(e*6+4)*4+j] * (u[3] + u[1]),
     322                 :      76509 :                          d[(e*6+5)*4+j] * (u[3] + u[2]) };
     323 [ +  - ][ +  - ]:      76509 :         G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
     324 [ +  - ][ +  - ]:      76509 :         G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
     325 [ +  - ][ +  - ]:      76509 :         G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
     326 [ +  - ][ +  - ]:      76509 :         G(N[3],i3+j) = G(N[3],i3+j) + f[3] + f[4] + f[5];
     327                 :            :       }
     328                 :            :     }
     329                 :            :   }
     330                 :            : 
     331                 :            :   // domain edge contributions: triangle superedges
     332         [ +  + ]:       8519 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
     333                 :       8134 :     const auto N = dsupedge[1].data() + e*3;
     334                 :       8134 :     const auto d = dsupint[1].data();
     335         [ +  + ]:      32536 :     for (std::size_t i=0; i<3; ++i) {
     336 [ +  - ][ +  - ]:      24402 :       tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
                 [ +  - ]
     337                 :      24402 :       auto i3 = i*3;
     338         [ +  + ]:      97608 :       for (std::size_t j=0; j<3; ++j) {
     339                 :      73206 :         tk::real f[] = { d[(e*3+0)*4+j] * (u[1] + u[0]),
     340                 :      73206 :                          d[(e*3+1)*4+j] * (u[2] + u[1]),
     341                 :      73206 :                          d[(e*3+2)*4+j] * (u[0] + u[2]) };
     342 [ +  - ][ +  - ]:      73206 :         G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
     343 [ +  - ][ +  - ]:      73206 :         G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
     344 [ +  - ][ +  - ]:      73206 :         G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
     345                 :            :       }
     346                 :            :     }
     347                 :            :   }
     348                 :            : 
     349                 :            :   // domain edge contributions: edges
     350         [ +  + ]:      23777 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     351                 :      23392 :     const auto N = dsupedge[2].data() + e*2;
     352                 :      23392 :     const auto d = dsupint[2].data() + e*4;
     353         [ +  + ]:      93568 :     for (std::size_t i=0; i<3; ++i) {
     354 [ +  - ][ +  - ]:      70176 :       tk::real u[] = { U(N[0],i+1), U(N[1],i+1) };
     355                 :      70176 :       auto i3 = i*3;
     356         [ +  + ]:     280704 :       for (std::size_t j=0; j<3; ++j) {
     357                 :     210528 :         tk::real f = d[j] * (u[1] + u[0]);
     358         [ +  - ]:     210528 :         G(N[0],i3+j) -= f;
     359         [ +  - ]:     210528 :         G(N[1],i3+j) += f;
     360                 :            :       }
     361                 :            :     }
     362                 :            :   }
     363                 :            : 
     364                 :            :   // boundary integral
     365                 :            : 
     366                 :        385 :   const auto& x = coord[0];
     367                 :        385 :   const auto& y = coord[1];
     368                 :        385 :   const auto& z = coord[2];
     369                 :            : 
     370         [ +  + ]:      33537 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
     371                 :      33152 :     const auto N = triinpoel.data() + e*3;
     372                 :            :     tk::real n[3];
     373                 :      33152 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
     374                 :      33152 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
     375                 :            :                   n[0], n[1], n[2] );
     376         [ +  + ]:     132608 :     for (std::size_t i=0; i<3; ++i) {
     377 [ +  - ][ +  - ]:      99456 :       tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
                 [ +  - ]
     378                 :      99456 :       auto i3 = i*3;
     379                 :      99456 :       auto f = (6.0*u[0] + u[1] + u[2])/8.0;
     380         [ +  - ]:      99456 :       G(N[0],i3+0) += f * n[0];
     381         [ +  - ]:      99456 :       G(N[0],i3+1) += f * n[1];
     382         [ +  - ]:      99456 :       G(N[0],i3+2) += f * n[2];
     383                 :      99456 :       f = (u[0] + 6.0*u[1] + u[2])/8.0;
     384         [ +  - ]:      99456 :       G(N[1],i3+0) += f * n[0];
     385         [ +  - ]:      99456 :       G(N[1],i3+1) += f * n[1];
     386         [ +  - ]:      99456 :       G(N[1],i3+2) += f * n[2];
     387                 :      99456 :       f = (u[0] + u[1] + 6.0*u[2])/8.0;
     388         [ +  - ]:      99456 :       G(N[2],i3+0) += f * n[0];
     389         [ +  - ]:      99456 :       G(N[2],i3+1) += f * n[1];
     390         [ +  - ]:      99456 :       G(N[2],i3+2) += f * n[2];
     391                 :            :     }
     392                 :            :   }
     393                 :            : 
     394                 :            :   #if defined(__clang__)
     395                 :            :     #pragma clang diagnostic pop
     396                 :            :   #elif defined(STRICT_GNUC)
     397                 :            :     #pragma GCC diagnostic pop
     398                 :            :   #endif
     399                 :        385 : }
     400                 :            : 
     401                 :            : static tk::real
     402                 :     889200 : flux( const tk::Fields& U,
     403                 :            :       const tk::Fields& G,
     404                 :            :       std::size_t i,
     405                 :            :       std::size_t j,
     406                 :            :       std::size_t p,
     407                 :            :       std::size_t q )
     408                 :            : // *****************************************************************************
     409                 :            : //! Compute momentum flux over edge of points p-q
     410                 :            : //! \param[in] U Velocity vector
     411                 :            : //! \param[in] G Velocity gradients
     412                 :            : //! \param[in] i Tensor component, 1st index
     413                 :            : //! \param[in] j Tensor component, 2nd index
     414                 :            : //! \param[in] p Left node index of edge
     415                 :            : //! \param[in] q Right node index of edge
     416                 :            : //! \return Momentum flux contribution for edge p-q
     417                 :            : // *****************************************************************************
     418                 :            : {
     419                 :     889200 :   auto inv = U(p,i+1)*U(p,j+1) + U(q,i+1)*U(q,j+1);
     420                 :            : 
     421                 :     889200 :   auto eps = std::numeric_limits< tk::real >::epsilon();
     422                 :     889200 :   auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
     423         [ +  + ]:     889200 :   if (mu < eps) return -inv;
     424                 :            : 
     425                 :     655731 :   auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
     426         [ +  + ]:     655731 :   if (i == j) {
     427                 :     218577 :     vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) + G(q,0) + G(q,4) + G(q,8) );
     428                 :            :   }
     429                 :     655731 :   return mu*vis - inv;
     430                 :            : }
     431                 :            : 
     432                 :            : static tk::real
     433                 :     895104 : flux( const tk::Fields& U,
     434                 :            :       const tk::Fields& G,
     435                 :            :       std::size_t i,
     436                 :            :       std::size_t j,
     437                 :            :       std::size_t p )
     438                 :            : // *****************************************************************************
     439                 :            : //! Compute momentum flux in point p
     440                 :            : //! \param[in] U Velocity vector
     441                 :            : //! \param[in] G Velocity gradients
     442                 :            : //! \param[in] i Tensor component, 1st index
     443                 :            : //! \param[in] j Tensor component, 2nd index
     444                 :            : //! \param[in] p Node index of point
     445                 :            : //! \return Momentum flux contribution for point p
     446                 :            : // *****************************************************************************
     447                 :            : {
     448                 :     895104 :   auto inv = U(p,i+1)*U(p,j+1);
     449                 :            : 
     450                 :     895104 :   auto eps = std::numeric_limits< tk::real >::epsilon();
     451                 :     895104 :   auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
     452         [ +  + ]:     895104 :   if (mu < eps) return -inv;
     453                 :            : 
     454                 :     636120 :   auto vis = G(p,i*3+j) + G(p,j*3+i);
     455         [ +  + ]:     636120 :   if (i == j) {
     456                 :     212040 :     vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
     457                 :            :   }
     458                 :     636120 :   return mu*vis - inv;
     459                 :            : }
     460                 :            : 
     461                 :            : void
     462                 :        385 : flux( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     463                 :            :       const std::array< std::vector< tk::real >, 3 >& dsupint,
     464                 :            :       const std::array< std::vector< tk::real >, 3 >& coord,
     465                 :            :       const std::vector< std::size_t >& triinpoel,
     466                 :            :       const tk::Fields& U,
     467                 :            :       const tk::Fields& G,
     468                 :            :       tk::Fields& F )
     469                 :            : // *****************************************************************************
     470                 :            : //  Compute momentum flux in all points
     471                 :            : //! \param[in] dsupedge Domain superedges
     472                 :            : //! \param[in] dsupint Domain superedge integrals
     473                 :            : //! \param[in] coord Mesh node coordinates
     474                 :            : //! \param[in] triinpoel Boundary face connectivity
     475                 :            : //! \param[in] U Velocity field
     476                 :            : //! \param[in] G Velocity gradients, dui/dxj, 9 components
     477                 :            : //! \param[in,out] F Momentum flux, Fi = d[ sij - uiuj ] / dxj, where
     478                 :            : //!    s_ij = mu[ dui/dxj + duj/dxi - 2/3 duk/dxk delta_ij ]
     479                 :            : // *****************************************************************************
     480                 :            : {
     481 [ -  + ][ -  - ]:        385 :   Assert( F.nunk() == U.nunk(), "Size mismatch" );
         [ -  - ][ -  - ]
     482 [ -  + ][ -  - ]:        385 :   Assert( F.nprop() == 3, "Size mismatch" );
         [ -  - ][ -  - ]
     483 [ -  + ][ -  - ]:        385 :   Assert( G.nunk() == U.nunk(), "Size mismatch" );
         [ -  - ][ -  - ]
     484 [ -  + ][ -  - ]:        385 :   Assert( G.nprop() == 9, "Size mismatch" );
         [ -  - ][ -  - ]
     485                 :            : 
     486                 :            :   #if defined(__clang__)
     487                 :            :     #pragma clang diagnostic push
     488                 :            :     #pragma clang diagnostic ignored "-Wvla"
     489                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     490                 :            :   #elif defined(STRICT_GNUC)
     491                 :            :     #pragma GCC diagnostic push
     492                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     493                 :            :   #endif
     494                 :            : 
     495                 :            :   // domain integral
     496                 :            : 
     497                 :            :   // domain edge contributions: tetrahedron superedges
     498         [ +  + ]:       8886 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
     499                 :       8501 :     const auto N = dsupedge[0].data() + e*4;
     500                 :       8501 :     const auto d = dsupint[0].data();
     501         [ +  + ]:      34004 :     for (std::size_t i=0; i<3; ++i) {
     502         [ +  + ]:     102012 :       for (std::size_t j=0; j<3; ++j) {
     503         [ +  - ]:      76509 :         tk::real f[] = { d[(e*6+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
     504                 :     153018 :                          d[(e*6+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
     505                 :     153018 :                          d[(e*6+2)*4+j] * flux(U,G,i,j,N[0],N[2]),
     506                 :     153018 :                          d[(e*6+3)*4+j] * flux(U,G,i,j,N[3],N[0]),
     507                 :     153018 :                          d[(e*6+4)*4+j] * flux(U,G,i,j,N[3],N[1]),
     508 [ +  - ][ +  - ]:      76509 :                          d[(e*6+5)*4+j] * flux(U,G,i,j,N[3],N[2]) };
         [ +  - ][ +  - ]
                 [ +  - ]
     509 [ +  - ][ +  - ]:      76509 :         F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
     510 [ +  - ][ +  - ]:      76509 :         F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
     511 [ +  - ][ +  - ]:      76509 :         F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
     512 [ +  - ][ +  - ]:      76509 :         F(N[3],i) = F(N[3],i) + f[3] + f[4] + f[5];
     513                 :            :       }
     514                 :            :     }
     515                 :            :   }
     516                 :            : 
     517                 :            :   // domain edge contributions: triangle superedges
     518         [ +  + ]:       8519 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
     519                 :       8134 :     const auto N = dsupedge[1].data() + e*3;
     520                 :       8134 :     const auto d = dsupint[1].data();
     521         [ +  + ]:      32536 :     for (std::size_t i=0; i<3; ++i) {
     522         [ +  + ]:      97608 :       for (std::size_t j=0; j<3; ++j) {
     523         [ +  - ]:      73206 :         tk::real f[] = { d[(e*3+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
     524                 :     146412 :                          d[(e*3+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
     525 [ +  - ][ +  - ]:      73206 :                          d[(e*3+2)*4+j] * flux(U,G,i,j,N[0],N[2]), };
     526 [ +  - ][ +  - ]:      73206 :         F(N[0],i) = F(N[0],i) - f[0] + f[2];
     527 [ +  - ][ +  - ]:      73206 :         F(N[1],i) = F(N[1],i) + f[0] - f[1];
     528 [ +  - ][ +  - ]:      73206 :         F(N[2],i) = F(N[2],i) + f[1] - f[2];
     529                 :            :       }
     530                 :            :     }
     531                 :            :   }
     532                 :            : 
     533                 :            :   // domain edge contributions: edges
     534         [ +  + ]:      23777 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     535                 :      23392 :     const auto N = dsupedge[2].data() + e*2;
     536                 :      23392 :     const auto d = dsupint[2].data() + e*4;
     537         [ +  + ]:      93568 :     for (std::size_t i=0; i<3; ++i) {
     538         [ +  + ]:     280704 :       for (std::size_t j=0; j<3; ++j) {
     539                 :     210528 :         tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
     540                 :     210528 :         F(N[0],i) -= f;
     541                 :     210528 :         F(N[1],i) += f;
     542                 :            :       }
     543                 :            :     }
     544                 :            :   }
     545                 :            : 
     546                 :            :   // boundary integral
     547                 :            : 
     548                 :        385 :   const auto& x = coord[0];
     549                 :        385 :   const auto& y = coord[1];
     550                 :        385 :   const auto& z = coord[2];
     551                 :            : 
     552         [ +  + ]:      33537 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
     553                 :      33152 :     const auto N = triinpoel.data() + e*3;
     554                 :            :     tk::real n[3];
     555                 :      33152 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
     556                 :      33152 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
     557                 :            :                   n[0], n[1], n[2] );
     558         [ +  + ]:     132608 :     for (std::size_t i=0; i<3; ++i) {
     559         [ +  - ]:      99456 :       auto fxA = flux(U,G,i,0,N[0]);
     560         [ +  - ]:      99456 :       auto fyA = flux(U,G,i,1,N[0]);
     561         [ +  - ]:      99456 :       auto fzA = flux(U,G,i,2,N[0]);
     562         [ +  - ]:      99456 :       auto fxB = flux(U,G,i,0,N[1]);
     563         [ +  - ]:      99456 :       auto fyB = flux(U,G,i,1,N[1]);
     564         [ +  - ]:      99456 :       auto fzB = flux(U,G,i,2,N[1]);
     565         [ +  - ]:      99456 :       auto fxC = flux(U,G,i,0,N[2]);
     566         [ +  - ]:      99456 :       auto fyC = flux(U,G,i,1,N[2]);
     567         [ +  - ]:      99456 :       auto fzC = flux(U,G,i,2,N[2]);
     568                 :      99456 :       auto fx = (6.0*fxA + fxB + fxC)/8.0;
     569                 :      99456 :       auto fy = (6.0*fyA + fyB + fyC)/8.0;
     570                 :      99456 :       auto fz = (6.0*fzA + fzB + fzC)/8.0;
     571         [ +  - ]:      99456 :       F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
     572                 :      99456 :       fx = (fxA + 6.0*fxB + fxC)/8.0;
     573                 :      99456 :       fy = (fyA + 6.0*fyB + fyC)/8.0;
     574                 :      99456 :       fz = (fzA + 6.0*fzB + fzC)/8.0;
     575         [ +  - ]:      99456 :       F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
     576                 :      99456 :       fx = (fxA + fxB + 6.0*fxC)/8.0;
     577                 :      99456 :       fy = (fyA + fyB + 6.0*fyC)/8.0;
     578                 :      99456 :       fz = (fzA + fzB + 6.0*fzC)/8.0;
     579         [ +  - ]:      99456 :       F(N[2],i) += fx*n[0] + fy*n[1] + fz*n[2];
     580                 :            :     }
     581                 :            :   }
     582                 :            : 
     583                 :            :   #if defined(__clang__)
     584                 :            :     #pragma clang diagnostic pop
     585                 :            :   #elif defined(STRICT_GNUC)
     586                 :            :     #pragma GCC diagnostic pop
     587                 :            :   #endif
     588                 :        385 : }
     589                 :            : 
     590                 :            : void
     591                 :       1040 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     592                 :            :       const std::array< std::vector< tk::real >, 3 >& dsupint,
     593                 :            :       const std::array< std::vector< tk::real >, 3 >& coord,
     594                 :            :       const std::vector< std::size_t >& triinpoel,
     595                 :            :       const tk::Fields& U,
     596                 :            :       tk::Fields& G )
     597                 :            : // *****************************************************************************
     598                 :            : //  Compute gradients of scalar in all points
     599                 :            : //! \param[in] dsupedge Domain superedges
     600                 :            : //! \param[in] dsupint Domain superedge integrals
     601                 :            : //! \param[in] coord Mesh node coordinates
     602                 :            : //! \param[in] triinpoel Boundary face connectivity
     603                 :            : //! \param[in] U Unknown vector whose gradient to compute
     604                 :            : //! \param[in,out] G Nodal gradients in all points
     605                 :            : // *****************************************************************************
     606                 :            : {
     607         [ -  + ]:       1040 :   if (G.nprop() == 0) return;
     608                 :            : 
     609 [ -  + ][ -  - ]:       1040 :   Assert( G.nunk() == U.nunk(), "Size mismatch" );
         [ -  - ][ -  - ]
     610 [ -  + ][ -  - ]:       1040 :   Assert( G.nprop() > 2, "Size mismatch" );
         [ -  - ][ -  - ]
     611 [ -  + ][ -  - ]:       1040 :   Assert( G.nprop() % 3 == 0, "Size mismatch" );
         [ -  - ][ -  - ]
     612 [ -  + ][ -  - ]:       1040 :   Assert( G.nprop() == U.nprop()*3, "Size mismatch" );
         [ -  - ][ -  - ]
     613                 :            : 
     614                 :       1040 :   const auto ncomp = U.nprop();
     615                 :            : 
     616                 :            :   #if defined(__clang__)
     617                 :            :     #pragma clang diagnostic push
     618                 :            :     #pragma clang diagnostic ignored "-Wvla"
     619                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     620                 :            :   #elif defined(STRICT_GNUC)
     621                 :            :     #pragma GCC diagnostic push
     622                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     623                 :            :   #endif
     624                 :            : 
     625                 :            :   // domain integral
     626                 :            : 
     627                 :            :   // domain edge contributions: tetrahedron superedges
     628         [ +  + ]:     175120 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
     629                 :     174080 :     const auto N = dsupedge[0].data() + e*4;
     630                 :     174080 :     const auto d = dsupint[0].data();
     631         [ +  + ]:     910880 :     for (std::size_t c=0; c<ncomp; ++c) {
     632 [ +  - ][ +  - ]:     736800 :       tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c), U(N[3],c) };
         [ +  - ][ +  - ]
     633                 :     736800 :       auto g = c*3;
     634         [ +  + ]:    2947200 :       for (std::size_t j=0; j<3; ++j) {
     635                 :            :         tk::real f[] = {
     636                 :    2210400 :           d[(e*6+0)*4+j] * (u[1] + u[0]),
     637                 :    2210400 :           d[(e*6+1)*4+j] * (u[2] + u[1]),
     638                 :    2210400 :           d[(e*6+2)*4+j] * (u[0] + u[2]),
     639                 :    2210400 :           d[(e*6+3)*4+j] * (u[3] + u[0]),
     640                 :    2210400 :           d[(e*6+4)*4+j] * (u[3] + u[1]),
     641                 :    2210400 :           d[(e*6+5)*4+j] * (u[3] + u[2]) };
     642 [ +  - ][ +  - ]:    2210400 :         G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2] - f[3];
     643 [ +  - ][ +  - ]:    2210400 :         G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1] - f[4];
     644 [ +  - ][ +  - ]:    2210400 :         G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2] - f[5];
     645 [ +  - ][ +  - ]:    2210400 :         G(N[3],g+j) = G(N[3],g+j) + f[3] + f[4] + f[5];
     646                 :            :       }
     647                 :            :     }
     648                 :            :   }
     649                 :            : 
     650                 :            :   // domain edge contributions: triangle superedges
     651         [ +  + ]:     149440 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
     652                 :     148400 :     const auto N = dsupedge[1].data() + e*3;
     653                 :     148400 :     const auto d = dsupint[1].data();
     654         [ +  + ]:     779600 :     for (std::size_t c=0; c<ncomp; ++c) {
     655 [ +  - ][ +  - ]:     631200 :       tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
                 [ +  - ]
     656                 :     631200 :       auto g = c*3;
     657         [ +  + ]:    2524800 :       for (std::size_t j=0; j<3; ++j) {
     658                 :            :         tk::real f[] = {
     659                 :    1893600 :           d[(e*3+0)*4+j] * (u[1] + u[0]),
     660                 :    1893600 :           d[(e*3+1)*4+j] * (u[2] + u[1]),
     661                 :    1893600 :           d[(e*3+2)*4+j] * (u[0] + u[2]) };
     662 [ +  - ][ +  - ]:    1893600 :         G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2];
     663 [ +  - ][ +  - ]:    1893600 :         G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1];
     664 [ +  - ][ +  - ]:    1893600 :         G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2];
     665                 :            :       }
     666                 :            :     }
     667                 :            :   }
     668                 :            : 
     669                 :            :   // domain edge contributions: edges
     670         [ +  + ]:     455440 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     671                 :     454400 :     const auto N = dsupedge[2].data() + e*2;
     672                 :     454400 :     const auto d = dsupint[2].data() + e*4;
     673         [ +  + ]:    2399680 :     for (std::size_t c=0; c<ncomp; ++c) {
     674 [ +  - ][ +  - ]:    1945280 :       tk::real u[] = { U(N[0],c), U(N[1],c) };
     675                 :    1945280 :       auto g = c*3;
     676         [ +  + ]:    7781120 :       for (std::size_t j=0; j<3; ++j) {
     677                 :    5835840 :         tk::real f = d[j] * (u[1] + u[0]);
     678         [ +  - ]:    5835840 :         G(N[0],g+j) -= f;
     679         [ +  - ]:    5835840 :         G(N[1],g+j) += f;
     680                 :            :       }
     681                 :            :     }
     682                 :            :   }
     683                 :            : 
     684                 :            :   // boundary integral
     685                 :            : 
     686                 :       1040 :   const auto& x = coord[0];
     687                 :       1040 :   const auto& y = coord[1];
     688                 :       1040 :   const auto& z = coord[2];
     689                 :            : 
     690         [ +  + ]:     585200 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
     691                 :     584160 :     const auto N = triinpoel.data() + e*3;
     692                 :            :     tk::real n[3];
     693                 :     584160 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
     694                 :     584160 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
     695                 :            :                   n[0], n[1], n[2] );
     696         [ +  + ]:    3112640 :     for (std::size_t c=0; c<ncomp; ++c) {
     697                 :    2528480 :       auto g = c*3;
     698         [ +  - ]:    2528480 :       auto uA = U(N[0],c);
     699         [ +  - ]:    2528480 :       auto uB = U(N[1],c);
     700         [ +  - ]:    2528480 :       auto uC = U(N[2],c);
     701                 :    2528480 :       auto f = (6.0*uA + uB + uC)/8.0;
     702         [ +  - ]:    2528480 :       G(N[0],g+0) += f * n[0];
     703         [ +  - ]:    2528480 :       G(N[0],g+1) += f * n[1];
     704         [ +  - ]:    2528480 :       G(N[0],g+2) += f * n[2];
     705                 :    2528480 :       f = (uA + 6.0*uB + uC)/8.0;
     706         [ +  - ]:    2528480 :       G(N[1],g+0) += f * n[0];
     707         [ +  - ]:    2528480 :       G(N[1],g+1) += f * n[1];
     708         [ +  - ]:    2528480 :       G(N[1],g+2) += f * n[2];
     709                 :    2528480 :       f = (uA + uB + 6.0*uC)/8.0;
     710         [ +  - ]:    2528480 :       G(N[2],g+0) += f * n[0];
     711         [ +  - ]:    2528480 :       G(N[2],g+1) += f * n[1];
     712         [ +  - ]:    2528480 :       G(N[2],g+2) += f * n[2];
     713                 :            :     }
     714                 :            :   }
     715                 :            : 
     716                 :            :   #if defined(__clang__)
     717                 :            :     #pragma clang diagnostic pop
     718                 :            :   #elif defined(STRICT_GNUC)
     719                 :            :     #pragma GCC diagnostic pop
     720                 :            :   #endif
     721                 :            : }
     722                 :            : 
     723                 :            : static void
     724                 :    2992380 : adv_damp2( const tk::real supint[],
     725                 :            :            const tk::Fields& U,
     726                 :            :            const tk::Fields&,
     727                 :            :            const std::array< std::vector< tk::real >, 3 >&,
     728                 :            :            std::size_t p,
     729                 :            :            std::size_t q,
     730                 :            :            tk::real f[] )
     731                 :            : // *****************************************************************************
     732                 :            : //! Compute advection fluxes on a single edge using 2nd-order damping
     733                 :            : //! \param[in] supint Edge integral
     734                 :            : //! \param[in] U Velocity and transported scalars at recent time step
     735                 :            : //! \param[in] p Left node index of edge
     736                 :            : //! \param[in] q Right node index of edge
     737                 :            : //! \param[in,out] f Flux computed
     738                 :            : // *****************************************************************************
     739                 :            : {
     740                 :    2992380 :   auto ncomp = U.nprop();
     741                 :            : 
     742                 :    2992380 :   auto nx = supint[0];
     743                 :    2992380 :   auto ny = supint[1];
     744                 :    2992380 :   auto nz = supint[2];
     745                 :            : 
     746                 :            :   // left state
     747                 :    2992380 :   auto pL = U(p,0);
     748                 :    2992380 :   auto uL = U(p,1);
     749                 :    2992380 :   auto vL = U(p,2);
     750                 :    2992380 :   auto wL = U(p,3);
     751                 :    2992380 :   auto vnL = uL*nx + vL*ny + wL*nz;
     752                 :            : 
     753                 :            :   // right state
     754                 :    2992380 :   auto pR = U(q,0);
     755                 :    2992380 :   auto uR = U(q,1);
     756                 :    2992380 :   auto vR = U(q,2);
     757                 :    2992380 :   auto wR = U(q,3);
     758                 :    2992380 :   auto vnR = uR*nx + vR*ny + wR*nz;
     759                 :            : 
     760                 :    2992380 :   auto s = g_cfg.get< tag::soundspeed >();
     761                 :    2992380 :   auto s2 = s*s;
     762                 :            : 
     763                 :            :   // viscosity
     764                 :    2992380 :   auto v = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
     765                 :            : 
     766                 :            :   // stabilization
     767                 :    2992380 :   tk::real aw = 0.0;
     768         [ +  - ]:    2992380 :   if (g_cfg.get< tag::stab >()) {
     769                 :    2992380 :     aw = std::abs( vnL + vnR ) / 2.0;
     770                 :            :   }
     771                 :            : 
     772                 :            :   // artificial viscosity
     773         [ +  + ]:    2992380 :   if (g_cfg.get< tag::stab2 >()) {
     774                 :    1193940 :     auto len = tk::length( nx, ny, nz );
     775                 :    1193940 :     auto sl = std::abs(vnL) + s*len;
     776                 :    1193940 :     auto sr = std::abs(vnR) + s*len;
     777                 :    1193940 :     aw += g_cfg.get< tag::stab2coef >() * std::max(sl,sr);
     778                 :            :   }
     779                 :            : 
     780                 :            :   // flow
     781                 :    2992380 :   auto pf = pL + pR;
     782                 :    2992380 :   f[0] = (vnL + vnR + aw*(pR - pL))*s2;
     783                 :    2992380 :   f[1] = uL*vnL + uR*vnR + pf*nx + (aw-v)*(uR - uL);
     784                 :    2992380 :   f[2] = vL*vnL + vR*vnR + pf*ny + (aw-v)*(vR - vL);
     785                 :    2992380 :   f[3] = wL*vnL + wR*vnR + pf*nz + (aw-v)*(wR - wL);
     786                 :            : 
     787                 :            :   // diffusion
     788                 :    2992380 :   auto d = supint[3] * g_cfg.get< tag::mat_dyn_diffusivity >();
     789                 :            : 
     790                 :            :   // scalar
     791         [ +  + ]:    4186320 :   for (std::size_t c=4; c<ncomp; ++c) {
     792                 :    1193940 :     f[c] = U(p,c)*vnL + U(q,c)*vnR + (aw-d)*(U(q,c) - U(p,c));
     793                 :            :   }
     794                 :    2992380 : }
     795                 :            : 
     796                 :            : static void
     797                 :    1944080 : adv_damp4( const tk::real supint[],
     798                 :            :            const tk::Fields& U,
     799                 :            :            const tk::Fields& G,
     800                 :            :            const std::array< std::vector< tk::real >, 3 >& coord,
     801                 :            :            std::size_t p,
     802                 :            :            std::size_t q,
     803                 :            :            tk::real f[] )
     804                 :            : // *****************************************************************************
     805                 :            : //! Compute advection fluxes on a single edge using 4th-order damping
     806                 :            : //! \param[in] supint Edge integral
     807                 :            : //! \param[in] U Velocity and transported scalars at recent time step
     808                 :            : //! \param[in] G Gradients of velocity and transported scalars
     809                 :            : //! \param[in] coord Mesh node coordinates
     810                 :            : //! \param[in] p Left node index of edge
     811                 :            : //! \param[in] q Right node index of edge
     812                 :            : //! \param[in,out] f Flux computed
     813                 :            : // *****************************************************************************
     814                 :            : {
     815                 :    1944080 :   const auto ncomp = U.nprop();
     816                 :            : 
     817                 :    1944080 :   const auto& x = coord[0];
     818                 :    1944080 :   const auto& y = coord[1];
     819                 :    1944080 :   const auto& z = coord[2];
     820                 :            : 
     821                 :            :   // edge vector
     822                 :    1944080 :   auto dx = x[q] - x[p];
     823                 :    1944080 :   auto dy = y[q] - y[p];
     824                 :    1944080 :   auto dz = z[q] - z[p];
     825                 :            : 
     826                 :            :   #if defined(__clang__)
     827                 :            :     #pragma clang diagnostic push
     828                 :            :     #pragma clang diagnostic ignored "-Wvla"
     829                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     830                 :            :   #elif defined(STRICT_GNUC)
     831                 :            :     #pragma GCC diagnostic push
     832                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     833                 :            :   #endif
     834                 :            : 
     835                 :    1944080 :   tk::real uL[ncomp], uR[ncomp];
     836         [ +  + ]:   10203760 :   for (std::size_t i=0; i<ncomp; ++i) {
     837         [ +  - ]:    8259680 :     uL[i] = U(p,i);
     838         [ +  - ]:    8259680 :     uR[i] = U(q,i);
     839                 :            :   }
     840                 :            : 
     841                 :            :   // MUSCL reconstruction in edge-end points
     842         [ +  + ]:   10203760 :   for (std::size_t c=0; c<ncomp; ++c) {
     843 [ +  - ][ +  - ]:    8259680 :     auto g1 = G(p,c*3+0)*dx + G(p,c*3+1)*dy + G(p,c*3+2)*dz;
                 [ +  - ]
     844 [ +  - ][ +  - ]:    8259680 :     auto g2 = G(q,c*3+0)*dx + G(q,c*3+1)*dy + G(q,c*3+2)*dz;
                 [ +  - ]
     845                 :    8259680 :     auto delta2 = uR[c] - uL[c];
     846                 :    8259680 :     auto delta1 = 2.0 * g1 - delta2;
     847                 :    8259680 :     auto delta3 = 2.0 * g2 - delta2;
     848                 :            : 
     849                 :            :     // van Leer limiter
     850                 :    8259680 :     auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
     851                 :    8259680 :     auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
     852                 :    8259680 :     auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
     853                 :    8259680 :     auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
     854                 :    8259680 :     auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
     855                 :    8259680 :     auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
     856                 :    8259680 :     auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
     857                 :    8259680 :     auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
     858                 :            :     // update unknowns with reconstructed unknowns
     859                 :    8259680 :     uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
     860                 :    8259680 :                    delta2*(1.0+muscl_const)*phi_L_inv);
     861                 :    8259680 :     uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
     862                 :    8259680 :                    delta2*(1.0+muscl_const)*phi_R_inv);
     863                 :            :   }
     864                 :            : 
     865                 :    1944080 :   auto nx = supint[0];
     866                 :    1944080 :   auto ny = supint[1];
     867                 :    1944080 :   auto nz = supint[2];
     868                 :            : 
     869                 :            :   // normal velocities
     870                 :    1944080 :   auto vnL = uL[1]*nx + uL[2]*ny + uL[3]*nz;
     871                 :    1944080 :   auto vnR = uR[1]*nx + uR[2]*ny + uR[3]*nz;
     872                 :            : 
     873                 :    1944080 :   auto s = g_cfg.get< tag::soundspeed >();
     874                 :    1944080 :   auto s2 = s*s;
     875                 :            : 
     876                 :            :   // viscosity
     877                 :    1944080 :   auto v = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
     878                 :            : 
     879                 :            :   // stabilization
     880                 :    1944080 :   tk::real aw = 0.0;
     881         [ +  - ]:    1944080 :   if (g_cfg.get< tag::stab >()) {
     882                 :    1944080 :     aw = std::abs( vnL + vnR ) / 2.0;
     883                 :            :   }
     884                 :            : 
     885                 :            :   // artificial viscosity
     886         [ +  + ]:    1944080 :   if (g_cfg.get< tag::stab2 >()) {
     887                 :    1510400 :     auto len = tk::length( nx, ny, nz );
     888                 :    1510400 :     auto sl = std::abs(vnL) + s*len;
     889                 :    1510400 :     auto sr = std::abs(vnR) + s*len;
     890                 :    1510400 :     aw += g_cfg.get< tag::stab2coef >() * std::max(sl,sr);
     891                 :            :   }
     892                 :            : 
     893                 :            :   // flow
     894                 :    1944080 :   auto pf = uL[0] + uR[0];
     895                 :    1944080 :   f[0] = (vnL + vnR + aw*(uR[0]-uL[0]))*s2;
     896 [ +  - ][ +  - ]:    1944080 :   f[1] = uL[1]*vnL + uR[1]*vnR + pf*nx + aw*(uR[1]-uL[1]) - v*(U(q,1)-U(p,1));
     897 [ +  - ][ +  - ]:    1944080 :   f[2] = uL[2]*vnL + uR[2]*vnR + pf*ny + aw*(uR[2]-uL[2]) - v*(U(q,2)-U(p,2));
     898 [ +  - ][ +  - ]:    1944080 :   f[3] = uL[3]*vnL + uR[3]*vnR + pf*nz + aw*(uR[3]-uL[3]) - v*(U(q,3)-U(p,3));
     899                 :            : 
     900                 :            :   // diffusion
     901                 :    1944080 :   auto d = supint[3] * g_cfg.get< tag::mat_dyn_diffusivity >();
     902                 :            : 
     903                 :            :   // scalar
     904         [ +  + ]:    2427440 :   for (std::size_t c=4; c<ncomp; ++c) {
     905 [ +  - ][ +  - ]:     483360 :     f[c] = uL[c]*vnL + uR[c]*vnR + aw*(uR[c]-uL[c]) - d*(U(q,c)-U(p,c));
     906                 :            :   }
     907                 :            : 
     908                 :            :   #if defined(__clang__)
     909                 :            :     #pragma clang diagnostic pop
     910                 :            :   #elif defined(STRICT_GNUC)
     911                 :            :     #pragma GCC diagnostic pop
     912                 :            :   #endif
     913                 :    3888160 : }
     914                 :            : 
     915                 :            : static void
     916                 :       8600 : adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     917                 :            :      const std::array< std::vector< tk::real >, 3 >& dsupint,
     918                 :            :      const std::array< std::vector< tk::real >, 3 >& coord,
     919                 :            :      const std::vector< std::size_t >& triinpoel,
     920                 :            :      const tk::Fields& U,
     921                 :            :      const tk::Fields& G,
     922                 :            :      // cppcheck-suppress constParameter
     923                 :            :      tk::Fields& R )
     924                 :            : // *****************************************************************************
     925                 :            : //! Add advection to rhs
     926                 :            : //! \param[in] dsupedge Domain superedges
     927                 :            : //! \param[in] dsupint Domain superedge integrals
     928                 :            : //! \param[in] coord Mesh node coordinates
     929                 :            : //! \param[in] triinpoel Boundary face connectivity
     930                 :            : //! \param[in] U Velocity and transported scalars at recent time step
     931                 :            : //! \param[in] G Gradients of velocity and transported scalars
     932                 :            : //! \param[in,out] R Right-hand side vector added to
     933                 :            : // *****************************************************************************
     934                 :            : {
     935                 :       8600 :   const auto ncomp = U.nprop();
     936                 :            : 
     937                 :            :   // configure advection
     938                 :          0 :   auto adv = [](){
     939                 :       8600 :     const auto& flux = g_cfg.get< tag::flux >();
     940         [ +  + ]:       8600 :          if (flux == "damp2") return adv_damp2;
     941         [ +  - ]:       1040 :     else if (flux == "damp4") return adv_damp4;
     942 [ -  - ][ -  - ]:          0 :     else Throw( "Flux not correctly configured" );
                 [ -  - ]
     943         [ +  - ]:       8600 :   }();
     944                 :            : 
     945                 :            :   #if defined(__clang__)
     946                 :            :     #pragma clang diagnostic push
     947                 :            :     #pragma clang diagnostic ignored "-Wvla"
     948                 :            :     #pragma clang diagnostic ignored "-Wvla-extension"
     949                 :            :   #elif defined(STRICT_GNUC)
     950                 :            :     #pragma GCC diagnostic push
     951                 :            :     #pragma GCC diagnostic ignored "-Wvla"
     952                 :            :   #endif
     953                 :            : 
     954                 :            :   // domain integral
     955                 :            : 
     956                 :            :   // domain edge contributions: tetrahedron superedges
     957         [ +  + ]:     440090 :   for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
     958                 :     431490 :     const auto N = dsupedge[0].data() + e*4;
     959                 :     431490 :     tk::real f[6][ncomp];
     960                 :     431490 :     const auto d = dsupint[0].data();
     961         [ +  - ]:     431490 :     adv( d+(e*6+0)*4, U, G, coord, N[0], N[1], f[0] );
     962         [ +  - ]:     431490 :     adv( d+(e*6+1)*4, U, G, coord, N[1], N[2], f[1] );
     963         [ +  - ]:     431490 :     adv( d+(e*6+2)*4, U, G, coord, N[2], N[0], f[2] );
     964         [ +  - ]:     431490 :     adv( d+(e*6+3)*4, U, G, coord, N[0], N[3], f[3] );
     965         [ +  - ]:     431490 :     adv( d+(e*6+4)*4, U, G, coord, N[1], N[3], f[4] );
     966         [ +  - ]:     431490 :     adv( d+(e*6+5)*4, U, G, coord, N[2], N[3], f[5] );
     967         [ +  + ]:    2296690 :     for (std::size_t c=0; c<ncomp; ++c) {
     968 [ +  - ][ +  - ]:    1865200 :       R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
     969 [ +  - ][ +  - ]:    1865200 :       R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
     970 [ +  - ][ +  - ]:    1865200 :       R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
     971 [ +  - ][ +  - ]:    1865200 :       R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
     972                 :            :     }
     973                 :     431490 :   }
     974                 :            : 
     975                 :            :   // domain edge contributions: triangle superedges
     976         [ +  + ]:     412120 :   for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
     977                 :     403520 :     const auto N = dsupedge[1].data() + e*3;
     978                 :     403520 :     tk::real f[3][ncomp];
     979                 :     403520 :     const auto d = dsupint[1].data();
     980         [ +  - ]:     403520 :     adv( d+(e*3+0)*4, U, G, coord, N[0], N[1], f[0] );
     981         [ +  - ]:     403520 :     adv( d+(e*3+1)*4, U, G, coord, N[1], N[2], f[1] );
     982         [ +  - ]:     403520 :     adv( d+(e*3+2)*4, U, G, coord, N[2], N[0], f[2] );
     983         [ +  + ]:    2149820 :     for (std::size_t c=0; c<ncomp; ++c) {
     984 [ +  - ][ +  - ]:    1746300 :       R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
     985 [ +  - ][ +  - ]:    1746300 :       R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
     986 [ +  - ][ +  - ]:    1746300 :       R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
     987                 :            :     }
     988                 :     403520 :   }
     989                 :            : 
     990                 :            :   // domain edge contributions: edges
     991         [ +  + ]:    1145560 :   for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
     992                 :    1136960 :     const auto N = dsupedge[2].data() + e*2;
     993                 :    1136960 :     tk::real f[ncomp];
     994                 :    1136960 :     const auto d = dsupint[2].data();
     995         [ +  - ]:    1136960 :     adv( d+e*4, U, G, coord, N[0], N[1], f );
     996         [ +  + ]:    6130000 :     for (std::size_t c=0; c<ncomp; ++c) {
     997         [ +  - ]:    4993040 :       R(N[0],c) -= f[c];
     998         [ +  - ]:    4993040 :       R(N[1],c) += f[c];
     999                 :            :     }
    1000                 :    1136960 :   }
    1001                 :            : 
    1002                 :            :   // boundary integral
    1003                 :            : 
    1004                 :       8600 :   const auto& x = coord[0];
    1005                 :       8600 :   const auto& y = coord[1];
    1006                 :       8600 :   const auto& z = coord[2];
    1007                 :            : 
    1008                 :       8600 :   auto s = g_cfg.get< tag::soundspeed >();
    1009                 :       8600 :   auto s2 = s * s;
    1010                 :            : 
    1011         [ +  + ]:    1743920 :   for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    1012                 :    3470640 :     const auto N = triinpoel.data() + e*3;
    1013                 :            :     tk::real n[3];
    1014                 :    1735320 :     tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
    1015                 :    1735320 :                   x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
    1016                 :            :                   n[0], n[1], n[2] );
    1017                 :    1735320 :     tk::real f[ncomp][3];
    1018                 :            : 
    1019         [ +  - ]:    1735320 :     auto p = U(N[0],0);
    1020         [ +  - ]:    1735320 :     auto u = U(N[0],1);
    1021         [ +  - ]:    1735320 :     auto v = U(N[0],2);
    1022         [ +  - ]:    1735320 :     auto w = U(N[0],3);
    1023                 :    1735320 :     auto vn = n[0]*u + n[1]*v + n[2]*w;
    1024                 :            :     // flow
    1025                 :    1735320 :     f[0][0] = vn * s2;
    1026                 :    1735320 :     f[1][0] = u*vn + p*n[0];
    1027                 :    1735320 :     f[2][0] = v*vn + p*n[1];
    1028                 :    1735320 :     f[3][0] = w*vn + p*n[2];
    1029                 :            :     // scalar
    1030 [ +  - ][ +  + ]:    2358800 :     for (std::size_t c=4; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
    1031                 :            : 
    1032         [ +  - ]:    1735320 :     p = U(N[1],0);
    1033         [ +  - ]:    1735320 :     u = U(N[1],1);
    1034         [ +  - ]:    1735320 :     v = U(N[1],2);
    1035         [ +  - ]:    1735320 :     w = U(N[1],3);
    1036                 :    1735320 :     vn = n[0]*u + n[1]*v + n[2]*w;
    1037                 :            :     // flow
    1038                 :    1735320 :     f[0][1] = vn * s2;
    1039                 :    1735320 :     f[1][1] = u*vn + p*n[0];
    1040                 :    1735320 :     f[2][1] = v*vn + p*n[1];
    1041                 :    1735320 :     f[3][1] = w*vn + p*n[2];
    1042                 :            :     // scalar
    1043 [ +  - ][ +  + ]:    2358800 :     for (std::size_t c=4; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
    1044                 :            : 
    1045         [ +  - ]:    1735320 :     p = U(N[2],0);
    1046         [ +  - ]:    1735320 :     u = U(N[2],1);
    1047         [ +  - ]:    1735320 :     v = U(N[2],2);
    1048         [ +  - ]:    1735320 :     w = U(N[2],3);
    1049                 :    1735320 :     vn = n[0]*u + n[1]*v + n[2]*w;
    1050                 :            :     // flow
    1051                 :    1735320 :     f[0][2] = vn * s2;
    1052                 :    1735320 :     f[1][2] = u*vn + p*n[0];
    1053                 :    1735320 :     f[2][2] = v*vn + p*n[1];
    1054                 :    1735320 :     f[3][2] = w*vn + p*n[2];
    1055                 :            :     // scalar
    1056 [ +  - ][ +  + ]:    2358800 :     for (std::size_t c=4; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
    1057                 :            : 
    1058         [ +  + ]:    9300080 :     for (std::size_t c=0; c<ncomp; ++c) {
    1059         [ +  - ]:    7564760 :       R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
    1060         [ +  - ]:    7564760 :       R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
    1061         [ +  - ]:    7564760 :       R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
    1062                 :            :     }
    1063                 :    1735320 :   }
    1064                 :            : 
    1065                 :            :   #if defined(__clang__)
    1066                 :            :     #pragma clang diagnostic pop
    1067                 :            :   #elif defined(STRICT_GNUC)
    1068                 :            :     #pragma GCC diagnostic pop
    1069                 :            :   #endif
    1070                 :       8600 : }
    1071                 :            : 
    1072                 :            : static void
    1073                 :       8600 : src( const std::array< std::vector< tk::real >, 3 >& coord,
    1074                 :            :      const std::vector< tk::real >& v,
    1075                 :            :      tk::real t,
    1076                 :            :      tk::Fields& R )
    1077                 :            : // *****************************************************************************
    1078                 :            : //  Compute source integral
    1079                 :            : //! \param[in] coord Mesh node coordinates
    1080                 :            : //! \param[in] v Nodal mesh volumes without contributions from other chares
    1081                 :            : //! \param[in] t Physical time
    1082                 :            : //! \param[in,out] R Right-hand side vector computed
    1083                 :            : // *****************************************************************************
    1084                 :            : {
    1085         [ +  - ]:       8600 :   auto src = problems::SRC();
    1086         [ +  + ]:       8600 :   if (!src) return;
    1087                 :            : 
    1088                 :       3020 :   const auto& x = coord[0];
    1089                 :       3020 :   const auto& y = coord[1];
    1090                 :       3020 :   const auto& z = coord[2];
    1091                 :            : 
    1092         [ +  + ]:     372160 :   for (std::size_t p=0; p<R.nunk(); ++p) {
    1093         [ +  - ]:     369140 :     auto s = src( x[p], y[p], z[p], t );
    1094 [ +  - ][ +  + ]:    2214840 :     for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
    1095                 :     369140 :   }
    1096         [ +  + ]:       8600 : }
    1097                 :            : 
    1098                 :            : void
    1099                 :       8600 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
    1100                 :            :      const std::array< std::vector< tk::real >, 3 >& dsupint,
    1101                 :            :      const std::array< std::vector< tk::real >, 3 >& coord,
    1102                 :            :      const std::vector< std::size_t >& triinpoel,
    1103                 :            :      const std::vector< tk::real >& v,
    1104                 :            :      tk::real t,
    1105                 :            :      const tk::Fields& U,
    1106                 :            :      const tk::Fields& G,
    1107                 :            :      tk::Fields& R )
    1108                 :            : // *****************************************************************************
    1109                 :            : //  Compute right hand side
    1110                 :            : //! \param[in] dsupedge Domain superedges
    1111                 :            : //! \param[in] dsupint Domain superedge integrals
    1112                 :            : //! \param[in] coord Mesh node coordinates
    1113                 :            : //! \param[in] triinpoel Boundary face connectivity
    1114                 :            : //! \param[in] v Nodal mesh volumes without contributions from other chares
    1115                 :            : //! \param[in] t Physical time
    1116                 :            : //! \param[in] U Solution vector of primitive variables at recent time step
    1117                 :            : //! \param[in] G Gradients of velocity and transported scalars
    1118                 :            : //! \param[in,out] R Right-hand side vector computed
    1119                 :            : // *****************************************************************************
    1120                 :            : {
    1121 [ -  + ][ -  - ]:       8600 :   Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
    1122                 :            :           "vector at recent time step incorrect" );
    1123 [ -  + ][ -  - ]:       8600 :   Assert( R.nunk() == coord[0].size(),
         [ -  - ][ -  - ]
    1124                 :            :           "Number of unknowns and/or number of components in right-hand "
    1125                 :            :           "side vector incorrect" );
    1126                 :            : 
    1127                 :       8600 :   R.fill( 0.0 );
    1128                 :       8600 :   adv( dsupedge, dsupint, coord, triinpoel, U, G, R );
    1129                 :       8600 :   src( coord, v, t, R );
    1130                 :       8600 : }
    1131                 :            : 
    1132                 :            : } // lohner::

Generated by: LCOV version 1.16