Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Walker/Integrator.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 Integrator advances differential equations
9 : : \details Integrator advances differential equations. There are a potentially
10 : : large number of Integrator Charm++ chares created by Distributor. Each
11 : : integrator gets a chunk of the full load and does the same: initializes and
12 : : advances multiple ordinary or stochastic differential equations in time.
13 : : Note that there is no spatial dependence, these equations describe spatially
14 : : homogeneous processes.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "Integrator.hpp"
19 : : #include "Collector.hpp"
20 : :
21 : : namespace walker {
22 : :
23 : : extern std::vector< DiffEq > g_diffeqs;
24 : :
25 : : }
26 : :
27 : : using walker::Integrator;
28 : :
29 : 956 : Integrator::Integrator( CProxy_Distributor hostproxy,
30 : : CProxy_Collector collproxy,
31 : : tk::CProxy_ParticleWriter particlewriterproxy,
32 : 956 : uint64_t npar ) :
33 : : m_host( hostproxy ),
34 : : m_coll( collproxy ),
35 : : m_particlewriter( particlewriterproxy ),
36 : : m_particles( npar, g_inputdeck.get< tag::component >().nprop() ),
37 : : m_stat( m_particles,
38 [ - - ]: 0 : g_inputdeck.get< tag::component >().offsetmap( g_inputdeck ),
39 : : g_inputdeck.get< tag::stat >(),
40 : : g_inputdeck.get< tag::pdf >(),
41 : : g_inputdeck.get< tag::discr, tag::binsize >() ),
42 : : m_dt( 0.0 ),
43 : : m_t( 0.0 ),
44 : : m_it( 0 ),
45 [ + - ][ + - ]: 2868 : m_itp( 0 )
[ + - ][ + - ]
[ + - ]
46 : : // *****************************************************************************
47 : : // Constructor
48 : : //! \param[in] hostproxy Host proxy to call back to
49 : : //! \param[in] collproxy Collector proxy to send results to
50 : : //! \param[in] particlewriterproxy Particle writer proxy to use
51 : : //! \param[in] npar Number of particles this integrator advances
52 : : // *****************************************************************************
53 : : {
54 : : // register with the local branch of the statistics collector
55 : : m_coll.ckLocalBranch()->checkin();
56 : : // Tell the Charm++ runtime system to call back to
57 : : // Distributor::registered() once all Integrator chares have registered
58 : : // themselves, i.e., checked in, with their local branch of the statistics
59 : : // merger group, Collector. The reduction is done via creating a callback
60 : : // that invokes the typed reduction client, where m_host is the proxy
61 : : // on which the reduction target method, registered(), is called upon
62 : : // completion of the reduction.
63 [ + - ][ + - ]: 1912 : contribute( CkCallback(CkReductionTarget(Distributor, registered), m_host) );
64 : 956 : }
65 : :
66 : : void
67 : 956 : Integrator::setup( tk::real dt,
68 : : tk::real t,
69 : : uint64_t it,
70 : : const std::map< tk::ctr::Product, tk::real >& moments )
71 : : // *****************************************************************************
72 : : // Perform setup: set initial conditions and advance a time step
73 : : //! \param[in] dt Size of time step
74 : : //! \param[in] t Physical time
75 : : //! \param[in] it Iteration count
76 : : //! \param[in] moments Map of statistical moments
77 : : // *****************************************************************************
78 : : {
79 : 956 : ic(); // set initial conditions for all equations
80 : 956 : advance( dt, t, it, moments ); // start time stepping all equations
81 : 956 : }
82 : :
83 : : void
84 : 956 : Integrator::ic()
85 : : // *****************************************************************************
86 : : // Set initial conditions
87 : : // *****************************************************************************
88 : : {
89 [ + + ]: 2092 : for (const auto& eq : g_diffeqs) eq.initialize( CkMyPe(), m_particles );
90 : 956 : }
91 : :
92 : : void
93 : 5509065 : Integrator::advance( tk::real dt,
94 : : tk::real t,
95 : : uint64_t it,
96 : : const std::map< tk::ctr::Product, tk::real >& moments )
97 : : // *****************************************************************************
98 : : // Advance all particles owned by this integrator
99 : : //! \param[in] dt Size of time step
100 : : //! \param[in] t Physical time
101 : : //! \param[in] it Iteration count
102 : : //! \param[in] moments Map of statistical moments
103 : : // *****************************************************************************
104 : : {
105 : : // Advance all equations one step in time. At the 0th iteration skip advance
106 : : // but estimate statistics and (potentially) PDFs (at the interval given by
107 : : // the user).
108 [ + + ]: 5509065 : if (it > 0)
109 [ + + ]: 11029718 : for (const auto& e : g_diffeqs)
110 : 5521609 : e.advance( m_particles, CkMyPe(), dt, t, moments );
111 : :
112 : : // Save time stepping data
113 : 5509065 : m_dt = dt;
114 : 5509065 : m_t = t;
115 : 5509065 : m_it = it;
116 : :
117 : : // Contribute number of particles we hit the particles output frequency
118 : : auto poseq =
119 : : !g_inputdeck.get< tag::param, tag::position, tag::depvar >().empty();
120 : : const auto parfreq =
121 : 5509065 : g_inputdeck.get< tag::output, tag::iter, tag::particles >();
122 : :
123 [ + - ]: 11018130 : CkCallback c( CkIndex_Integrator::out(), thisProxy[thisIndex] );
124 : :
125 [ + + ][ - + ]: 5509065 : if (poseq && !((m_it+1) % parfreq))
126 [ - - ][ - - ]: 0 : m_particlewriter[ CkMyNode() ].npar( m_particles.nunk(), c );
[ - - ][ - - ]
127 : : else
128 [ + - ]: 5509065 : c.send();
129 : 5509065 : }
130 : :
131 : : void
132 : 5509065 : Integrator::out()
133 : : // *****************************************************************************
134 : : // Output particle positions to file
135 : : // *****************************************************************************
136 : : {
137 : : auto poseq =
138 : : !g_inputdeck.get< tag::param, tag::position, tag::depvar >().empty();
139 : : const auto parfreq =
140 : 5509065 : g_inputdeck.get< tag::output, tag::iter, tag::particles >();
141 : :
142 [ + - ]: 11018130 : CkCallback c( CkIndex_Integrator::accumulate(), thisProxy[thisIndex] );
143 : :
144 : : // Output particles data to file if we hit the particles output frequency
145 [ + + ][ - + ]: 5509065 : if (poseq && !((m_it+1) % parfreq)) {
146 : : // query position eq offset in particle array (0: only first particle pos)
147 : 0 : auto po = g_inputdeck.get< tag::component >().offset< tag::position >( 0 );
148 : : // output particle positions to file
149 [ - - ]: 0 : m_particlewriter[ CkMyNode() ].
150 [ - - ][ - - ]: 0 : writeCoords( m_itp++, m_particles.extract(0,po),
[ - - ][ - - ]
[ - - ][ - - ]
151 [ - - ][ - - ]: 0 : m_particles.extract(1,po), m_particles.extract(2,po), c );
[ - - ][ - - ]
[ - - ][ - - ]
152 : : } else {
153 [ + - ]: 5509065 : c.send();
154 : : }
155 : 5509065 : }
156 : :
157 : : void
158 [ + + ]: 5509065 : Integrator::accumulate()
159 : : // *****************************************************************************
160 : : // Start collecting statistics
161 : : // *****************************************************************************
162 : : {
163 : : if (!g_inputdeck.stat()) {// if no stats to estimate, skip to end of time step
164 [ - - ]: 0 : contribute( CkCallback(CkReductionTarget(Distributor, nostat), m_host) );
165 : : } else {
166 : : // Accumulate sums for ordinary moments (every time step)
167 : 5509065 : accumulateOrd( m_it, m_t, m_dt );
168 : : }
169 : 5509065 : }
170 : :
171 : : void
172 : 5509065 : Integrator::accumulateOrd( uint64_t it, tk::real t, tk::real dt )
173 : : // *****************************************************************************
174 : : // Accumulate sums for ordinary moments and ordinary PDFs
175 : : //! \param[in] it Iteration count
176 : : //! \param[in] t Physical time
177 : : //! \param[in] dt Time step size
178 : : // *****************************************************************************
179 : : {
180 : 5509065 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
181 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
182 : 5509065 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
183 : 5509065 : const auto pdffreq = g_inputdeck.get< tag::output, tag::iter, tag::pdf >();
184 : :
185 : : // Accumulate partial sums for ordinary moments
186 : 5509065 : m_stat.accumulateOrd();
187 : : // Accumulate sums for ordinary PDFs at first and last iterations and at
188 : : // select times
189 [ + + ][ + + ]: 5509065 : if ( g_inputdeck.pdf() &&
190 : 253572 : ( it == 0 ||
191 [ + + ]: 253572 : !((it+1) % pdffreq) ||
192 [ - + ][ - - ]: 253532 : (std::fabs(t+dt-term) < eps && (it+1) >= nstep) ) )
193 : 104 : m_stat.accumulateOrdPDF();
194 : :
195 : : // Send accumulated ordinary moments and ordinary PDFs to collector for
196 : : // estimation
197 : 5509065 : m_coll.ckLocalBranch()->chareOrd( m_stat.ord(),
198 : : m_stat.oupdf(),
199 : : m_stat.obpdf(),
200 : : m_stat.otpdf() );
201 : 5509065 : }
202 : :
203 : : void
204 : 5509065 : Integrator::accumulateCen( uint64_t it,
205 : : tk::real t,
206 : : tk::real dt,
207 : : const std::vector< tk::real >& ord )
208 : : // *****************************************************************************
209 : : // Accumulate sums for central moments and central PDFs
210 : : //! \param[in] it Iteration count
211 : : //! \param[in] t Physical time
212 : : //! \param[in] dt Time step size
213 : : //! \param[in] ord Estimated ordinary moments (collected from all PEs)
214 : : // *****************************************************************************
215 : : {
216 : 5509065 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
217 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
218 : 5509065 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
219 : 5509065 : const auto pdffreq = g_inputdeck.get< tag::output, tag::iter, tag::pdf >();
220 : :
221 : : // Accumulate partial sums for central moments
222 : 5509065 : m_stat.accumulateCen( ord );
223 : : // Accumulate partial sums for central PDFs at first and last iteraions and
224 : : // at select times
225 [ + + ][ + + ]: 5509065 : if ( g_inputdeck.pdf() &&
226 : 253572 : ( it == 0 ||
227 [ + + ]: 253572 : !((it+1) % pdffreq) ||
228 [ - + ][ - - ]: 253532 : (std::fabs(t+dt-term) < eps && (it+1) >= nstep) ) )
229 : 104 : m_stat.accumulateCenPDF( ord );
230 : :
231 : : // Send accumulated central moments to host for estimation
232 : 5509065 : m_coll.ckLocalBranch()->chareCen( m_stat.ctr(),
233 : : m_stat.cupdf(),
234 : : m_stat.cbpdf(),
235 : : m_stat.ctpdf() );
236 : 5509065 : }
237 : :
238 : : #include "NoWarning/integrator.def.h"
|
