diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt index a7fba32..e7bd2a3 100644 --- a/src/physics/CMakeLists.txt +++ b/src/physics/CMakeLists.txt @@ -1,18 +1,19 @@ set(PHYSICS_INCLUDE_LIST constants.hpp laws.hpp units.hpp + maths.hpp ) add_custom_target(physics_inc SOURCES ${PHYSICS_INCLUDE_LIST} ) install(FILES ${PHYSICS_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/physics ) install(FILES io/units.hpp DESTINATION ${INCLUDE_INSTALL_DIR}/physics/io ) diff --git a/src/physics/maths.hpp b/src/physics/maths.hpp new file mode 100644 index 0000000..866c78e --- /dev/null +++ b/src/physics/maths.hpp @@ -0,0 +1,114 @@ +/*------------------------------------------------------------------------------ + +Copyright (c)2015 F. Georget , Princeton University +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this +list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation and/or +other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its contributors +may be used to endorse or promote products derived from this software without +specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +------------------------------------------------------------------------------*/ + +#ifndef SPECMICP_UTILS_MATHS_HPP +#define SPECMICP_UTILS_MATHS_HPP + +/*! +\file maths.hpp +\brief Simple math functions + */ + +#include "../types.hpp" + +namespace specmicp { + +//! \brief Average methods available +enum class Average { + arithmetic, + geometric, + harmonic +}; + +//! \brief Average two numbers +//! +//! \tparam average the average method +template +scalar_t SPECMICP_CONST_F average(scalar_t a, scalar_t b); + +template<> +inline +scalar_t SPECMICP_CONST_F average(scalar_t a, scalar_t b) +{ + return (a + b)/2.0; +} + +template<> +inline +scalar_t SPECMICP_CONST_F average(scalar_t a, scalar_t b) +{ + return 2.0/(1.0/a + 1.0/b); +} + +template<> +inline +scalar_t SPECMICP_CONST_F average(scalar_t a, scalar_t b) +{ + return std::sqrt(a*b); +} + + +//! \brief Average a vector +//! +//! \tparam average the average method +//! \param vector the vector to average +template +scalar_t SPECMICP_CONST_F average(const Vector& vector); + +template<> +inline +scalar_t SPECMICP_CONST_F average(const Vector& vector) +{ + return vector.sum()/vector.rows(); +} + +template<> +inline +scalar_t SPECMICP_CONST_F average(const Vector& vector) +{ + return std::pow(vector.prod(), 1.0/vector.rows()); +} + +template<> +inline +scalar_t SPECMICP_CONST_F average(const Vector& vector) +{ + scalar_t sum = 0; + for (index_t ind=0; ind using namespace specmicp; // TOC // ======================= // // 1 - Logger // 2 - Timer -// 3 - Moving Average -// 4 - Performance handler -// 5 - Options handler +// 3 - Average +// 4 - Moving Average +// 5 - Performance handler +// 6 - Options handler + +// Logger +// ======================= TEST_CASE("Logger", "[io],[log]") { SECTION("Non init logging") { init_logger(nullptr, logger::LogLevel::Warning); // compilation test, and runtime test of not failling nullptr SPAM << "Spam"; DEBUG << "Debug"; INFO << "Info"; WARNING << "Warning"; ERROR << "Error"; CRITICAL << "Critical"; } SECTION("Level test") { std::stringstream stream_log; init_logger(&stream_log, logger::LogLevel::Warning); SPAM << "Spam"; REQUIRE(stream_log.str().size() == 0); DEBUG << "Debug"; REQUIRE(stream_log.str().size() == 0); INFO << "Info"; REQUIRE(stream_log.str().size() == 0); WARNING << "Warning"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); ERROR << "Error"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); CRITICAL << "Critical"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); init_logger(&stream_log, logger::LogLevel::Debug); SPAM << "Spam"; REQUIRE(stream_log.str().size() == 0); DEBUG << "Debug"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); INFO << "Info"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); WARNING << "Warning"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); ERROR << "Error"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); CRITICAL << "Critical"; REQUIRE(stream_log.str().size() != 0); stream_log.str(""); } } -// Logger -// ======================= - // Timer // ======================= TEST_CASE("Timer", "[time],[CPU]") { SECTION("Timer test") { Timer timer; timer.stop(); CHECK(timer.elapsed_time() >= 0.0); timer.start(); timer.stop(); CHECK(timer.elapsed_time() >= 0.0); CHECK(timer.get_stop() >= timer.get_start()); CHECK(timer.get_ctime_start() != nullptr); CHECK(timer.get_ctime_stop() != nullptr); } } +// Average +// ==================== + +#define test_average(method, a, b, sol) \ + CHECK(average(a, b) == Approx(sol).epsilon(1e-10)); + +#define test_average_vector(method, vector, sol) \ + CHECK(average(vector) == Approx(sol).epsilon(1e-10)); + + +TEST_CASE("Average", "[average],[maths]") { + + Vector test_vector(4); + test_vector << 1.0, 2.0, 3.0, 4.0; + + SECTION("Arithmetic average") { + + test_average(Average::arithmetic, 1.0, 1.0, 1.0); + test_average(Average::arithmetic, 2.0, 1.0, 1.5); + test_average(Average::arithmetic, -1.0, 1.0, 0.0); + + test_average_vector(Average::arithmetic, test_vector, 10.0/4); + } + + SECTION("Harmonic average") { + + test_average(Average::harmonic, 1.0, 1.0, 1.0); + test_average(Average::harmonic, 2.0, 1.0, 2.0/(1+0.5)); + + test_average_vector(Average::harmonic, test_vector, 1.92); + } + + SECTION("Geometric average") { + + test_average(Average::geometric, 1.0, 1.0, 1.0); + test_average(Average::geometric, 1.0, 2.0, std::sqrt(2.0)); + + test_average_vector(Average::geometric, test_vector, std::pow(24.0, 1.0/4.0)); + + } +} + +#undef test_average +#undef test_average_vector + // Moving Average // ======================= TEST_CASE("Moving average", "[average],[timestep]") { SECTION("Moving average test") { utils::ExponentialMovingAverage moving_average(0.1, 1.0); REQUIRE(moving_average.current_value() == 1.0); CHECK(moving_average.add_point(1.0) == 1.0); CHECK(moving_average.add_point(2.0) == 1.1); REQUIRE(moving_average.add_point(3.0) == 0.9*1.1+0.3); moving_average.reset(1.0); REQUIRE(moving_average.current_value() == 1.0); moving_average.set_alpha(0.2); REQUIRE(moving_average.add_point(2.0) == Approx(0.8+0.4)); } } // Performance handler // ======================= struct MockPerf { scalar_t residuals {-1}; index_t nb_iterations {-1}; }; class MockSolverPerf: public PerformanceHandler { public: MockSolverPerf() {} void do_stuff() { get_perfs().residuals = 1e-6; get_perfs().nb_iterations = 10; } void reset_solver() { reset_perfs(); } }; TEST_CASE("PerformanceHandler", "[performance],[base]") { SECTION("Performance handler test") { auto my_solver = MockSolverPerf(); const auto& perfs = my_solver.get_perfs(); CHECK(perfs.nb_iterations == -1); CHECK(perfs.residuals == -1); my_solver.do_stuff(); CHECK(perfs.nb_iterations == 10); CHECK(perfs.residuals == 1e-6); my_solver.reset_solver(); CHECK(perfs.nb_iterations == -1); CHECK(perfs.residuals == -1); } } // Options handler // ======================= struct MockOptions { MockOptions() {} MockOptions(scalar_t res, index_t max_iter): residuals(res), max_iterations(max_iter) {} scalar_t residuals {1e-6}; index_t max_iterations {100}; }; class MockSolverOptions: public OptionsHandler { public: MockSolverOptions() {} MockSolverOptions(scalar_t res, index_t max_iter): OptionsHandler(res, max_iter) {} }; TEST_CASE("Options handler", "[options],[base]") { SECTION("Options handler test") { auto my_solver_default = MockSolverOptions(); const auto& ro_options = my_solver_default.get_options(); CHECK(ro_options.max_iterations == 100); CHECK(ro_options.residuals == 1e-6); auto& rw_options = my_solver_default.get_options(); rw_options.max_iterations = 10; rw_options.residuals = 1e-4; CHECK(ro_options.max_iterations == 10); CHECK(ro_options.residuals == 1e-4); } SECTION("Options handler - non default initialization") { auto my_solver = MockSolverOptions(1e-4, 10); const auto& ro_options = my_solver.get_options(); CHECK(ro_options.max_iterations == 10); CHECK(ro_options.residuals == 1e-4); } }