Line data Source code
1 : // Copyright (C) 2023 Helmut Schmidt University
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder
4 :
5 : #pragma once
6 :
7 : namespace coupling {
8 : namespace services {
9 : template <unsigned int dim> class ParallelTimeIntegrationService;
10 : } // namespace services
11 : } // namespace coupling
12 :
13 : #include "coupling/configurations/MaMiCoConfiguration.h"
14 : #include "coupling/interface/PintableMacroSolver.h"
15 : #include "coupling/scenario/Scenario.h"
16 : #include <functional>
17 :
18 : // Convenience operators, to be able to write parareal iterations in a readable notation
19 : using st_ptr = std::unique_ptr<coupling::interface::PintableMacroSolverState>;
20 : inline st_ptr operator+(const st_ptr& lhs, const st_ptr& rhs) { return *lhs + *rhs; }
21 : inline st_ptr operator-(const st_ptr& lhs, const st_ptr& rhs) { return *lhs - *rhs; }
22 :
23 : /**
24 : * Service to manage timeloop of a coupled simulation scenario. Supports sequential or parallel-in-time integration using a Parareal variant,
25 : * as described in "Blumers, A. L., Li, Z., & Karniadakis, G. E. (2019). Supervised parallel-in-time algorithm for long-time Lagrangian
26 : * simulations of stochastic dynamics: Application to hydrodynamics. Journal of Computational Physics, 393, 214-228".
27 : *
28 : * @author Piet Jarmatz
29 : */
30 : template <unsigned int dim> class coupling::services::ParallelTimeIntegrationService {
31 : public:
32 : using State = coupling::interface::PintableMacroSolverState;
33 : using Solver = coupling::interface::PintableMacroSolver;
34 :
35 : ParallelTimeIntegrationService(coupling::configurations::MaMiCoConfiguration<dim> mamicoConfig, Scenario* scenario)
36 : : _pint_domain(0), _rank(0), _ranks_per_domain(1), _world_rank(0), _cfg(mamicoConfig.getTimeIntegrationConfiguration()), _scenario(scenario) {
37 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
38 : int world_size;
39 : MPI_Comm_size(MPI_COMM_WORLD, &world_size);
40 : MPI_Comm_rank(MPI_COMM_WORLD, &_world_rank);
41 : if (_cfg.isPintEnabled()) {
42 : if (world_size % _cfg.getPintDomains() != 0) {
43 : if (_world_rank == 0) {
44 : std::cout << "ERROR coupling::services::ParallelTimeIntegrationService: " << "MPI ranks not divisible by number of PinT subdomains!" << std::endl;
45 : std::cout << "When PinT is used, the number of required MPI ranks increases by a factor of number-subdomains." << std::endl;
46 : std::cout << "Check your configuration." << std::endl;
47 : }
48 : exit(EXIT_FAILURE);
49 : }
50 : _ranks_per_domain = world_size / _cfg.getPintDomains();
51 : _pint_domain = _world_rank / _ranks_per_domain;
52 : // This initializes _local_pint_comm by splitting MPI_COMM_WORLD into getPintDomains() disjoint communicators
53 : MPI_Comm_split(MPI_COMM_WORLD, _pint_domain, _world_rank, &_local_pint_comm);
54 : } else {
55 : _local_pint_comm = MPI_COMM_WORLD;
56 : }
57 : MPI_Comm_rank(_local_pint_comm, &_rank);
58 :
59 : #ifdef PINT_DEBUG
60 : std::cout << "PINT_DEBUG: world_rank " << _world_rank << " is rank " << _rank << " in pint domain " << _pint_domain << std::endl;
61 : #endif
62 :
63 : #endif
64 : }
65 :
66 : void run(int num_cycles) {
67 : if (!_cfg.isPintEnabled())
68 : run_cycles(0, num_cycles);
69 : else {
70 : PintDomain domain = setup_domain(num_cycles);
71 : setup_solvers(domain);
72 : init_parareal();
73 : run_parareal(_cfg.getPintIterations());
74 : }
75 : }
76 :
77 : int getPintDomain() const { return _pint_domain; }
78 : int getRank() const { return _rank; }
79 0 : bool isPintEnabled() const { return _cfg.isPintEnabled(); }
80 0 : int getIteration() const { return _iteration; }
81 :
82 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
83 : MPI_Comm getPintComm() const { return _local_pint_comm; }
84 : #endif
85 :
86 : private:
87 : void run_cycles(int start, int end) {
88 : for (int cycle = start; cycle < end; cycle++)
89 : _scenario->runOneCouplingCycle(cycle);
90 : }
91 :
92 : bool isFirst() const { return _pint_domain == 0; }
93 : bool isLast() const { return _pint_domain == _cfg.getPintDomains() - 1; }
94 :
95 : struct PintDomain {
96 : /** @brief number of time steps (coupling cycles) in this temporal domain */
97 : int size;
98 : /** @brief number of the first coupling cycle in this temporal domain (inclusive) */
99 : int minCycle;
100 : /** @brief number of the last coupling cycle of this temporal domain (exclusive) */
101 : int maxCycle;
102 : };
103 :
104 : PintDomain setup_domain(int num_cycles) const {
105 : PintDomain res;
106 : res.size = (int)(num_cycles / _cfg.getPintDomains());
107 : res.minCycle = _pint_domain * res.size;
108 : res.maxCycle = res.minCycle + res.size;
109 : // In case num_cycles is not divisible by _cfg.getPintDomains(), the last domain gets the remainder
110 : if (isLast())
111 : res.maxCycle = num_cycles;
112 : #ifdef PINT_DEBUG
113 : if (_rank == 0) {
114 : std::cout << "PINT_DEBUG: _pint_domain " << _pint_domain << " has minCycle " << res.minCycle << " and maxCycle " << res.maxCycle << std::endl;
115 : }
116 : #endif
117 : return res;
118 : }
119 :
120 : void setup_solvers(PintDomain domain) {
121 : auto solver = dynamic_cast<Solver*>(_scenario->getSolver());
122 : if (solver == nullptr) {
123 : std::cout << "ERROR coupling::services::ParallelTimeIntegrationService: "
124 : << "macroscopic solver is not pintable (= not compatible with parallel in time coupling)" << std::endl;
125 : exit(EXIT_FAILURE);
126 : }
127 : #ifdef PINT_DEBUG
128 : if (_cfg.getViscMultiplier() != 1.0)
129 : if (_world_rank == 0) {
130 : std::cout << "PINT_DEBUG: Starting supervisor with viscosity modified by " << _cfg.getViscMultiplier() << std::endl;
131 : }
132 : #endif
133 : _supervisor = solver->getSupervisor(domain.size, _cfg.getViscMultiplier());
134 : _F = [this, solver, domain](const std::unique_ptr<State>& s) {
135 : solver->setState(s, domain.minCycle);
136 : if (domain.minCycle > 0)
137 : _scenario->equilibrateMicro();
138 : run_cycles(domain.minCycle, domain.maxCycle);
139 : return solver->getState();
140 : };
141 : _G = [this, domain](const std::unique_ptr<State>& s) {
142 : auto& G = *_supervisor;
143 : return G(s, domain.minCycle);
144 : };
145 : _u_0 = solver->getState();
146 : }
147 :
148 : void receive(std::unique_ptr<State>& state) const {
149 : if (!state)
150 : state = _u_0->clone();
151 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
152 : int source_rank = _world_rank - _ranks_per_domain;
153 : MPI_Recv(state->getData(), state->getSizeBytes(), MPI_BYTE, source_rank, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
154 : #endif
155 : }
156 :
157 : void send(std::unique_ptr<State>& state) const {
158 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
159 : int destination_rank = _world_rank + _ranks_per_domain;
160 : MPI_Send(state->getData(), state->getSizeBytes(), MPI_BYTE, destination_rank, 0, MPI_COMM_WORLD);
161 : #endif
162 : }
163 :
164 : void get_past_state(std::unique_ptr<State>& state) const {
165 : if (isFirst())
166 : state = _u_0->clone();
167 : else
168 : receive(state);
169 : }
170 :
171 : void init_parareal() {
172 : get_past_state(_u_last_past);
173 : _u_last_future = _G(_u_last_past);
174 : if (!isLast())
175 : send(_u_last_future);
176 : }
177 :
178 : void run_parareal(int iterations) {
179 : while (_iteration < iterations) {
180 : // Correction step
181 : auto delta = _F(_u_last_past) - _G(_u_last_past);
182 : // Alternative variant, together with second _u_last_future line below.
183 : // Better scalability, should yield the same results. TODO test and verify.
184 : // auto delta = _F(_u_last_past) - _u_last_future;
185 :
186 : _iteration++;
187 :
188 : // Prediction step
189 : get_past_state(_u_next_past);
190 : auto prediction = _G(_u_next_past);
191 :
192 : // Refinement step
193 : _u_next_future = prediction + delta;
194 : if (!isLast())
195 : send(_u_next_future);
196 :
197 : // move for next iteration
198 : _u_last_past = std::move(_u_next_past);
199 : _u_last_future = std::move(_u_next_future);
200 : // Alternative variant
201 : //_u_last_future = std::move(prediction);
202 : }
203 :
204 : #ifdef PINT_DEBUG
205 : if (_world_rank == 0) {
206 : std::cout << "PINT_DEBUG: Finished all PinT iterations. " << std::endl;
207 : }
208 : #endif
209 : }
210 :
211 : int _pint_domain; // the index of the time domain to which this process belongs to
212 : int _rank; // rank of current process in _local_pint_comm
213 : int _ranks_per_domain; // number of MPI ranks in each time domain
214 : int _world_rank; // rank of current process in MPI_COMM_WORLD
215 : #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES)
216 : MPI_Comm _local_pint_comm; // the communicator of the local time domain of this rank
217 : #endif
218 : coupling::configurations::TimeIntegrationConfiguration _cfg;
219 : Scenario* _scenario;
220 : std::unique_ptr<Solver> _supervisor; // The supervisor, i.e. the coarse predictor
221 : std::function<std::unique_ptr<State>(const std::unique_ptr<State>&)> _F;
222 : std::function<std::unique_ptr<State>(const std::unique_ptr<State>&)> _G;
223 :
224 : // These objects represent coupled simulation states. There are used by the supervised parallel in time algorithm for operations
225 : // "last" and "next" describe two consecutive parareal iterations
226 : // "past" and "future" describe two points in simulation time, one pint time domain apart
227 : std::unique_ptr<State> _u_0; // initial state
228 : std::unique_ptr<State> _u_last_past;
229 : std::unique_ptr<State> _u_last_future;
230 : std::unique_ptr<State> _u_next_past;
231 : std::unique_ptr<State> _u_next_future;
232 :
233 : int _iteration{0};
234 : };
|
