From 3c5fdc3fea947f4d9e89edc947f16ccde6f633cb Mon Sep 17 00:00:00 2001 From: Stephen Berry Date: Thu, 21 Nov 2019 08:10:41 -0600 Subject: [PATCH] More modular integrators and added to ChaiEngine New integrators: - Heun - RK3 - NCRK4 Added threading tools to ChaiEngine Add const qualifier to Sampler's operator()() --- include/ascent/ChaiEngine.h | 63 ++++++---- include/ascent/integrators_modular/Heun.h | 114 +++++++++++++++++ include/ascent/integrators_modular/NCRK4.h | 137 +++++++++++++++++++++ include/ascent/integrators_modular/RK3.h | 129 +++++++++++++++++++ include/ascent/threading/Queue.h | 1 + include/ascent/timing/Sampler.h | 2 +- 6 files changed, 418 insertions(+), 28 deletions(-) create mode 100644 include/ascent/integrators_modular/Heun.h create mode 100644 include/ascent/integrators_modular/NCRK4.h create mode 100644 include/ascent/integrators_modular/RK3.h diff --git a/include/ascent/ChaiEngine.h b/include/ascent/ChaiEngine.h index 62299be..700e7b6 100644 --- a/include/ascent/ChaiEngine.h +++ b/include/ascent/ChaiEngine.h @@ -15,6 +15,7 @@ #pragma once #include "ascent/Ascent.h" +#include "ascent/threading/Pool.h" #include "chaiscript/chaiscript.hpp" #include "chaiscript/chaiscript_stdlib.hpp" @@ -57,49 +58,57 @@ namespace asc ChaiEngine() { - add(chaiscript::vector_conversion()); - add(chaiscript::vector_conversion>()); - add(chaiscript::bootstrap::standard_library::vector_type("state_t")); + using namespace chaiscript; + + add(vector_conversion()); + add(vector_conversion>()); + add(bootstrap::standard_library::vector_type("state_t")); scriptRecorder(*this, "Recorder"); - add(chaiscript::fun([](RecorderT& rec, int& x) { rec.record(x); }), "record"); - add(chaiscript::fun([](RecorderT& rec, size_t& x) { rec.record(x); }), "record"); - add(chaiscript::fun([](RecorderT& rec, int& x, const std::string& title) { rec.record(x, title); }), "record"); - add(chaiscript::fun([](RecorderT& rec, size_t& x, const std::string& title) { rec.record(x, title); }), "record"); + add(fun([](RecorderT& rec, int& x) { rec.record(x); }), "record"); + add(fun([](RecorderT& rec, size_t& x) { rec.record(x); }), "record"); + add(fun([](RecorderT& rec, int& x, const std::string& title) { rec.record(x, title); }), "record"); + add(fun([](RecorderT& rec, size_t& x, const std::string& title) { rec.record(x, title); }), "record"); // This allows for more generic mixing of various data types, all they need is to be converted to a string prior to being passed into the recorder // The RecorderString is primarily for data that is going to be output to a file and not operated on during the simulation scriptRecorder(*this, "RecorderString"); - add(chaiscript::user_type(), "Param"); - add(chaiscript::constructor(), "Param"); - add(chaiscript::constructor(), "Param"); - add(chaiscript::fun([](asc::Param& param, const value_t y) { return param = y; }), "="); - add(chaiscript::fun([](asc::Param& param) { std::cout << param << '\n'; }), "print"); - add(chaiscript::fun([](asc::Param& param) { return static_cast(param); }), "value"); + add(user_type(), "Param"); + add(constructor(), "Param"); + add(constructor(), "Param"); + add(fun([](asc::Param& param, const value_t y) { return param = y; }), "="); + add(fun([](asc::Param& param) { std::cout << param << '\n'; }), "print"); + add(fun([](asc::Param& param) { return static_cast(param); }), "value"); - add(chaiscript::user_type(), "ParamV"); - add(chaiscript::constructor(), "ParamV"); - add(chaiscript::fun([](ParamV& lhs, const std::vector& rhs) { return lhs = rhs; }), "="); - add(chaiscript::fun([](ParamV& v) { v.zero(); }), "zero"); + add(user_type(), "ParamV"); + add(constructor(), "ParamV"); + add(fun([](ParamV& lhs, const std::vector& rhs) { return lhs = rhs; }), "="); + add(fun([](ParamV& v) { v.zero(); }), "zero"); - add(chaiscript::user_type(), "Sampler"); - add(chaiscript::constructor(), "Sampler"); - add(chaiscript::fun(&Sampler::operator()), "eval"); - add(chaiscript::fun(&Sampler::event), "event"); - add(chaiscript::fun(&Sampler::reset), "reset"); + add(user_type(), "Sampler"); + add(constructor(), "Sampler"); + add(fun(&Sampler::operator()), "eval"); + add(fun(&Sampler::event), "event"); + add(fun(&Sampler::reset), "reset"); // System - add(chaiscript::user_type(), "System"); - add(chaiscript::constructor(), "System"); - add(chaiscript::fun([](System& sys, const asc::state_t& x, asc::state_t& D, const asc::value_t t) { sys(x, D, t); }), "eval"); - add(chaiscript::type_conversion()); - add(chaiscript::fun([](System& sys, const system_t& func) { sys.push_back(func); }), "push_back"); + add(user_type(), "System"); + add(constructor(), "System"); + add(fun([](System& sys, const asc::state_t& x, asc::state_t& D, const asc::value_t t) { sys(x, D, t); }), "eval"); + add(type_conversion()); + add(fun([](System& sys, const system_t& func) { sys.push_back(func); }), "push_back"); // Integrators integrator("ascEuler"); integrator("RK2"); integrator("RK4"); + + // threading + Queue::script(*this, "Queue"); + Pool::script(*this, "Pool"); + add(fun([] { return std::thread::hardware_concurrency(); }), "hardware_concurrency"); + add(fun([](asc::Recorder& rec, const int sig_digits) { rec.precision = sig_digits; }), "precision"); } ChaiEngine(const ChaiEngine&) = default; ChaiEngine(ChaiEngine&&) = default; diff --git a/include/ascent/integrators_modular/Heun.h b/include/ascent/integrators_modular/Heun.h new file mode 100644 index 0000000..5d29876 --- /dev/null +++ b/include/ascent/integrators_modular/Heun.h @@ -0,0 +1,114 @@ +// Copyright (c) 2016-2019 Anyar, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "ascent/modular/Module.h" +#include "ascent/integrators_modular/ModularIntegrators.h" + +// Heun's Method + +namespace asc +{ + namespace modular + { + template + struct Heunprop : public Propagator + { + void operator()(State& state, const value_t dt) override + { + auto& x = *state.x; + auto& xd = *state.xd; + state.memory.resize(2); + auto& x0 = state.memory[0]; + auto& xd0 = state.memory[1]; + + switch (Propagator::pass) + { + case 0: + x0 = x; + xd0 = xd; //k1 + x = x0 + dt * xd; + break; + case 1: + x = x0 + dt * 0.5 * ( xd0 + xd ); + break; + } + } + }; + + template + struct Heunstepper : public TimeStepper + { + value_t t0{}; + + void operator()(const size_t pass, value_t& t, const value_t dt) override + { + switch (pass) + { + case 0: + t0 = t; + t += dt; + break; + case 1: + t = t0 + dt; + break; + default: + break; + } + } + }; + + template + struct Heun + { + Heunprop propagator; + Heunstepper stepper; + + template + void operator()(modules_t& modules, value_t& t, const value_t dt) + { + auto& pass = propagator.pass; + pass = 0; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + ++pass; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + } + + template + void update(modules_t& modules) + { + for (auto& module : modules) + { + (*module)(); + } + } + + template + void propagate(modules_t& modules, const value_t dt) + { + for (auto& module : modules) + { + module->propagate(propagator, dt); + } + } + }; + } +} \ No newline at end of file diff --git a/include/ascent/integrators_modular/NCRK4.h b/include/ascent/integrators_modular/NCRK4.h new file mode 100644 index 0000000..a177328 --- /dev/null +++ b/include/ascent/integrators_modular/NCRK4.h @@ -0,0 +1,137 @@ +// Copyright (c) 2016-2019 Anyar, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "ascent/modular/Module.h" +#include "ascent/integrators_modular/ModularIntegrators.h" + +// Non-confluent RK4 method ( "3/8th's Rule" ) + +namespace asc +{ + namespace modular + { + template + struct NCRK4prop : public Propagator + { + void operator()(State& state, const value_t dt) override + { + auto& x = *state.x; + auto& xd = *state.xd; + state.memory.resize(5); + auto& x0 = state.memory[0]; + auto& xd0 = state.memory[1]; + auto& xd1 = state.memory[2]; + auto& xd2 = state.memory[3]; + auto& xd3 = state.memory[4]; + + switch (Propagator::pass) + { + case 0: + x0 = x; + xd0 = xd; // k1 + x = x0 + (1.0/3.0) * dt * xd0; + break; + case 1: + xd1 = xd; // k2 + x = x0 + dt * ( (-1.0/3.0) * xd0 + xd1) ; + break; + case 2: + xd2 = xd; // k3 + x = x0 + dt * ( xd0 - xd1 + xd2); + break; + case 3: + xd3 = xd; // k4 + x = x0 + dt / 8.0 * (xd0 + 3.0 * xd1 + 3.0 * xd2 + xd3); + break; + } + } + }; + + template + struct NCRK4stepper : public TimeStepper + { + value_t t0{}; + + void operator()(const size_t pass, value_t& t, const value_t dt) override + { + switch (pass) + { + case 0: + t0 = t; + t += 1.0/3.0 * dt; + break; + case 1: + t = t0 + 2.0/3.0 * dt; + break; + case 2: + t = t0 + dt; + break; + default: + break; + } + } + }; + + template + struct NCRK4 + { + NCRK4prop propagator; + NCRK4stepper stepper; + + template + void operator()(modules_t& modules, value_t& t, const value_t dt) + { + auto& pass = propagator.pass; + pass = 0; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + ++pass; + + update(modules); + propagate(modules, dt); + ++pass; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + ++pass; + + update(modules); + propagate(modules, dt); + } + + template + void update(modules_t& modules) + { + for (auto& module : modules) + { + (*module)(); + } + } + + template + void propagate(modules_t& modules, const value_t dt) + { + for (auto& module : modules) + { + module->propagate(propagator, dt); + } + } + }; + } +} diff --git a/include/ascent/integrators_modular/RK3.h b/include/ascent/integrators_modular/RK3.h new file mode 100644 index 0000000..456aa0c --- /dev/null +++ b/include/ascent/integrators_modular/RK3.h @@ -0,0 +1,129 @@ +// Copyright (c) 2016-2019 Anyar, Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include "ascent/modular/Module.h" +#include "ascent/integrators_modular/ModularIntegrators.h" + +// Third order, Three stage Runge Kutta (Heun's 3rd Order). + +namespace asc +{ + namespace modular + { + template + struct RK3prop : public Propagator + { + void operator()(State& state, const value_t dt) override + { + auto& x = *state.x; + auto& xd = *state.xd; + state.memory.resize(4); + auto& x0 = state.memory[0]; + auto& xd0 = state.memory[1]; + auto& xd1 = state.memory[2]; + auto& xd2 = state.memory[3]; + + switch (Propagator::pass) + { + case 0: + x0 = x; + xd0 = xd; // k1 + x = x0 + dt * (1.0/3.0) * xd0; + break; + case 1: + xd1 = xd; // k2 + x = x0 + dt * (2.0/3.0) * xd1; + break; + case 2: + xd2 = xd; // k3 + x = x0 + dt * ( (1.0/4.0)*xd0 + (3.0/4.0)*xd2 ); + break; + } + } + }; + + template + struct RK3stepper : public TimeStepper + { + value_t t0{}; + + void operator()(const size_t pass, value_t& t, const value_t dt) override + { + switch (pass) + { + case 0: + t0 = t; + t += (1.0/3.0) * dt; + break; + case 1: + t = t0 + (2.0/3.0) * dt; + break; + case 2: + t = t0 + dt; + break; + default: + break; + } + } + }; + + template + struct RK3 + { + RK3prop propagator; + RK3stepper stepper; + + template + void operator()(modules_t& modules, value_t& t, const value_t dt) + { + auto& pass = propagator.pass; + pass = 0; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + ++pass; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + ++pass; + + update(modules); + propagate(modules, dt); + stepper(pass, t, dt); + } + + template + void update(modules_t& modules) + { + for (auto& module : modules) + { + (*module)(); + } + } + + template + void propagate(modules_t& modules, const value_t dt) + { + for (auto& module : modules) + { + module->propagate(propagator, dt); + } + } + }; + } +} diff --git a/include/ascent/threading/Queue.h b/include/ascent/threading/Queue.h index d14f3fa..26ab129 100644 --- a/include/ascent/threading/Queue.h +++ b/include/ascent/threading/Queue.h @@ -23,6 +23,7 @@ #include #include #include +#include // A threaded worker queue diff --git a/include/ascent/timing/Sampler.h b/include/ascent/timing/Sampler.h index 57de22e..3ae99af 100644 --- a/include/ascent/timing/Sampler.h +++ b/include/ascent/timing/Sampler.h @@ -37,7 +37,7 @@ namespace asc SamplerT& operator=(SamplerT&&) = default; ~SamplerT() noexcept { dt = dt_base; } - bool operator()(const T sample_rate) noexcept + bool operator()(const T sample_rate) const noexcept { const size_t n = static_cast((t + eps) / sample_rate); // the number of sample time steps that have occured const T sample_time = (n + 1) * sample_rate; // the next sample time