1 // Copyright (c) 2016-2019 Anyar, Inc.
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 //      http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14 
15 #pragma once
16 
17 #include "ascent/modular/Module.h"
18 #include "ascent/integrators_modular/ModularIntegrators.h"
19 
20 namespace asc
21 {
22    namespace modular
23    {
24       template <class value_t>
25       struct RK4prop : public Propagator<value_t>
26       {
operatorRK4prop27          void operator()(State& state, const value_t dt) override
28          {
29             auto& x = *state.x;
30             auto& xd = *state.xd;
31             if (state.memory.size() < 5)
32             {
33                state.memory.resize(5);
34             }
35             auto& x0 = state.memory[0];
36             auto& xd0 = state.memory[1];
37             auto& xd1 = state.memory[2];
38             auto& xd2 = state.memory[3];
39             auto& xd3 = state.memory[4];
40 
41             switch (Propagator<value_t>::pass)
42             {
43             case 0:
44                x0 = x;
45                xd0 = xd;
46                x = x0 + 0.5 * dt * xd0;
47                break;
48             case 1:
49                xd1 = xd;
50                x = x0 + 0.5 * dt * xd1;
51                break;
52             case 2:
53                xd2 = xd;
54                x = x0 + dt * xd2;
55                break;
56             case 3:
57                xd3 = xd;
58                x = x0 + dt / 6.0 * (xd0 + 2 * xd1 + 2 * xd2 + xd3);
59                break;
60             }
61          }
62       };
63 
64       template <class value_t>
65       struct RK4stepper : public TimeStepper<value_t>
66       {
67          value_t t0{};
68 
operatorRK4stepper69          void operator()(const size_t pass, value_t& t, const value_t dt) override
70          {
71             switch (pass)
72             {
73             case 0:
74                t0 = t;
75                t += 0.5 * dt;
76                break;
77             case 2:
78                t = t0 + dt;
79                break;
80             default:
81                break;
82             }
83          }
84       };
85 
86       template <class value_t>
87       struct RK4
88       {
89          static constexpr size_t n_substeps = 4;
90 
91          asc::Module* run_first{};
92 
93          RK4prop<value_t> propagator;
94          RK4stepper<value_t> stepper;
95 
96          template <class modules_t>
operatorRK497          void operator()(modules_t& blocks, value_t& t, const value_t dt)
98          {
99             auto& pass = propagator.pass;
100             pass = 0;
101 
102             update(blocks, run_first);
103             apply(blocks);
104             propagate(blocks, propagator, dt);
105             stepper(pass, t, dt);
106             postprop(blocks);
107             ++pass;
108 
109             update(blocks, run_first);
110             apply(blocks);
111             propagate(blocks, propagator, dt);
112             postprop(blocks);
113             ++pass;
114 
115             update(blocks, run_first);
116             apply(blocks);
117             propagate(blocks, propagator, dt);
118             stepper(pass, t, dt);
119             postprop(blocks);
120             ++pass;
121 
122             update(blocks, run_first);
123             apply(blocks);
124             propagate(blocks, propagator, dt);
125             postprop(blocks);
126          }
127       };
128    }
129 }
130