Cppcheck report - [Xyst]: /home/jbakosi/code/xyst/src/Physics/Chorin.cpp

// *****************************************************************************
/*!
  \file      src/Physics/Chorin.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.,
             2022-2024 J. Bakosi
             All rights reserved. See the LICENSE file for details.
  \brief     ChoCG: Projection-based solver for incompressible flow
*/
// *****************************************************************************

#include "Vector.hpp"
#include "Around.hpp"
#include "DerivedData.hpp"
#include "EOS.hpp"
#include "Chorin.hpp"
#include "Problems.hpp"
#include "InciterConfig.hpp"

namespace inciter {

extern ctr::Config g_cfg;

} // ::inciter

namespace chorin {

static const tk::real muscl_eps = 1.0e-9;
static const tk::real muscl_const = 1.0/3.0;

using inciter::g_cfg;

static tk::real
div( const std::array< std::vector< tk::real >, 3 >& coord,
     const tk::real d[],
     tk::real dt,
     const std::vector< tk::real >& P,
     const tk::Fields& G,
     const tk::Fields& U,
     std::size_t p,
     std::size_t q,
     bool stab )
// *****************************************************************************
//! Compute divergence over edge of points p-q with 4th-order pressure damping
//! \param[in] coord Mesh node coordinates
//! \param[in] d Edge integral
//! \param[in] dt Physical time step size
//! \param[in] P Pressure
//! \param[in] G Pressure gradients
//! \param[in] U Velocity vector
//! \param[in] p Left node index of edge
//! \param[in] q Right node index of edge
//! \param[in] stab True to stabilize
//! \return Divergence contribution for edge p-q
// *****************************************************************************
{
  tk::real div = d[0] * (U(p,0) + U(q,0)) +
                 d[1] * (U(p,1) + U(q,1)) +
                 d[2] * (U(p,2) + U(q,2));

  if (!stab) return div;

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  auto dx = x[p] - x[q];
  auto dy = y[p] - y[q];
  auto dz = z[p] - z[q];
  auto dl = std::sqrt( dx*dx + dy*dy + dz*dz );
  auto p2 = P[q] - P[p];
  auto D = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] );

  auto dpx = G(p,0) + G(q,0);
  auto dpy = G(p,1) + G(q,1);
  auto dpz = G(p,2) + G(q,2);
  auto p4 = 0.5 * (dx*dpx + dy*dpy + dz*dpz);

  div += D*dt/dl*(p2 + p4);

  return div;
}

void
div( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     const std::array< std::vector< tk::real >, 3 >& dsupint,
     const std::array< std::vector< tk::real >, 3 >& coord,
     const std::vector< std::size_t >& triinpoel,
     tk::real dt,
     const std::vector< tk::real >& P,
     const tk::Fields& G,
     const tk::Fields& U,
     std::vector< tk::real >& D,
     bool stab )
// *****************************************************************************
//  Compute divergence of a vector in all points
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] dt Physical time size
//! \param[in] P Pressure
//! \param[in] G Pressure gradients
//! \param[in] U Vector whose divergence to compute
//! \param[in,out] D Nodal divergence of vector in all points
//! \param[in] stab True to stabilize
// *****************************************************************************
{
  Assert( G.nunk() == U.nunk(), "Size mismatch" );
  Assert( G.nprop() > 2, "Size mismatch" );

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // domain integral

  // domain edge contributions: tetrahedron superedges
  for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
    const auto N = dsupedge[0].data() + e*4;
    const auto d = dsupint[0].data();
    // edge fluxes
    tk::real f[] = {<--- Variable 'f' can be declared as const array
      div( coord, d+(e*6+0)*5, dt, P, G, U, N[0], N[1], stab ),
      div( coord, d+(e*6+1)*5, dt, P, G, U, N[1], N[2], stab ),
      div( coord, d+(e*6+2)*5, dt, P, G, U, N[2], N[0], stab ),
      div( coord, d+(e*6+3)*5, dt, P, G, U, N[0], N[3], stab ),
      div( coord, d+(e*6+4)*5, dt, P, G, U, N[1], N[3], stab ),
      div( coord, d+(e*6+5)*5, dt, P, G, U, N[2], N[3], stab ) };
    // edge flux contributions
    D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
    D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
    D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
    D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
  }

  // domain edge contributions: triangle superedges
  for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
    const auto N = dsupedge[1].data() + e*3;
    const auto d = dsupint[1].data();
    // edge fluxes
    tk::real f[] = {<--- Variable 'f' can be declared as const array
      div( coord, d+(e*3+0)*5, dt, P, G, U, N[0], N[1], stab ),
      div( coord, d+(e*3+1)*5, dt, P, G, U, N[1], N[2], stab ),
      div( coord, d+(e*3+2)*5, dt, P, G, U, N[2], N[0], stab ) };
    // edge flux contributions
    D[N[0]] = D[N[0]] - f[0] + f[2];
    D[N[1]] = D[N[1]] + f[0] - f[1];
    D[N[2]] = D[N[2]] + f[1] - f[2];
  }

  // domain edge contributions: edges
  for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
    const auto N = dsupedge[2].data() + e*2;
    const auto d = dsupint[2].data();
    // edge flux
    tk::real f = div( coord, d+e*5, dt, P, G, U, N[0], N[1], stab );
    // edge flux contributions
    D[N[0]] -= f;
    D[N[1]] += f;
  }

  // boundary integral

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    const auto N = triinpoel.data() + e*3;
    tk::real n[3];
    tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
                  x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
                  n[0], n[1], n[2] );
    auto uxA = U(N[0],0);
    auto uyA = U(N[0],1);
    auto uzA = U(N[0],2);
    auto uxB = U(N[1],0);
    auto uyB = U(N[1],1);
    auto uzB = U(N[1],2);
    auto uxC = U(N[2],0);
    auto uyC = U(N[2],1);
    auto uzC = U(N[2],2);
    auto ux = (6.0*uxA + uxB + uxC)/8.0;
    auto uy = (6.0*uyA + uyB + uyC)/8.0;
    auto uz = (6.0*uzA + uzB + uzC)/8.0;
    D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
    ux = (uxA + 6.0*uxB + uxC)/8.0;
    uy = (uyA + 6.0*uyB + uyC)/8.0;
    uz = (uzA + 6.0*uzB + uzC)/8.0;
    D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
    ux = (uxA + uxB + 6.0*uxC)/8.0;
    uy = (uyA + uyB + 6.0*uyC)/8.0;
    uz = (uzA + uzB + 6.0*uzC)/8.0;
    D[N[2]] += ux*n[0] + uy*n[1] + uz*n[2];
  }

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

void
vgrad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
       const std::array< std::vector< tk::real >, 3 >& dsupint,
       const std::array< std::vector< tk::real >, 3 >& coord,
       const std::vector< std::size_t >& triinpoel,
       const tk::Fields& U,
       tk::Fields& G )
// *****************************************************************************
//  Compute velocity gradients in all points
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] U Velocity whose gradient to compute
//! \param[in,out] G Nodal velocity gradients (9 components) in all points
// *****************************************************************************
{
  Assert( G.nunk() == U.nunk(), "Size mismatch" );
  Assert( G.nprop() == 9, "Size mismatch" );

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // domain integral

  // domain edge contributions: tetrahedron superedges
  for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
    const auto N = dsupedge[0].data() + e*4;
    const auto d = dsupint[0].data();
    for (std::size_t i=0; i<3; ++i) {
      tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i), U(N[3],i) };<--- Variable 'u' can be declared as const array
      auto i3 = i*3;
      for (std::size_t j=0; j<3; ++j) {
        tk::real f[] = { d[(e*6+0)*5+j] * (u[1] + u[0]),<--- Variable 'f' can be declared as const array
                         d[(e*6+1)*5+j] * (u[2] + u[1]),
                         d[(e*6+2)*5+j] * (u[0] + u[2]),
                         d[(e*6+3)*5+j] * (u[3] + u[0]),
                         d[(e*6+4)*5+j] * (u[3] + u[1]),
                         d[(e*6+5)*5+j] * (u[3] + u[2]) };
        G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
        G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
        G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
        G(N[3],i3+j) = G(N[3],i3+j) + f[3] + f[4] + f[5];
      }
    }
  }

  // domain edge contributions: triangle superedges
  for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
    const auto N = dsupedge[1].data() + e*3;
    const auto d = dsupint[1].data();
    for (std::size_t i=0; i<3; ++i) {
      tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };<--- Variable 'u' can be declared as const array
      auto i3 = i*3;
      for (std::size_t j=0; j<3; ++j) {
        tk::real f[] = { d[(e*3+0)*5+j] * (u[1] + u[0]),<--- Variable 'f' can be declared as const array
                         d[(e*3+1)*5+j] * (u[2] + u[1]),
                         d[(e*3+2)*5+j] * (u[0] + u[2]) };
        G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
        G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
        G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
      }
    }
  }

  // domain edge contributions: edges
  for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
    const auto N = dsupedge[2].data() + e*2;
    const auto d = dsupint[2].data() + e*5;
    for (std::size_t i=0; i<3; ++i) {
      tk::real u[] = { U(N[0],i), U(N[1],i) };<--- Variable 'u' can be declared as const array
      auto i3 = i*3;
      for (std::size_t j=0; j<3; ++j) {
        tk::real f = d[j] * (u[1] + u[0]);
        G(N[0],i3+j) -= f;
        G(N[1],i3+j) += f;
      }
    }
  }

  // boundary integral

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    const auto N = triinpoel.data() + e*3;
    tk::real n[3];
    tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
                  x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
                  n[0], n[1], n[2] );
    for (std::size_t i=0; i<3; ++i) {
      tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };<--- Variable 'u' can be declared as const array
      auto i3 = i*3;
      auto f = (6.0*u[0] + u[1] + u[2])/8.0;
      G(N[0],i3+0) += f * n[0];
      G(N[0],i3+1) += f * n[1];
      G(N[0],i3+2) += f * n[2];
      f = (u[0] + 6.0*u[1] + u[2])/8.0;
      G(N[1],i3+0) += f * n[0];
      G(N[1],i3+1) += f * n[1];
      G(N[1],i3+2) += f * n[2];
      f = (u[0] + u[1] + 6.0*u[2])/8.0;
      G(N[2],i3+0) += f * n[0];
      G(N[2],i3+1) += f * n[1];
      G(N[2],i3+2) += f * n[2];
    }
  }

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

void
grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
      const std::array< std::vector< tk::real >, 3 >& dsupint,
      const std::array< std::vector< tk::real >, 3 >& coord,
      const std::vector< std::size_t >& triinpoel,
      const std::vector< tk::real >& U,
      tk::Fields& G )
// *****************************************************************************
//  Compute gradients of scalar in all points
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] U Scalar whose gradient to compute
//! \param[in,out] G Nodal gradient of scalar in all points
// *****************************************************************************
{
  Assert( G.nunk() == U.size(), "Size mismatch" );
  Assert( G.nprop() > 2, "Size mismatch" );

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // domain integral

  // domain edge contributions: tetrahedron superedges
  for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
    const auto N = dsupedge[0].data() + e*4;
    const auto d = dsupint[0].data();
    tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };<--- Variable 'u' can be declared as const array
    for (std::size_t j=0; j<3; ++j) {
      tk::real f[] = {<--- Variable 'f' can be declared as const array
        d[(e*6+0)*5+j] * (u[1] + u[0]),
        d[(e*6+1)*5+j] * (u[2] + u[1]),
        d[(e*6+2)*5+j] * (u[0] + u[2]),
        d[(e*6+3)*5+j] * (u[3] + u[0]),
        d[(e*6+4)*5+j] * (u[3] + u[1]),
        d[(e*6+5)*5+j] * (u[3] + u[2]) };
      G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
      G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
      G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
      G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
    }
  }

  // domain edge contributions: triangle superedges
  for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
    const auto N = dsupedge[1].data() + e*3;
    const auto d = dsupint[1].data();
    tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };<--- Variable 'u' can be declared as const array
    for (std::size_t j=0; j<3; ++j) {
      tk::real f[] = {<--- Variable 'f' can be declared as const array
        d[(e*3+0)*5+j] * (u[1] + u[0]),
        d[(e*3+1)*5+j] * (u[2] + u[1]),
        d[(e*3+2)*5+j] * (u[0] + u[2]) };
      G(N[0],j) = G(N[0],j) - f[0] + f[2];
      G(N[1],j) = G(N[1],j) + f[0] - f[1];
      G(N[2],j) = G(N[2],j) + f[1] - f[2];
    }
  }

  // domain edge contributions: edges
  for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
    const auto N = dsupedge[2].data() + e*2;
    const auto d = dsupint[2].data() + e*5;
    tk::real u[] = { U[N[0]], U[N[1]] };<--- Variable 'u' can be declared as const array
    for (std::size_t j=0; j<3; ++j) {
      tk::real f = d[j] * (u[1] + u[0]);
      G(N[0],j) -= f;
      G(N[1],j) += f;
    }
  }

  // boundary integral

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    const auto N = triinpoel.data() + e*3;
    tk::real n[3];
    tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
                  x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
                  n[0], n[1], n[2] );
    auto uA = U[N[0]];
    auto uB = U[N[1]];
    auto uC = U[N[2]];
    auto f = (6.0*uA + uB + uC)/8.0;
    G(N[0],0) += f * n[0];
    G(N[0],1) += f * n[1];
    G(N[0],2) += f * n[2];
    f = (uA + 6.0*uB + uC)/8.0;
    G(N[1],0) += f * n[0];
    G(N[1],1) += f * n[1];
    G(N[1],2) += f * n[2];
    f = (uA + uB + 6.0*uC)/8.0;
    G(N[2],0) += f * n[0];
    G(N[2],1) += f * n[1];
    G(N[2],2) += f * n[2];
  }

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

static tk::real
flux( const tk::Fields& U,
      const tk::Fields& G,
      std::size_t i,
      std::size_t j,
      std::size_t p,
      std::size_t q )
// *****************************************************************************
//! Compute momentum flux over edge of points p-q
//! \param[in] U Velocity vector
//! \param[in] G Velocity gradients
//! \param[in] i Tensor component, 1st index
//! \param[in] j Tensor component, 2nd index
//! \param[in] p Left node index of edge
//! \param[in] q Right node index of edge
//! \return Momentum flux contribution for edge p-q
// *****************************************************************************
{
  auto inv = U(p,i)*U(p,j) + U(q,i)*U(q,j);

  auto eps = std::numeric_limits< tk::real >::epsilon();
  auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
  if (mu < eps) return -inv;

  auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
  if (i == j) {
    vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) + G(q,0) + G(q,4) + G(q,8) );
  }
  return mu*vis - inv;
}

static tk::real
flux( const tk::Fields& U,
      const tk::Fields& G,
      std::size_t i,
      std::size_t j,
      std::size_t p )
// *****************************************************************************
//! Compute momentum flux in point p
//! \param[in] U Velocity vector
//! \param[in] G Velocity gradients
//! \param[in] i Tensor component, 1st index
//! \param[in] j Tensor component, 2nd index
//! \param[in] p Node index of point
//! \return Momentum flux contribution for point p
// *****************************************************************************
{
  auto inv = U(p,i)*U(p,j);

  auto eps = std::numeric_limits< tk::real >::epsilon();
  auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
  if (mu < eps) return -inv;

  auto vis = G(p,i*3+j) + G(p,j*3+i);
  if (i == j) {
    vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
  }
  return mu*vis - inv;
}

void
flux( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
      const std::array< std::vector< tk::real >, 3 >& dsupint,
      const std::array< std::vector< tk::real >, 3 >& coord,
      const std::vector< std::size_t >& triinpoel,
      const tk::Fields& U,
      const tk::Fields& G,
      tk::Fields& F )
// *****************************************************************************
//  Compute momentum flux in all points
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] U Velocity field
//! \param[in] G Velocity gradients, dui/dxj, 9 components
//! \param[in,out] F Momentum flux, Fi = d[ sij - uiuj ] / dxj, where
//!    s_ij = mu[ dui/dxj + duj/dxi - 2/3 duk/dxk delta_ij ]
// *****************************************************************************
{
  Assert( F.nunk() == U.nunk(), "Size mismatch" );
  Assert( F.nprop() == 3, "Size mismatch" );
  Assert( G.nunk() == U.nunk(), "Size mismatch" );
  Assert( G.nprop() == 9, "Size mismatch" );

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // domain integral

  // domain edge contributions: tetrahedron superedges
  for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
    const auto N = dsupedge[0].data() + e*4;
    const auto d = dsupint[0].data();
    for (std::size_t i=0; i<3; ++i) {
      for (std::size_t j=0; j<3; ++j) {
        tk::real f[] = { d[(e*6+0)*5+j] * flux(U,G,i,j,N[1],N[0]),<--- Variable 'f' can be declared as const array
                         d[(e*6+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
                         d[(e*6+2)*5+j] * flux(U,G,i,j,N[0],N[2]),
                         d[(e*6+3)*5+j] * flux(U,G,i,j,N[3],N[0]),
                         d[(e*6+4)*5+j] * flux(U,G,i,j,N[3],N[1]),
                         d[(e*6+5)*5+j] * flux(U,G,i,j,N[3],N[2]) };
        F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
        F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
        F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
        F(N[3],i) = F(N[3],i) + f[3] + f[4] + f[5];
      }
    }
  }

  // domain edge contributions: triangle superedges
  for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
    const auto N = dsupedge[1].data() + e*3;
    const auto d = dsupint[1].data();
    for (std::size_t i=0; i<3; ++i) {
      for (std::size_t j=0; j<3; ++j) {
        tk::real f[] = { d[(e*3+0)*5+j] * flux(U,G,i,j,N[1],N[0]),<--- Variable 'f' can be declared as const array
                         d[(e*3+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
                         d[(e*3+2)*5+j] * flux(U,G,i,j,N[0],N[2]), };
        F(N[0],i) = F(N[0],i) - f[0] + f[2];
        F(N[1],i) = F(N[1],i) + f[0] - f[1];
        F(N[2],i) = F(N[2],i) + f[1] - f[2];
      }
    }
  }

  // domain edge contributions: edges
  for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
    const auto N = dsupedge[2].data() + e*2;
    const auto d = dsupint[2].data() + e*5;
    for (std::size_t i=0; i<3; ++i) {
      for (std::size_t j=0; j<3; ++j) {
        tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
        F(N[0],i) -= f;
        F(N[1],i) += f;
      }
    }
  }

  // boundary integral

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    const auto N = triinpoel.data() + e*3;
    tk::real n[3];
    tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
                  x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
                  n[0], n[1], n[2] );
    for (std::size_t i=0; i<3; ++i) {
      auto fxA = flux(U,G,i,0,N[0]);
      auto fyA = flux(U,G,i,1,N[0]);
      auto fzA = flux(U,G,i,2,N[0]);
      auto fxB = flux(U,G,i,0,N[1]);
      auto fyB = flux(U,G,i,1,N[1]);
      auto fzB = flux(U,G,i,2,N[1]);
      auto fxC = flux(U,G,i,0,N[2]);
      auto fyC = flux(U,G,i,1,N[2]);
      auto fzC = flux(U,G,i,2,N[2]);
      auto fx = (6.0*fxA + fxB + fxC)/8.0;
      auto fy = (6.0*fyA + fyB + fyC)/8.0;
      auto fz = (6.0*fzA + fzB + fzC)/8.0;
      F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
      fx = (fxA + 6.0*fxB + fxC)/8.0;
      fy = (fyA + 6.0*fyB + fyC)/8.0;
      fz = (fzA + 6.0*fzB + fzC)/8.0;
      F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
      fx = (fxA + fxB + 6.0*fxC)/8.0;
      fy = (fyA + fyB + 6.0*fyC)/8.0;
      fz = (fzA + fzB + 6.0*fzC)/8.0;
      F(N[2],i) += fx*n[0] + fy*n[1] + fz*n[2];
    }
  }

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

static void
adv_tg( const tk::real supint[],<--- The function 'adv_tg' is never used.
        const tk::Fields& U,
        const tk::Fields&,
        const std::vector< tk::real >& P,
        const tk::Fields&,
        const std::array< std::vector< tk::real >, 3 >& coord,
        tk::real dt,
        std::size_t p,
        std::size_t q,
        tk::real f[] )
// *****************************************************************************
//! Compute advection fluxes on a single edge using Taylor-Galerkin
//! \param[in] supint Edge integral
//! \param[in] U Velocity and transported scalars at recent time step
//! \param[in] P Pressure
//! \param[in] coord Mesh node coordinates
//! \param[in] dt Physical time step size
//! \param[in] p Left node index of edge
//! \param[in] q Right node index of edge
//! \param[in,out] f Flux computed
// *****************************************************************************
{
  const auto ncomp = U.nprop();
  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  // edge vector
  auto dx = x[p] - x[q];
  auto dy = y[p] - y[q];
  auto dz = z[p] - z[q];
  auto dl = dx*dx + dy*dy + dz*dz;
  dx /= dl;
  dy /= dl;
  dz /= dl;

  // left state
  auto uL = U(p,0);
  auto vL = U(p,1);
  auto wL = U(p,2);
  auto pL = P[p];
  auto dnL = uL*dx + vL*dy + wL*dz;

  // right state
  auto uR = U(q,0);
  auto vR = U(q,1);
  auto wR = U(q,2);
  auto pR = P[q];
  auto dnR = uR*dx + vR*dy + wR*dz;

  auto nx = supint[0];
  auto ny = supint[1];
  auto nz = supint[2];

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // Taylor-Galerkin first half step

  tk::real ue[ncomp];

  // flow
  auto dp = pL - pR;
  ue[0] = 0.5*(uL + uR - dt*(uL*dnL - uR*dnR + dp*dx));
  ue[1] = 0.5*(vL + vR - dt*(vL*dnL - vR*dnR + dp*dy));
  ue[2] = 0.5*(wL + wR - dt*(wL*dnL - wR*dnR + dp*dz));

  // Taylor-Galerkin second half step

  auto uh = ue[0];
  auto vh = ue[1];
  auto wh = ue[2];
  auto ph = (pL + pR)/2.0;
  auto vn = uh*nx + vh*ny + wh*nz;

  // viscosity
  auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();

  // flow
  f[0] = 2.0*(uh*vn + ph*nx) - d*(uR - uL);
  f[1] = 2.0*(vh*vn + ph*ny) - d*(vR - vL);
  f[2] = 2.0*(wh*vn + ph*nz) - d*(wR - wL);

  // artificial viscosity

  const auto stab2 = g_cfg.get< tag::stab2 >();
  if (!stab2) return;

  auto stab2coef = g_cfg.get< tag::stab2coef >();
  auto vnL = uL*nx + vL*ny + wL*nz;
  auto vnR = uR*nx + vR*ny + wR*nz;
  auto sl = std::abs(vnL);
  auto sr = std::abs(vnR);
  auto fw = stab2coef * std::max( sl, sr );

  // flow
  f[0] += fw*(uR - uL);
  f[1] += fw*(vR - vL);
  f[2] += fw*(wR - wL);

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

static void
adv_damp2( const tk::real supint[],<--- The function 'adv_damp2' is never used.
           const tk::Fields& U,
           const tk::Fields&,
           const std::vector< tk::real >& P,
           const tk::Fields&,
           const std::array< std::vector< tk::real >, 3 >&,
           tk::real,
           std::size_t p,
           std::size_t q,
           tk::real f[] )
// *****************************************************************************
//! Compute advection fluxes on a single edge using 2nd-order damping
//! \param[in] supint Edge integral
//! \param[in] U Velocity and transported scalars at recent time step
//! \param[in] P Pressure
//! \param[in] p Left node index of edge
//! \param[in] q Right node index of edge
//! \param[in,out] f Flux computed
// *****************************************************************************
{
  auto nx = supint[0];
  auto ny = supint[1];
  auto nz = supint[2];

  // left state
  auto uL = U(p,0);
  auto vL = U(p,1);
  auto wL = U(p,2);
  auto vnL = uL*nx + vL*ny + wL*nz;

  // right state
  auto uR = U(q,0);
  auto vR = U(q,1);
  auto wR = U(q,2);
  auto vnR = uR*nx + vR*ny + wR*nz;

  // stabilization
  auto aw = std::abs( vnL + vnR ) / 2.0 * tk::length( nx, ny, nz );

  // viscosity
  auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();

  // flow
  auto pf = P[p] + P[q];
  f[0] = uL*vnL + uR*vnR + pf*nx + (aw-d)*(uR-uL);
  f[1] = vL*vnL + vR*vnR + pf*ny + (aw-d)*(vR-vL);
  f[2] = wL*vnL + wR*vnR + pf*nz + (aw-d)*(wR-wL);
}

static void
adv_damp4( const tk::real supint[],<--- The function 'adv_damp4' is never used.
           const tk::Fields& U,
           const tk::Fields& G,
           const std::vector< tk::real >& P,
           const tk::Fields& W,
           const std::array< std::vector< tk::real >, 3 >& coord,
           tk::real,
           std::size_t p,
           std::size_t q,
           tk::real f[] )
// *****************************************************************************
//! Compute advection fluxes on a single edge using 4th-order damping
//! \param[in] supint Edge integral
//! \param[in] U Velocity and transported scalars at recent time step
//! \param[in] G Gradients of velocity and transported scalars
//! \param[in] P Pressure
//! \param[in] W Pressure gradient
//! \param[in] coord Mesh node coordinates
//! \param[in] p Left node index of edge
//! \param[in] q Right node index of edge
//! \param[in,out] f Flux computed
// *****************************************************************************
{
  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  // edge vector
  auto dx = x[p] - x[q];
  auto dy = y[p] - y[q];
  auto dz = z[p] - z[q];

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  tk::real uL[] = { U(p,0), U(p,1), U(p,2), P[p] };
  tk::real uR[] = { U(q,0), U(q,1), U(q,2), P[q] };
  tk::real gL[] = { G(p,0), G(p,1), G(p,2),<--- Variable 'gL' can be declared as const array
                    G(p,3), G(p,4), G(p,5),
                    G(p,6), G(p,7), G(p,8),
                    W(p,0), W(p,1), W(p,2) };
  tk::real gR[] = { G(q,0), G(q,1), G(q,2),<--- Variable 'gR' can be declared as const array
                    G(q,3), G(q,4), G(q,5),
                    G(q,6), G(q,7), G(q,8),
                    W(q,0), W(q,1), W(q,2) };

  // MUSCL reconstruction in edge-end points
  for (std::size_t c=0; c<4; ++c) {
    auto g = c*3;
    auto g1 = gL[g+0]*dx + gL[g+1]*dy + gL[g+2]*dz;
    auto g2 = gR[g+0]*dx + gR[g+1]*dy + gR[g+2]*dz;
    auto delta2 = uR[c] - uL[c];
    auto delta1 = 2.0 * g1 - delta2;
    auto delta3 = 2.0 * g2 - delta2;

    // van Leer limiter
    auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
    auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
    auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
    auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
    auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
    auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
    auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
    auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
    // update unknowns with reconstructed unknowns
    uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
                   delta2*(1.0+muscl_const)*phi_L_inv);
    uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
                   delta2*(1.0+muscl_const)*phi_R_inv);
  }

  auto nx = supint[0];
  auto ny = supint[1];
  auto nz = supint[2];

  // normal velocities
  auto vnL = uL[0]*nx + uL[1]*ny + uL[2]*nz;
  auto vnR = uR[0]*nx + uR[1]*ny + uR[2]*nz;

  // stabilization
  auto aw = std::abs( vnL + vnR ) / 2.0 * tk::length( nx, ny, nz );

  // viscosity
  auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();

  // flow
  auto pf = uL[3] + uR[3];
  f[0] = uL[0]*vnL + uR[0]*vnR + pf*nx + aw*(uR[0]-uL[0]) - d*(U(q,0)-U(p,0));
  f[1] = uL[1]*vnL + uR[1]*vnR + pf*ny + aw*(uR[1]-uL[1]) - d*(U(q,1)-U(p,1));
  f[2] = uL[2]*vnL + uR[2]*vnR + pf*nz + aw*(uR[2]-uL[2]) - d*(U(q,2)-U(p,2));

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

static void
adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     const std::array< std::vector< tk::real >, 3 >& dsupint,
     const std::array< std::vector< tk::real >, 3 >& coord,
     const std::vector< std::size_t >& triinpoel,
     tk::real dt,
     const tk::Fields& U,
     const tk::Fields& G,
     const std::vector< tk::real >& P,
     const tk::Fields& W,
     // cppcheck-suppress constParameter
     tk::Fields& R )
// *****************************************************************************
//! Add advection to rhs
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] dt Physical time step size
//! \param[in] U Velocity and transported scalars at recent time step
//! \param[in] G Gradients of velocity and transported scalars
//! \param[in] P Pressure
//! \param[in] W Pressure gradient
//! \param[in,out] R Right-hand side vector added to
// *****************************************************************************
{
  // configure advection stabilization
  auto adv = [](){
    const auto& flux = g_cfg.get< tag::flux >();
         if (flux == "tg")    return adv_tg;
    else if (flux == "damp2") return adv_damp2;
    else if (flux == "damp4") return adv_damp4;
    else Throw( "Flux not correctly configured" );
  }();

  auto ncomp = U.nprop();

  #if defined(__clang__)
    #pragma clang diagnostic push
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wvla-extension"
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wvla"
  #endif

  // domain integral

  // domain edge contributions: tetrahedron superedges
  for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
    const auto N = dsupedge[0].data() + e*4;
    tk::real f[6][ncomp];
    const auto d = dsupint[0].data();
    adv( d+(e*6+0)*5, U, G, P, W, coord, dt, N[0], N[1], f[0] );
    adv( d+(e*6+1)*5, U, G, P, W, coord, dt, N[1], N[2], f[1] );
    adv( d+(e*6+2)*5, U, G, P, W, coord, dt, N[2], N[0], f[2] );
    adv( d+(e*6+3)*5, U, G, P, W, coord, dt, N[0], N[3], f[3] );
    adv( d+(e*6+4)*5, U, G, P, W, coord, dt, N[1], N[3], f[4] );
    adv( d+(e*6+5)*5, U, G, P, W, coord, dt, N[2], N[3], f[5] );
    for (std::size_t c=0; c<ncomp; ++c) {
      R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
      R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
      R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
      R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
    }
  }

  // domain edge contributions: triangle superedges
  for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
    const auto N = dsupedge[1].data() + e*3;
    tk::real f[3][ncomp];
    const auto d = dsupint[1].data();
    adv( d+(e*3+0)*5, U, G, P, W, coord, dt, N[0], N[1], f[0] );
    adv( d+(e*3+1)*5, U, G, P, W, coord, dt, N[1], N[2], f[1] );
    adv( d+(e*3+2)*5, U, G, P, W, coord, dt, N[2], N[0], f[2] );
    for (std::size_t c=0; c<ncomp; ++c) {
      R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
      R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
      R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
    }
  }

  // domain edge contributions: edges
  for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
    const auto N = dsupedge[2].data() + e*2;
    tk::real f[ncomp];
    const auto d = dsupint[2].data();
    adv( d+e*5, U, G, P, W, coord, dt, N[0], N[1], f );
    for (std::size_t c=0; c<ncomp; ++c) {
      R(N[0],c) -= f[c];
      R(N[1],c) += f[c];
    }
  }

  // boundary integral

  const auto& x = coord[0];
  const auto& y = coord[1];
  const auto& z = coord[2];

  for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
    const auto N = triinpoel.data() + e*3;
    tk::real n[3];
    tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
                  x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
                  n[0], n[1], n[2] );
    tk::real f[ncomp][3];

    auto u = U(N[0],0);
    auto v = U(N[0],1);
    auto w = U(N[0],2);
    auto p = P[N[0]];
    auto vn = n[0]*u + n[1]*v + n[2]*w;
    f[0][0] = u*vn + p*n[0];
    f[1][0] = v*vn + p*n[1];
    f[2][0] = w*vn + p*n[2];
    for (std::size_t c=3; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;

    u = U(N[1],0);
    v = U(N[1],1);
    w = U(N[1],2);
    p = P[N[1]];
    vn = n[0]*u + n[1]*v + n[2]*w;
    f[0][1] = u*vn + p*n[0];
    f[1][1] = v*vn + p*n[1];
    f[2][1] = w*vn + p*n[2];
    for (std::size_t c=3; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;

    u = U(N[2],0);
    v = U(N[2],1);
    w = U(N[2],2);
    p = P[N[2]];
    vn = n[0]*u + n[1]*v + n[2]*w;
    f[0][2] = u*vn + p*n[0];
    f[1][2] = v*vn + p*n[1];
    f[2][2] = w*vn + p*n[2];
    for (std::size_t c=3; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;

    for (std::size_t c=0; c<ncomp; ++c) {
      R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
      R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
      R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
    }
  }

  #if defined(__clang__)
    #pragma clang diagnostic pop
  #elif defined(STRICT_GNUC)
    #pragma GCC diagnostic pop
  #endif
}

void
rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
     const std::array< std::vector< tk::real >, 3 >& dsupint,
     const std::array< std::vector< tk::real >, 3 >& coord,
     const std::vector< std::size_t >& triinpoel,
     tk::real dt,
     const std::vector< tk::real >& P,
     const tk::Fields& U,
     const tk::Fields& G,
     const tk::Fields& W,
     tk::Fields& R )
// *****************************************************************************
//  Compute right hand side
//! \param[in] dsupedge Domain superedges
//! \param[in] dsupint Domain superedge integrals
//! \param[in] coord Mesh node coordinates
//! \param[in] triinpoel Boundary face connectivity
//! \param[in] dt Physical time step size
//! \param[in] P Pressure
//! \param[in] U Solution vector of primitive variables at recent time step
//! \param[in] G Gradients of velocity and transported scalars
//! \param[in] W Pressure gradient
//! \param[in,out] R Right-hand side vector computed
// *****************************************************************************
{
  Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
          "vector at recent time step incorrect" );
  Assert( R.nunk() == coord[0].size(),
          "Number of unknowns and/or number of components in right-hand "
          "side vector incorrect" );

  R.fill( 0.0 );
  adv( dsupedge, dsupint, coord, triinpoel, dt, U, G, P, W, R );
}

} // chorin::