Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Transport/Physics/CGAdvDiff.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Physics policy for the transport equations using continuous
       9                 :            :              Galerkin (CG)
      10                 :            :   \details   This file defines a Physics policy class for the transport
      11                 :            :     equations, defined in PDE/Transport/CGTransport.h implementing
      12                 :            :     node-centered continuous Galerkin (CG) discretizations.
      13                 :            :     See PDE/Transport/Physics/CG.h for general requirements on Physics policy
      14                 :            :     classes for cg::Transport.
      15                 :            : */
      16                 :            : // *****************************************************************************
      17                 :            : 
      18                 :            : #include "CGAdvDiff.hpp"
      19                 :            : #include "Inciter/InputDeck/InputDeck.hpp"
      20                 :            : 
      21                 :            : namespace inciter {
      22                 :            : 
      23                 :            : extern ctr::InputDeck g_inputdeck;
      24                 :            : 
      25                 :            : } // ::inciter
      26                 :            : 
      27                 :            : using inciter::cg::TransportPhysicsAdvDiff;
      28                 :            : 
      29                 :            : void
      30                 :     875420 : TransportPhysicsAdvDiff::diffusionRhs(
      31                 :            :   ncomp_t system,
      32                 :            :   ncomp_t ncomp,
      33                 :            :   tk::real deltat,
      34                 :            :   tk::real J,
      35                 :            :   const std::array< std::array< tk::real, 3 >, 4 >& grad,
      36                 :            :   const std::array< std::size_t, 4 >& N,
      37                 :            :   const std::vector< std::array< tk::real, 4 > >& u,
      38                 :            :   const std::vector< const tk::real* >& r,
      39                 :            :   tk::Fields& R ) const
      40                 :            : // *****************************************************************************
      41                 :            : //  Add diffusion contribution to rhs
      42                 :            : //! \param[in] system Equation system index, i.e., which transport equation
      43                 :            : //!   system we operate on among the systems of PDEs
      44                 :            : //! \param[in] ncomp Number of components in this PDE
      45                 :            : //! \param[in] deltat Size of time step
      46                 :            : //! \param[in] J Element Jacobi determinant
      47                 :            : //! \param[in] grad Shape function derivatives, nnode*ndim [4][3]
      48                 :            : //! \param[in] N Element node indices
      49                 :            : //! \param[in] u Solution at element nodes at recent time step
      50                 :            : //! \param[in] r Pointers to right hand side at component and offset
      51                 :            : //! \param[in,out] R Right-hand side vector contributing to
      52                 :            : // *****************************************************************************
      53                 :            : {
      54                 :            :   // diffusivities for all components
      55                 :            :   const auto& diff =
      56                 :     875420 :     g_inputdeck.get< tag::param, eq, tag::diffusivity >()[ system ];
      57                 :            : 
      58                 :            :   // add diffusion contribution to right hand side
      59                 :     875420 :   const auto d = deltat * J/6.0;
      60         [ +  + ]:    4377100 :   for (std::size_t a=0; a<4; ++a)
      61         [ +  + ]:    7503600 :     for (ncomp_t c=0; c<ncomp; ++c)
      62         [ +  + ]:   16007680 :       for (std::size_t k=0; k<3; ++k) {
      63                 :   12005760 :         const auto D = diff[ 3*c+k ];
      64         [ +  + ]:   60028800 :           for (std::size_t b=0; b<4; ++b)
      65                 :   48023040 :             R.var(r[c],N[a]) -= d * D * grad[a][k] * grad[b][k] * u[c][b];
      66                 :            :       }
      67                 :     875420 : }
      68                 :            : 
      69                 :            : tk::real
      70                 :     875420 : TransportPhysicsAdvDiff::diffusion_dt(
      71                 :            :   ncomp_t system,
      72                 :            :   ncomp_t ncomp,
      73                 :            :   tk::real L,
      74                 :            :   const std::vector< std::array< tk::real, 4 > >& ) const
      75                 :            : // *****************************************************************************
      76                 :            : //! Compute the minimum time step size based on the diffusion
      77                 :            : //! \param[in] system Equation system index, i.e., which transport equation
      78                 :            : //!   system we operate on among the systems of PDEs
      79                 :            : //! \param[in] ncomp Number of components in this PDE
      80                 :            : //! \param[in] L Characteristic length scale
      81                 :            : //! \return Minimum time step size based on diffusion
      82                 :            : // *****************************************************************************
      83                 :            : {
      84                 :            :   // diffusivities for all components
      85                 :            :   const auto& df =
      86                 :     875420 :     g_inputdeck.get< tag::param, eq, tag::diffusivity >()[ system ];
      87                 :            : 
      88                 :            :   // compute the minimum diffusion time step size across the four nodes
      89                 :            :   tk::real mindt = std::numeric_limits< tk::real >::max();
      90         [ +  + ]:    1875900 :   for (ncomp_t c=0; c<ncomp; ++c) {
      91                 :    1000480 :     const auto di = 3*c;
      92 [ +  + ][ +  + ]:    1125540 :     const auto d = std::max( df[di+2], std::max( df[di+0], df[di+1] ) );
      93                 :    1000480 :     const auto dt = L * L / (2.0*d);  // dt ~ dx^2/(2D)
      94         [ +  + ]:    1000480 :     if (dt < mindt) mindt = dt;
      95                 :            :   }
      96                 :            : 
      97                 :     875420 :   return mindt;
      98                 :            : }
      99                 :            :