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