Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/DiffEq/Position/Position.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 A position model for Lagrangian particles
9 : : \details This file implements the time integration of a system of
10 : : deterministic or stochastic differential equations to model fluctuating
11 : : particle postitions.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef Position_h
15 : : #define Position_h
16 : :
17 : : #include <array>
18 : :
19 : : #include "InitPolicy.hpp"
20 : : #include "PositionCoeffPolicy.hpp"
21 : : #include "RNG.hpp"
22 : : #include "Particles.hpp"
23 : : #include "CoupledEq.hpp"
24 : :
25 : : namespace walker {
26 : :
27 : : extern ctr::InputDeck g_inputdeck;
28 : : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
29 : :
30 : : //! \brief Position equation used polymorphically with DiffEq
31 : : //! \details The template arguments specify policies and are used to configure
32 : : //! the behavior of the class. The policies are:
33 : : //! - Init - initialization policy, see DiffEq/InitPolicy.h
34 : : //! - Coefficients - coefficients policy, see DiffEq/PositionCoeffPolicy.h
35 : : template< class Init, class Coefficients >
36 : : class Position {
37 : :
38 : : private:
39 : : using ncomp_t = tk::ctr::ncomp_t;
40 : : using eq = tag::position;
41 : :
42 : : public:
43 : : //! \brief Constructor
44 : : //! \param[in] c Index specifying which system of position SDEs to construct
45 : : //! There can be multiple position ... end blocks in a control file. This
46 : : //! index specifies which position SDE system to instantiate. The index
47 : : //! corresponds to the order in which the position ... end blocks are
48 : : //! given the control file.
49 : 18 : explicit Position( ncomp_t c ) :
50 : : m_c( c ),
51 : : m_depvar( g_inputdeck.get< tag::param, eq, tag::depvar >().at(c) ),
52 : : m_ncomp( g_inputdeck.get< tag::component, eq >().at(c) ),
53 : : m_offset(
54 : 18 : g_inputdeck.get< tag::component >().offset< eq >(c) ),
55 : 18 : m_rng( g_rng.at( tk::ctr::raw(
56 : : g_inputdeck.get< tag::param, eq, tag::rng >().at(c) ) ) ),
57 : : m_velocity_coupled( coupled< eq, tag::velocity >( c ) ),
58 : : m_velocity_depvar( depvar< eq, tag::velocity >( c ) ),
59 : 18 : m_velocity_offset( offset< eq, tag::velocity, tag::velocity_id >( c ) ),
60 [ - + ][ - + ]: 36 : m_coeff( m_dU )
[ - + ][ - + ]
61 : : {
62 : : // Zero prescribed mean velocity gradient if full variable is solved for
63 [ - + ]: 18 : if (g_inputdeck.get< tag::param, eq, tag::solve >().at(c) ==
64 : : ctr::DepvarType::FULLVAR) {
65 : : m_dU.fill( 0.0 );
66 : : }
67 : : Assert( m_ncomp == 3, "Position eq number of components must be 3" );
68 : 18 : }
69 : :
70 : : //! Initalize SDE, prepare for time integration
71 : : //! \param[in] stream Thread (or more precisely stream) ID
72 : : //! \param[in,out] particles Array of particle properties
73 : : void initialize( int stream, tk::Particles& particles ) {
74 : : //! Set initial conditions using initialization policy
75 : : Init::template init< eq >
76 : 90 : ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
77 : : }
78 : :
79 : : //! \brief Advance particles according to the system of position SDEs
80 : : //! \param[in,out] particles Array of particle properties
81 : : //! \param[in] dt Time step size
82 : 6750 : void advance( tk::Particles& particles,
83 : : int,
84 : : tk::real dt,
85 : : tk::real,
86 : : const std::map< tk::ctr::Product, tk::real >& )
87 : : {
88 : : const auto npar = particles.nunk();
89 [ + + ]: 31506750 : for (auto p=decltype(npar){0}; p<npar; ++p) {
90 : : // Access particle velocity
91 : 31500000 : tk::real u = particles( p, 0, m_velocity_offset );
92 : 31500000 : tk::real v = particles( p, 1, m_velocity_offset );
93 : 31500000 : tk::real w = particles( p, 2, m_velocity_offset );
94 : : // Advance all particle positions
95 : 31500000 : tk::real& Xp = particles( p, 0, m_offset );
96 : : tk::real& Yp = particles( p, 1, m_offset );
97 : : tk::real& Zp = particles( p, 2, m_offset );
98 : : // Advance particle position
99 : 31500000 : Xp += (m_dU[0]*Xp + m_dU[1]*Yp + m_dU[2]*Zp + u)*dt;
100 : 31500000 : Yp += (m_dU[3]*Xp + m_dU[4]*Yp + m_dU[5]*Zp + v)*dt;
101 : 31500000 : Zp += (m_dU[6]*Xp + m_dU[7]*Yp + m_dU[8]*Zp + w)*dt;
102 : : }
103 : 6750 : }
104 : :
105 : : private:
106 : : const ncomp_t m_c; //!< Equation system index
107 : : const char m_depvar; //!< Dependent variable
108 : : const ncomp_t m_ncomp; //!< Number of components
109 : : const ncomp_t m_offset; //!< Offset SDE operates from
110 : : const tk::RNG& m_rng; //!< Random number generator
111 : :
112 : : const bool m_velocity_coupled; //!< True if coupled to velocity
113 : : const char m_velocity_depvar; //!< Coupled velocity dependent variable
114 : : const ncomp_t m_velocity_offset; //!< Offset of coupled velocity eq
115 : :
116 : : //! Coefficients policy
117 : : Coefficients m_coeff;
118 : :
119 : : //! (Optionally) prescribed mean velocity gradient
120 : : std::array< tk::real, 9 > m_dU;
121 : : };
122 : :
123 : : } // walker::
124 : :
125 : : #endif // Position_h
|
