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                 :          0 : TransportPhysicsAdvDiff::diffusionRhs(
      31                 :            :   ncomp_t ncomp,
      32                 :            :   tk::real deltat,
      33                 :            :   tk::real J,
      34                 :            :   const std::array< std::array< tk::real, 3 >, 4 >& grad,
      35                 :            :   const std::array< std::size_t, 4 >& N,
      36                 :            :   const std::vector< std::array< tk::real, 4 > >& u,
      37                 :            :   const std::vector< const tk::real* >& r,
      38                 :            :   tk::Fields& R ) const
      39                 :            : // *****************************************************************************
      40                 :            : //  Add diffusion contribution to rhs
      41                 :            : //! \param[in] ncomp Number of components in this PDE
      42                 :            : //! \param[in] deltat Size of time step
      43                 :            : //! \param[in] J Element Jacobi determinant
      44                 :            : //! \param[in] grad Shape function derivatives, nnode*ndim [4][3]
      45                 :            : //! \param[in] N Element node indices
      46                 :            : //! \param[in] u Solution at element nodes at recent time step
      47                 :            : //! \param[in] r Pointers to right hand side at component
      48                 :            : //! \param[in,out] R Right-hand side vector contributing to
      49                 :            : // *****************************************************************************
      50                 :            : {
      51                 :            :   // diffusivities for all components
      52         [ -  - ]:          0 :   auto diff = g_inputdeck.get< eq, tag::diffusivity >();
      53                 :            : 
      54                 :            :   // add diffusion contribution to right hand side
      55                 :          0 :   const auto d = deltat * J/6.0;
      56         [ -  - ]:          0 :   for (std::size_t a=0; a<4; ++a)
      57         [ -  - ]:          0 :     for (ncomp_t c=0; c<ncomp; ++c)
      58         [ -  - ]:          0 :       for (std::size_t k=0; k<3; ++k) {
      59                 :          0 :         const auto D = diff[ 3*c+k ];
      60         [ -  - ]:          0 :           for (std::size_t b=0; b<4; ++b)
      61         [ -  - ]:          0 :             R.var(r[c],N[a]) -= d * D * grad[a][k] * grad[b][k] * u[c][b];
      62                 :            :       }
      63                 :          0 : }
      64                 :            : 
      65                 :            : tk::real
      66                 :          0 : TransportPhysicsAdvDiff::diffusion_dt(
      67                 :            :   ncomp_t ncomp,
      68                 :            :   tk::real L,
      69                 :            :   const std::vector< std::array< tk::real, 4 > >& ) const
      70                 :            : // *****************************************************************************
      71                 :            : //! Compute the minimum time step size based on the diffusion
      72                 :            : //! \param[in] ncomp Number of components in this PDE
      73                 :            : //! \param[in] L Characteristic length scale
      74                 :            : //! \return Minimum time step size based on diffusion
      75                 :            : // *****************************************************************************
      76                 :            : {
      77                 :            :   // diffusivities for all components
      78         [ -  - ]:          0 :   auto df = g_inputdeck.get< eq, tag::diffusivity >();
      79                 :            : 
      80                 :            :   // compute the minimum diffusion time step size across the four nodes
      81                 :          0 :   tk::real mindt = std::numeric_limits< tk::real >::max();
      82         [ -  - ]:          0 :   for (ncomp_t c=0; c<ncomp; ++c) {
      83                 :          0 :     const auto di = 3*c;
      84                 :          0 :     const auto d = std::max( df[di+2], std::max( df[di+0], df[di+1] ) );
      85                 :          0 :     const auto dt = L * L / (2.0*d);  // dt ~ dx^2/(2D)
      86         [ -  - ]:          0 :     if (dt < mindt) mindt = dt;
      87                 :            :   }
      88                 :            : 
      89                 :          0 :   return mindt;
      90                 :            : }
      91                 :            :