Quinoa all test code coverage report
Current view: top level - DiffEq/Dirichlet - Dirichlet.hpp (source / functions) Hit Total Coverage
Commit: Quinoa_v0.3-957-gb4f0efae0 Lines: 28 28 100.0 %
Date: 2021-11-11 18:25:50 Functions: 3 24 12.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 18 28 64.3 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/DiffEq/Dirichlet/Dirichlet.hpp
       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     Dirichlet SDE
       9                 :            :   \details   This file implements the time integration of a system of stochastic
      10                 :            :     differential equations (SDEs), whose invariant is the [Dirichlet
      11                 :            :     distribution](http://en.wikipedia.org/wiki/Dirichlet_distribution).
      12                 :            : 
      13                 :            :     In a nutshell, the equation integrated governs a set of scalars,
      14                 :            :     \f$0\!\le\!Y_\alpha\f$, \f$\alpha\!=\!1,\dots,N-1\f$,
      15                 :            :     \f$\sum_{\alpha=1}^{N-1}Y_\alpha\!\le\!1\f$, as
      16                 :            : 
      17                 :            :     @m_class{m-show-m}
      18                 :            : 
      19                 :            :     \f[
      20                 :            :        \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
      21                 :            :        (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
      22                 :            :        Y_N}\mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N-1
      23                 :            :     \f]
      24                 :            : 
      25                 :            :     @m_class{m-hide-m}
      26                 :            : 
      27                 :            :     \f[
      28                 :            :        \begin{split}
      29                 :            :        \mathrm{d}Y_\alpha(t) = \frac{b_\alpha}{2}\big[S_\alpha Y_N -
      30                 :            :        (1-S_\alpha)Y_\alpha\big]\mathrm{d}t + \sqrt{\kappa_\alpha Y_\alpha
      31                 :            :        Y_N}\mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N-1
      32                 :            :        \end{split}
      33                 :            :     \f]
      34                 :            : 
      35                 :            :     with parameter vectors \f$b_\alpha > 0\f$, \f$\kappa_\alpha > 0\f$, and \f$0
      36                 :            :     < S_\alpha < 1\f$, and \f$Y_N = 1-\sum_{\alpha=1}^{N-1}Y_\alpha\f$. Here
      37                 :            :     \f$\mathrm{d}W_\alpha(t)\f$ is an isotropic vector-valued [Wiener
      38                 :            :     process](http://en.wikipedia.org/wiki/Wiener_process) with independent
      39                 :            :     increments. The invariant distribution is the Dirichlet distribution,
      40                 :            :     provided the parameters of the drift and diffusion terms satisfy
      41                 :            :     \f[
      42                 :            :       (1-S_1) b_1 / \kappa_1 = \dots = (1-S_{N-1}) b_{N-1} / \kappa_{N-1}.
      43                 :            :     \f]
      44                 :            :     To keep the invariant distribution Dirichlet, the above constraint on the
      45                 :            :     coefficients must be satisfied. For more details on the Dirichlet SDE, see
      46                 :            :     https://doi.org/10.1155/2013/842981.
      47                 :            : */
      48                 :            : // *****************************************************************************
      49                 :            : #ifndef Dirichlet_h
      50                 :            : #define Dirichlet_h
      51                 :            : 
      52                 :            : #include <vector>
      53                 :            : #include <cmath>
      54                 :            : 
      55                 :            : #include "InitPolicy.hpp"
      56                 :            : #include "DirichletCoeffPolicy.hpp"
      57                 :            : #include "RNG.hpp"
      58                 :            : #include "Particles.hpp"
      59                 :            : 
      60                 :            : namespace walker {
      61                 :            : 
      62                 :            : extern ctr::InputDeck g_inputdeck;
      63                 :            : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
      64                 :            : 
      65                 :            : //! \brief Dirichlet SDE used polymorphically with DiffEq
      66                 :            : //! \details The template arguments specify policies and are used to configure
      67                 :            : //!   the behavior of the class. The policies are:
      68                 :            : //!   - Init - initialization policy, see DiffEq/InitPolicy.h
      69                 :            : //!   - Coefficients - coefficients policy, see DiffEq/DirCoeffPolicy.h
      70                 :            : template< class Init, class Coefficients >
      71                 :            : class Dirichlet {
      72                 :            : 
      73                 :            :   private:
      74                 :            :     using ncomp_t = tk::ctr::ncomp_t;
      75                 :            : 
      76                 :            :   public:
      77                 :            :     //! \brief Constructor
      78                 :            :     //! \param[in] c Index specifying which Dirichlet SDE to construct. There
      79                 :            :     //!   can be multiple dirichlet ... end blocks in a control file. This index
      80                 :            :     //!   specifies which Dirichlet SDE to instantiate. The index corresponds to
      81                 :            :     //!   the order in which the dirichlet ... end blocks are given the control
      82                 :            :     //!   file.
      83                 :         22 :     explicit Dirichlet( ncomp_t c ) :
      84                 :            :       m_c( c ),
      85                 :            :       m_depvar(
      86                 :         22 :         g_inputdeck.get< tag::param, tag::dirichlet, tag::depvar >().at(c) ),
      87                 :            :       m_ncomp(
      88                 :         22 :         g_inputdeck.get< tag::component >().get< tag::dirichlet >().at(c) ),
      89                 :            :       m_offset(
      90                 :         22 :         g_inputdeck.get< tag::component >().offset< tag::dirichlet >(c) ),
      91         [ +  - ]:         22 :       m_rng( g_rng.at( tk::ctr::raw(
      92                 :         22 :         g_inputdeck.get< tag::param, tag::dirichlet, tag::rng >().at(c) ) ) ),
      93                 :            :       m_b(),
      94                 :            :       m_S(),
      95                 :            :       m_k(),
      96                 :         22 :       coeff( m_ncomp,
      97         [ +  - ]:         22 :              g_inputdeck.get< tag::param, tag::dirichlet, tag::b >().at(c),
      98         [ +  - ]:         22 :              g_inputdeck.get< tag::param, tag::dirichlet, tag::S >().at(c),
      99         [ +  - ]:         22 :              g_inputdeck.get< tag::param, tag::dirichlet, tag::kappa >().at(c),
     100         [ +  - ]:        132 :              m_b, m_S, m_k ) {}
     101                 :            : 
     102                 :            :     //! Initalize SDE, prepare for time integration
     103                 :            :     //! \param[in] stream Thread (or more precisely stream) ID 
     104                 :            :     //! \param[in,out] particles Array of particle properties 
     105                 :         94 :     void initialize( int stream, tk::Particles& particles ) {
     106                 :            :       //! Set initial conditions using initialization policy
     107                 :            :       Init::template
     108                 :            :         init< tag::dirichlet >
     109                 :         94 :             ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
     110                 :         94 :     }
     111                 :            : 
     112                 :            :     //! \brief Advance particles according to the Dirichlet SDE
     113                 :            :     //! \param[in,out] particles Array of particle properties
     114                 :            :     //! \param[in] stream Thread (or more precisely stream) ID
     115                 :            :     //! \param[in] dt Time step size
     116                 :     263200 :     void advance( tk::Particles& particles,
     117                 :            :                   int stream,
     118                 :            :                   tk::real dt,
     119                 :            :                   tk::real,
     120                 :            :                   const std::map< tk::ctr::Product, tk::real >& )
     121                 :            :     {
     122                 :     263200 :       const auto npar = particles.nunk();
     123         [ +  + ]:  224263200 :       for (auto p=decltype(npar){0}; p<npar; ++p) {
     124                 :            :         // Compute Nth scalar
     125         [ +  - ]:  224000000 :         tk::real yn = 1.0 - particles(p, 0, m_offset);
     126         [ +  + ]:  448000000 :         for (ncomp_t i=1; i<m_ncomp; ++i)
     127         [ +  - ]:  224000000 :           yn -= particles( p, i, m_offset );
     128                 :            : 
     129                 :            :         // Generate Gaussian random numbers with zero mean and unit variance
     130         [ +  - ]:  448000000 :         std::vector< tk::real > dW( m_ncomp );
     131         [ +  - ]:  224000000 :         m_rng.gaussian( stream, m_ncomp, dW.data() );
     132                 :            : 
     133                 :            :         // Advance first m_ncomp (K=N-1) scalars
     134         [ +  + ]:  672000000 :         for (ncomp_t i=0; i<m_ncomp; ++i) {
     135         [ +  - ]:  448000000 :           tk::real& par = particles( p, i, m_offset );
     136                 :  448000000 :           tk::real d = m_k[i] * par * yn * dt;
     137         [ +  + ]:  448000000 :           d = (d > 0.0 ? std::sqrt(d) : 0.0);
     138                 :  448000000 :           par += 0.5*m_b[i]*( m_S[i]*yn - (1.0-m_S[i]) * par )*dt + d*dW[i];
     139                 :            :         }
     140                 :            :       }
     141                 :     263200 :     }
     142                 :            : 
     143                 :            :   private:
     144                 :            :     const ncomp_t m_c;                  //!< Equation system index
     145                 :            :     const char m_depvar;                //!< Dependent variable
     146                 :            :     const ncomp_t m_ncomp;              //!< Number of components
     147                 :            :     const ncomp_t m_offset;             //!< Offset SDE operates from
     148                 :            :     const tk::RNG& m_rng;               //!< Random number generator
     149                 :            : 
     150                 :            :     //! Coefficients
     151                 :            :     std::vector< kw::sde_b::info::expect::type > m_b;
     152                 :            :     std::vector< kw::sde_S::info::expect::type > m_S;
     153                 :            :     std::vector< kw::sde_kappa::info::expect::type > m_k;
     154                 :            : 
     155                 :            :     //! Coefficients policy
     156                 :            :     Coefficients coeff;
     157                 :            : };
     158                 :            : 
     159                 :            : } // walker::
     160                 :            : 
     161                 :            : #endif // Dirichlet_h

Generated by: LCOV version 1.14