commit e63239ed4710fada267c946a55642cbadbd245c5 Author: Jacob Date: Fri Oct 22 21:58:37 2021 +0200 Refactor diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..4ede542 --- /dev/null +++ b/.clang-format @@ -0,0 +1,98 @@ +AccessModifierOffset: -2 +AlignAfterOpenBracket: DontAlign +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlines: Left +AlignOperands: true +AlignTrailingComments: false +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: false +BinPackArguments: false +BinPackParameters: false +BraceWrapping: + AfterClass: true + AfterControlStatement: false + AfterEnum: false + AfterFunction: true + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: true + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false + SplitEmptyFunction: false + SplitEmptyNamespace: true + SplitEmptyRecord: true +BreakAfterJavaFieldAnnotations: true +BreakBeforeBinaryOperators: NonAssignment +BreakBeforeBraces: Custom +BreakBeforeInheritanceComma: true +BreakBeforeTernaryOperators: true +BreakConstructorInitializers: BeforeColon +BreakConstructorInitializersBeforeComma: false +BreakStringLiterals: true +ColumnLimit: 0 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 2 +ContinuationIndentWidth: 2 +Cpp11BracedListStyle: false +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: true +FixNamespaceComments: true +ForEachMacros: +- foreach +- Q_FOREACH +- BOOST_FOREACH +IncludeCategories: +- Priority: 2 + Regex: ^"(llvm|llvm-c|clang|clang-c)/ +- Priority: 3 + Regex: ^(<|"(gtest|gmock|isl|json)/) +- Priority: 1 + Regex: .* +IncludeIsMainRegex: (Test)?$ +IndentCaseLabels: false +IndentWidth: 2 +IndentWrappedFunctionNames: true +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +Language: Cpp +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 2 +NamespaceIndentation: Inner +ObjCBlockIndentWidth: 7 +ObjCSpaceAfterProperty: true +ObjCSpaceBeforeProtocolList: false +PointerAlignment: Right +ReflowComments: true +SortIncludes: false +SortUsingDeclarations: false +SpaceAfterCStyleCast: false +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 0 +SpacesInAngles: false +SpacesInCStyleCastParentheses: false +SpacesInContainerLiterals: true +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 8 +UseTab: Never + diff --git a/.clang-tidy b/.clang-tidy new file mode 100644 index 0000000..c5b7916 --- /dev/null +++ b/.clang-tidy @@ -0,0 +1,7 @@ +--- +Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers' +WarningsAsErrors: '' +HeaderFilterRegex: '' +FormatStyle: none + + diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d2475db --- /dev/null +++ b/.gitignore @@ -0,0 +1,28 @@ +# Build directories and binary files +build/ +out/ +cmake-build-*/ + +# IDE files +.vs/ +.idea/ +.vscode/ +!.vscode/settings.json +!.vscode/tasks.json +!.vscode/launch.json +!.vscode/extensions.json +*.swp +*~ + +# OS Generated Files +.DS_Store +.AppleDouble +.LSOverride +._* +.Spotlight-V100 +.Trashes +.Trash-* +$RECYCLE.BIN/ +.TemporaryItems +ehthumbs.db +Thumbs.db \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..b16f771 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,81 @@ +cmake_minimum_required(VERSION 3.15) + +# Set the project name to your project name, my project isn't very descriptive +project(myproject CXX) +include(cmake/StandardProjectSettings.cmake) +include(cmake/PreventInSourceBuilds.cmake) + +# Link this 'library' to set the c++ standard / compile-time options requested +add_library(project_options INTERFACE) +target_compile_features(project_options INTERFACE cxx_std_20) + +if(CMAKE_CXX_COMPILER_ID MATCHES ".*Clang") + option(ENABLE_BUILD_WITH_TIME_TRACE "Enable -ftime-trace to generate time tracing .json files on clang" OFF) + if(ENABLE_BUILD_WITH_TIME_TRACE) + target_compile_options(project_options INTERFACE -ftime-trace) + endif() +endif() + +# Link this 'library' to use the warnings specified in CompilerWarnings.cmake +add_library(project_warnings INTERFACE) + +# enable cache system +include(cmake/Cache.cmake) + +# Add linker configuration +include(cmake/Linker.cmake) +configure_linker(project_options) + +# standard compiler warnings +include(cmake/CompilerWarnings.cmake) +set_project_warnings(project_warnings) + +# sanitizer options if supported by compiler +include(cmake/Sanitizers.cmake) +enable_sanitizers(project_options) + +# enable doxygen +include(cmake/Doxygen.cmake) +enable_doxygen() + +# allow for static analysis options +include(cmake/StaticAnalyzers.cmake) + +option(BUILD_SHARED_LIBS "Enable compilation of shared libraries" OFF) +option(ENABLE_TESTING "Enable Test Builds" ON) +option(ENABLE_FUZZING "Enable Fuzzing Builds" OFF) + +# Very basic PCH example +option(ENABLE_PCH "Enable Precompiled Headers" OFF) +if(ENABLE_PCH) + # This sets a global PCH parameter, each project will build its own PCH, which is a good idea if any #define's change + # + # consider breaking this out per project as necessary + target_precompile_headers( + project_options + INTERFACE + + + + ) +endif() + +option(ENABLE_CONAN "Use Conan for dependency management" ON) +if(ENABLE_CONAN) + include(cmake/Conan.cmake) + run_conan() +endif() + +if(ENABLE_TESTING) + enable_testing() + message("Building Tests.") + add_subdirectory(test) +endif() + +if(ENABLE_FUZZING) + message("Building Fuzz Tests, using fuzzing sanitizer https://www.llvm.org/docs/LibFuzzer.html") + add_subdirectory(fuzz_test) +endif() + +add_subdirectory(src) + diff --git a/cmake/Cache.cmake b/cmake/Cache.cmake new file mode 100644 index 0000000..31f5e7e --- /dev/null +++ b/cmake/Cache.cmake @@ -0,0 +1,29 @@ +option(ENABLE_CACHE "Enable cache if available" ON) +if(NOT ENABLE_CACHE) + return() +endif() + +set(CACHE_OPTION + "ccache" + CACHE STRING "Compiler cache to be used") +set(CACHE_OPTION_VALUES "ccache" "sccache") +set_property(CACHE CACHE_OPTION PROPERTY STRINGS ${CACHE_OPTION_VALUES}) +list( + FIND + CACHE_OPTION_VALUES + ${CACHE_OPTION} + CACHE_OPTION_INDEX) + +if(${CACHE_OPTION_INDEX} EQUAL -1) + message( + STATUS + "Using custom compiler cache system: '${CACHE_OPTION}', explicitly supported entries are ${CACHE_OPTION_VALUES}") +endif() + +find_program(CACHE_BINARY ${CACHE_OPTION}) +if(CACHE_BINARY) + message(STATUS "${CACHE_OPTION} found and enabled") + set(CMAKE_CXX_COMPILER_LAUNCHER ${CACHE_BINARY}) +else() + message(WARNING "${CACHE_OPTION} is enabled but was not found. Not using it") +endif() diff --git a/cmake/CompilerWarnings.cmake b/cmake/CompilerWarnings.cmake new file mode 100644 index 0000000..c0531f1 --- /dev/null +++ b/cmake/CompilerWarnings.cmake @@ -0,0 +1,79 @@ +# from here: +# +# https://github.com/lefticus/cppbestpractices/blob/master/02-Use_the_Tools_Available.md + +function(set_project_warnings project_name) + option(WARNINGS_AS_ERRORS "Treat compiler warnings as errors" ON) + + set(MSVC_WARNINGS + /W4 # Baseline reasonable warnings + /w14242 # 'identifier': conversion from 'type1' to 'type1', possible loss of data + /w14254 # 'operator': conversion from 'type1:field_bits' to 'type2:field_bits', possible loss of data + /w14263 # 'function': member function does not override any base class virtual member function + /w14265 # 'classname': class has virtual functions, but destructor is not virtual instances of this class may not + # be destructed correctly + /w14287 # 'operator': unsigned/negative constant mismatch + /we4289 # nonstandard extension used: 'variable': loop control variable declared in the for-loop is used outside + # the for-loop scope + /w14296 # 'operator': expression is always 'boolean_value' + /w14311 # 'variable': pointer truncation from 'type1' to 'type2' + /w14545 # expression before comma evaluates to a function which is missing an argument list + /w14546 # function call before comma missing argument list + /w14547 # 'operator': operator before comma has no effect; expected operator with side-effect + /w14549 # 'operator': operator before comma has no effect; did you intend 'operator'? + /w14555 # expression has no effect; expected expression with side- effect + /w14619 # pragma warning: there is no warning number 'number' + /w14640 # Enable warning on thread un-safe static member initialization + /w14826 # Conversion from 'type1' to 'type_2' is sign-extended. This may cause unexpected runtime behavior. + /w14905 # wide string literal cast to 'LPSTR' + /w14906 # string literal cast to 'LPWSTR' + /w14928 # illegal copy-initialization; more than one user-defined conversion has been implicitly applied + /permissive- # standards conformance mode for MSVC compiler. + ) + + set(CLANG_WARNINGS + -Wall + -Wextra # reasonable and standard + -Wshadow # warn the user if a variable declaration shadows one from a parent context + -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. This helps + # catch hard to track down memory errors + -Wold-style-cast # warn for c-style casts + -Wcast-align # warn for potential performance problem casts + -Wunused # warn on anything being unused + -Woverloaded-virtual # warn if you overload (not override) a virtual function + -Wpedantic # warn if non-standard C++ is used + -Wconversion # warn on type conversions that may lose data + -Wsign-conversion # warn on sign conversions + -Wnull-dereference # warn if a null dereference is detected + -Wdouble-promotion # warn if float is implicit promoted to double + -Wformat=2 # warn on security issues around functions that format output (ie printf) + -Wimplicit-fallthrough # warn on statements that fallthrough without an explicit annotation + ) + + if(WARNINGS_AS_ERRORS) + set(CLANG_WARNINGS ${CLANG_WARNINGS} -Werror) + set(MSVC_WARNINGS ${MSVC_WARNINGS} /WX) + endif() + + set(GCC_WARNINGS + ${CLANG_WARNINGS} + -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist + -Wduplicated-cond # warn if if / else chain has duplicated conditions + -Wduplicated-branches # warn if if / else branches have duplicated code + -Wlogical-op # warn about logical operations being used where bitwise were probably wanted + -Wuseless-cast # warn if you perform a cast to the same type + ) + + if(MSVC) + set(PROJECT_WARNINGS ${MSVC_WARNINGS}) + elseif(CMAKE_CXX_COMPILER_ID MATCHES ".*Clang") + set(PROJECT_WARNINGS ${CLANG_WARNINGS}) + elseif(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + set(PROJECT_WARNINGS ${GCC_WARNINGS}) + else() + message(AUTHOR_WARNING "No compiler warnings set for '${CMAKE_CXX_COMPILER_ID}' compiler.") + endif() + + target_compile_options(${project_name} INTERFACE ${PROJECT_WARNINGS}) + +endfunction() diff --git a/cmake/Conan.cmake b/cmake/Conan.cmake new file mode 100644 index 0000000..fa4a0b7 --- /dev/null +++ b/cmake/Conan.cmake @@ -0,0 +1,58 @@ +macro(run_conan) + # Download automatically, you can also just copy the conan.cmake file + if(NOT EXISTS "${CMAKE_BINARY_DIR}/conan.cmake") + message( + STATUS + "Downloading conan.cmake from https://github.com/conan-io/cmake-conan") + file( + DOWNLOAD + "https://raw.githubusercontent.com/conan-io/cmake-conan/v0.16.1/conan.cmake" + "${CMAKE_BINARY_DIR}/conan.cmake" + EXPECTED_HASH + SHA256=396e16d0f5eabdc6a14afddbcfff62a54a7ee75c6da23f32f7a31bc85db23484 + TLS_VERIFY ON) + endif() + + set(ENV{CONAN_REVISIONS_ENABLED} 1) + list(APPEND CMAKE_MODULE_PATH ${CMAKE_BINARY_DIR}) + list(APPEND CMAKE_PREFIX_PATH ${CMAKE_BINARY_DIR}) + + include(${CMAKE_BINARY_DIR}/conan.cmake) + + # Add (or remove) remotes as needed + # conan_add_remote(NAME conan-center URL https://conan.bintray.com) + conan_add_remote(NAME cci URL https://center.conan.io INDEX 0) + conan_add_remote( + NAME bincrafters URL + https://bincrafters.jfrog.io/artifactory/api/conan/public-conan) + + # For multi configuration generators, like VS and XCode + if(NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Single configuration build!") + set(LIST_OF_BUILD_TYPES ${CMAKE_BUILD_TYPE}) + else() + message(STATUS "Multi-configuration build: '${CMAKE_CONFIGURATION_TYPES}'!") + set(LIST_OF_BUILD_TYPES ${CMAKE_CONFIGURATION_TYPES}) + endif() + + foreach(TYPE ${LIST_OF_BUILD_TYPES}) + message(STATUS "Running Conan for build type '${TYPE}'") + + # Detects current build settings to pass into conan + conan_cmake_autodetect(settings BUILD_TYPE ${TYPE}) + + # PATH_OR_REFERENCE ${CMAKE_SOURCE_DIR} is used to tell conan to process + # the external "conanfile.py" provided with the project + # Alternatively a conanfile.txt could be used + conan_cmake_install( + PATH_OR_REFERENCE + ${CMAKE_SOURCE_DIR} + BUILD + missing + # Pass compile-time configured options into conan + OPTIONS + SETTINGS + ${settings}) + endforeach() + +endmacro() diff --git a/cmake/Doxygen.cmake b/cmake/Doxygen.cmake new file mode 100644 index 0000000..4dad807 --- /dev/null +++ b/cmake/Doxygen.cmake @@ -0,0 +1,11 @@ +function(enable_doxygen) + option(ENABLE_DOXYGEN "Enable doxygen doc builds of source" OFF) + if(ENABLE_DOXYGEN) + set(DOXYGEN_CALLER_GRAPH YES) + set(DOXYGEN_CALL_GRAPH YES) + set(DOXYGEN_EXTRACT_ALL YES) + find_package(Doxygen REQUIRED dot) + doxygen_add_docs(doxygen-docs ${PROJECT_SOURCE_DIR}) + + endif() +endfunction() diff --git a/cmake/Linker.cmake b/cmake/Linker.cmake new file mode 100644 index 0000000..8b9471e --- /dev/null +++ b/cmake/Linker.cmake @@ -0,0 +1,33 @@ +option(ENABLE_USER_LINKER "Enable a specific linker if available" OFF) + +include(CheckCXXCompilerFlag) + +set(USER_LINKER_OPTION + "lld" + CACHE STRING "Linker to be used") +set(USER_LINKER_OPTION_VALUES "lld" "gold" "bfd") +set_property(CACHE USER_LINKER_OPTION PROPERTY STRINGS ${USER_LINKER_OPTION_VALUES}) +list( + FIND + USER_LINKER_OPTION_VALUES + ${USER_LINKER_OPTION} + USER_LINKER_OPTION_INDEX) + +if(${USER_LINKER_OPTION_INDEX} EQUAL -1) + message( + STATUS + "Using custom linker: '${USER_LINKER_OPTION}', explicitly supported entries are ${USER_LINKER_OPTION_VALUES}") +endif() + +function(configure_linker project_name) + if(NOT ENABLE_USER_LINKER) + return() + endif() + + set(LINKER_FLAG "-fuse-ld=${USER_LINKER_OPTION}") + + check_cxx_compiler_flag(${LINKER_FLAG} CXX_SUPPORTS_USER_LINKER) + if(CXX_SUPPORTS_USER_LINKER) + target_compile_options(${project_name} INTERFACE ${LINKER_FLAG}) + endif() +endfunction() diff --git a/cmake/PreventInSourceBuilds.cmake b/cmake/PreventInSourceBuilds.cmake new file mode 100644 index 0000000..57d9c59 --- /dev/null +++ b/cmake/PreventInSourceBuilds.cmake @@ -0,0 +1,18 @@ +# +# This function will prevent in-source builds +function(AssureOutOfSourceBuilds) + # make sure the user doesn't play dirty with symlinks + get_filename_component(srcdir "${CMAKE_SOURCE_DIR}" REALPATH) + get_filename_component(bindir "${CMAKE_BINARY_DIR}" REALPATH) + + # disallow in-source builds + if("${srcdir}" STREQUAL "${bindir}") + message("######################################################") + message("Warning: in-source builds are disabled") + message("Please create a separate build directory and run cmake from there") + message("######################################################") + message(FATAL_ERROR "Quitting configuration") + endif() +endfunction() + +assureoutofsourcebuilds() diff --git a/cmake/Sanitizers.cmake b/cmake/Sanitizers.cmake new file mode 100644 index 0000000..55cdf8c --- /dev/null +++ b/cmake/Sanitizers.cmake @@ -0,0 +1,67 @@ +function(enable_sanitizers project_name) + + if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU" OR CMAKE_CXX_COMPILER_ID MATCHES ".*Clang") + option(ENABLE_COVERAGE "Enable coverage reporting for gcc/clang" OFF) + + if(ENABLE_COVERAGE) + target_compile_options(${project_name} INTERFACE --coverage -O0 -g) + target_link_libraries(${project_name} INTERFACE --coverage) + endif() + + set(SANITIZERS "") + + option(ENABLE_SANITIZER_ADDRESS "Enable address sanitizer" OFF) + if(ENABLE_SANITIZER_ADDRESS) + list(APPEND SANITIZERS "address") + endif() + + option(ENABLE_SANITIZER_LEAK "Enable leak sanitizer" OFF) + if(ENABLE_SANITIZER_LEAK) + list(APPEND SANITIZERS "leak") + endif() + + option(ENABLE_SANITIZER_UNDEFINED_BEHAVIOR "Enable undefined behavior sanitizer" OFF) + if(ENABLE_SANITIZER_UNDEFINED_BEHAVIOR) + list(APPEND SANITIZERS "undefined") + endif() + + option(ENABLE_SANITIZER_THREAD "Enable thread sanitizer" OFF) + if(ENABLE_SANITIZER_THREAD) + if("address" IN_LIST SANITIZERS OR "leak" IN_LIST SANITIZERS) + message(WARNING "Thread sanitizer does not work with Address and Leak sanitizer enabled") + else() + list(APPEND SANITIZERS "thread") + endif() + endif() + + option(ENABLE_SANITIZER_MEMORY "Enable memory sanitizer" OFF) + if(ENABLE_SANITIZER_MEMORY AND CMAKE_CXX_COMPILER_ID MATCHES ".*Clang") + message(WARNING "Memory sanitizer requires all the code (including libc++) to be MSan-instrumented otherwise it reports false positives") + if("address" IN_LIST SANITIZERS + OR "thread" IN_LIST SANITIZERS + OR "leak" IN_LIST SANITIZERS) + message(WARNING "Memory sanitizer does not work with Address, Thread and Leak sanitizer enabled") + else() + list(APPEND SANITIZERS "memory") + endif() + endif() + + list( + JOIN + SANITIZERS + "," + LIST_OF_SANITIZERS) + + endif() + + if(LIST_OF_SANITIZERS) + if(NOT + "${LIST_OF_SANITIZERS}" + STREQUAL + "") + target_compile_options(${project_name} INTERFACE -fsanitize=${LIST_OF_SANITIZERS}) + target_link_options(${project_name} INTERFACE -fsanitize=${LIST_OF_SANITIZERS}) + endif() + endif() + +endfunction() diff --git a/cmake/StandardProjectSettings.cmake b/cmake/StandardProjectSettings.cmake new file mode 100644 index 0000000..59c36aa --- /dev/null +++ b/cmake/StandardProjectSettings.cmake @@ -0,0 +1,42 @@ +# Set a default build type if none was specified +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Setting build type to 'RelWithDebInfo' as none was specified.") + set(CMAKE_BUILD_TYPE + RelWithDebInfo + CACHE STRING "Choose the type of build." FORCE) + # Set the possible values of build type for cmake-gui, ccmake + set_property( + CACHE CMAKE_BUILD_TYPE + PROPERTY STRINGS + "Debug" + "Release" + "MinSizeRel" + "RelWithDebInfo") +endif() + +# Generate compile_commands.json to make it easier to work with clang based tools +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +option(ENABLE_IPO "Enable Interprocedural Optimization, aka Link Time Optimization (LTO)" OFF) + +if(ENABLE_IPO) + include(CheckIPOSupported) + check_ipo_supported( + RESULT + result + OUTPUT + output) + if(result) + set(CMAKE_INTERPROCEDURAL_OPTIMIZATION ON) + else() + message(SEND_ERROR "IPO is not supported: ${output}") + endif() +endif() +if(CMAKE_CXX_COMPILER_ID MATCHES ".*Clang") + add_compile_options(-fcolor-diagnostics) +elseif(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + add_compile_options(-fdiagnostics-color=always) +else() + message(STATUS "No colored compiler diagnostic set for '${CMAKE_CXX_COMPILER_ID}' compiler.") +endif() + diff --git a/cmake/StaticAnalyzers.cmake b/cmake/StaticAnalyzers.cmake new file mode 100644 index 0000000..6dffdf8 --- /dev/null +++ b/cmake/StaticAnalyzers.cmake @@ -0,0 +1,41 @@ +option(ENABLE_CPPCHECK "Enable static analysis with cppcheck" OFF) +option(ENABLE_CLANG_TIDY "Enable static analysis with clang-tidy" OFF) +option(ENABLE_INCLUDE_WHAT_YOU_USE "Enable static analysis with include-what-you-use" OFF) + +if(ENABLE_CPPCHECK) + find_program(CPPCHECK cppcheck) + if(CPPCHECK) + set(CMAKE_CXX_CPPCHECK + ${CPPCHECK} + --suppress=missingInclude + --enable=all + --inline-suppr + --inconclusive) + if(WARNINGS_AS_ERRORS) + list(APPEND CMAKE_CXX_CPPCHECK --error-exitcode=2) + endif() + else() + message(SEND_ERROR "cppcheck requested but executable not found") + endif() +endif() + +if(ENABLE_CLANG_TIDY) + find_program(CLANGTIDY clang-tidy) + if(CLANGTIDY) + set(CMAKE_CXX_CLANG_TIDY ${CLANGTIDY} -extra-arg=-Wno-unknown-warning-option) + if(WARNINGS_AS_ERRORS) + list(APPEND CMAKE_CXX_CLANG_TIDY -warnings-as-errors=*) + endif() + else() + message(SEND_ERROR "clang-tidy requested but executable not found") + endif() +endif() + +if(ENABLE_INCLUDE_WHAT_YOU_USE) + find_program(INCLUDE_WHAT_YOU_USE include-what-you-use) + if(INCLUDE_WHAT_YOU_USE) + set(CMAKE_CXX_INCLUDE_WHAT_YOU_USE ${INCLUDE_WHAT_YOU_USE}) + else() + message(SEND_ERROR "include-what-you-use requested but executable not found") + endif() +endif() diff --git a/conanfile.py b/conanfile.py new file mode 100644 index 0000000..34f4d7d --- /dev/null +++ b/conanfile.py @@ -0,0 +1,14 @@ +from conans import ConanFile + + +class CppStarterProject(ConanFile): + name = "IntegratorBD" + version = "0.1" + requires = ( + "catch2/2.13.7", + "docopt.cpp/0.6.2", + "fmt/8.0.1", + "spdlog/1.9.2", + "eigen/3.4.0" + ) + generators = "cmake", "gcc", "txt", "cmake_find_package" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..9e77618 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,25 @@ +find_package(fmt) +find_package(spdlog) +find_package(Eigen3) + +# Generic test that uses conan libs +add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp + Simulation.h Integratoren2d_forceless.cpp + Calculation.cpp Calculation.h Compute.cpp Compute.h) +target_include_directories(integratoren PUBLIC .) +target_link_libraries( + integratoren + PRIVATE project_options + project_warnings + fmt::fmt + spdlog::spdlog + Eigen3::Eigen3) + +add_executable(main main.cpp) +target_link_libraries( + main + PRIVATE project_options + project_warnings + integratoren + Eigen3::Eigen3 +) \ No newline at end of file diff --git a/src/Calculation.cpp b/src/Calculation.cpp new file mode 100644 index 0000000..dc107e4 --- /dev/null +++ b/src/Calculation.cpp @@ -0,0 +1,37 @@ +// +// Created by jholder on 22.10.21. +// + +#include "Calculation.h" + +#include + + +void Calculation::run(size_t steps) { + for (size_t step = 0; step < steps; ++step) { + m_integrator(rod, sim); + for (auto &comp: computes) { comp.eval(rod); } + } +} + +const Rod2d &Calculation::getRod() const { + return rod; +} + +Calculation::Calculation(std::function t_integrator, + std::initializer_list> t_computes, double deltaT, + size_t seed) + : sim(deltaT, seed), m_integrator( + std::move(t_integrator)) { + for (const auto &pair: t_computes) { + computes.emplace_back(rod, pair.first, pair.second); + } +} + +const std::vector &Calculation::getComputes() const { + return computes; +} + +const Simulation &Calculation::getSim() const { + return sim; +} diff --git a/src/Calculation.h b/src/Calculation.h new file mode 100644 index 0000000..cb4e135 --- /dev/null +++ b/src/Calculation.h @@ -0,0 +1,36 @@ +// +// Created by jholder on 22.10.21. +// + +#ifndef MYPROJECT_CALCULATION_H +#define MYPROJECT_CALCULATION_H + +#include "Rod2d.hpp" +#include "Simulation.h" +#include "functional" +#include "Compute.h" + + +class Calculation { +private: + Rod2d rod = Rod2d(1.0); + Simulation sim; + std::function m_integrator; + std::vector computes; +public: + const Simulation &getSim() const; + +public: + explicit Calculation(std::function t_integrator, + std::initializer_list> t_computes, double deltaT, + size_t seed); + + [[nodiscard]] const std::vector &getComputes() const; + + void run(size_t steps); + + [[nodiscard]] const Rod2d &getRod() const; +}; + + +#endif //MYPROJECT_CALCULATION_H diff --git a/src/Compute.cpp b/src/Compute.cpp new file mode 100644 index 0000000..28e6aab --- /dev/null +++ b/src/Compute.cpp @@ -0,0 +1,69 @@ +// +// Created by jholder on 22.10.21. +// + +#include "Compute.h" + +#include +#include +#include "Rod2d.hpp" + +void Compute::evalMSD(const Rod2d &rod2D) { + const auto &new_pos = rod2D.getPos(); + auto old_pos = start_rod.getPos(); + auto msd = (new_pos - old_pos).squaredNorm(); + agg.feed(msd); +} + +void Compute::evalOAF(const Rod2d &rod2D) { + const auto &new_e = rod2D.getE(); + auto old_e = start_rod.getE(); + auto oaf = old_e.dot(new_e); + agg.feed(oaf); +} + +void Compute::eval_empXX(const Rod2d & /*rod2D*/) { + //TODO +} + +void Compute::eval_empYY(const Rod2d & /*rod2D*/) { + //TODO +} + +void Compute::eval(const Rod2d &rod2D) { + time_step++; + if (time_step % every == 0){ + switch (type) { + case Type::msd: + evalMSD(rod2D); + break; + case Type::oaf: + evalOAF(rod2D); + break; + case Type::empxx: + eval_empXX(rod2D); + break; + case Type::empyy: + eval_empYY(rod2D); + break; + } + start_rod = rod2D; + } + +} + +Compute::Compute(Rod2d rod, Type t_type, size_t t_every) + : start_rod(std::move(rod)), agg({}), every(t_every), time_step(0), type(t_type) { +} + +const LiveAgg &Compute::getAgg() const { + return agg; +} + +Compute::Type Compute::getType() const { + return type; +} + +double Compute::getTime() const { + return static_cast(every); +} diff --git a/src/Compute.h b/src/Compute.h new file mode 100644 index 0000000..bdd24f2 --- /dev/null +++ b/src/Compute.h @@ -0,0 +1,46 @@ +// +// Created by jholder on 22.10.21. +// + +#ifndef MYPROJECT_COMPUTE_H +#define MYPROJECT_COMPUTE_H + + +#include "LiveAgg.hpp" +#include "Rod2d.hpp" + +class Compute { + +public: + enum class Type { + msd, oaf, empxx, empyy + }; + + explicit Compute(Rod2d rod, Type t_type, size_t t_every); + + void eval(const Rod2d &rod2D); + + void evalMSD(const Rod2d &rod2D); + + void evalOAF(const Rod2d &rod2D); + + void eval_empXX(const Rod2d &rod2D); + + void eval_empYY(const Rod2d &rod2D); + + [[nodiscard]] const LiveAgg &getAgg() const; + + [[nodiscard]] Type getType() const; + + double getTime() const; + +private: + Rod2d start_rod; + LiveAgg agg; + size_t every; + size_t time_step; + Type type; +}; + + +#endif //MYPROJECT_COMPUTE_H diff --git a/src/Integratoren2d_forceless.cpp b/src/Integratoren2d_forceless.cpp new file mode 100644 index 0000000..520785b --- /dev/null +++ b/src/Integratoren2d_forceless.cpp @@ -0,0 +1,141 @@ +// +// Created by jholder on 21.10.21. +// + +#include "Integratoren2d_forceless.h" + +Eigen::Vector2d ortho(Eigen::Vector2d e) { + return Eigen::Vector2d{-e[1], e[0]}; +} + +Eigen::Matrix2d mat(Eigen::Vector2d e) { + Eigen::Matrix2d matrix; + matrix << e, ortho(e); + + return matrix; +} + +void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { + auto trans_lab = Eigen::Vector2d({0, 0.01}); + rod2D.setPos(rod2D.getPos() + trans_lab); + Eigen::Rotation2D rot(0.1); + rod2D.setE(rot * rod2D.getE()); +} + +void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) { + auto std = sim.getSTD(); // sqrt(2*delta_T) + //translation + Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise + Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; + + //rotation + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + Eigen::Vector2d e = rod2D.getE(); + Eigen::Vector2d new_e = e + rot * ortho(e); + new_e.normalize(); + //apply + rod2D.setPos(rod2D.getPos() + trans_lab); + rod2D.setE(new_e); + +} + +void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { + // Predict with Euler + auto std = sim.getSTD(); + + Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; + Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab; + //rotation + double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + Eigen::Vector2d e = rod2D.getE(); + Eigen::Vector2d delta_e_predict = rot_predict * ortho(e); + Eigen::Vector2d e_predict = (e + delta_e_predict).normalized(); + + + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); + Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); + //apply + rod2D.setPos(pos_predictor); + rod2D.setE(e_integrated.normalized()); +} + +constexpr double kBT = 1.0; + +void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { + auto std = sim.getSTD(); + //translation + Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise + Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; + + //rotation + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); + Eigen::Vector2d e = rod2D.getE(); + auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); + + Eigen::Vector2d delta_e = correction * e + rot * ortho(e); + + //apply + rod2D.setPos(rod2D.getPos() + trans_lab); + rod2D.setE((e + delta_e).normalized()); + +} + +const Eigen::Vector2d Integratoren2d_forceless::unitVec = {1, 0}; + +//this is a slow implementation to keep the same data structure for performance analysis this must be altered +void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { + auto std = sim.getSTD(); + //translation + auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE()))); + auto trans_lab = rotMat * trans_part; + + auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); + Eigen::Rotation2Dd rotation(rot); + auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); + + // Normalisation should not be necessary if a proper angular representation is used. + // But with vector e it is done in case of numerical errors + + rod2D.setPos(rod2D.getPos() + trans_lab); + rod2D.setE(e_new); +} + +Eigen::Matrix2d Integratoren2d_forceless::e_2_matrix(Eigen::Vector2d m_e) { + Eigen::Matrix2d mat; + mat << m_e[0], -m_e[1], + m_e[1], m_e[0]; + return mat; +} + +void Integratoren2d_forceless::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) { + /* + auto std = sim.getSTD(); + + auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + auto trans_lab = rod2D.getE_Base_matrix() * trans_part; + auto pos_predictor = rod2D.getPos() + trans_lab; + + //rotation + auto rot_predict = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt());// rotationsmatrix verwenden. + auto e = rod2D.getE(); + auto delta_e_predict = rot_predict * ortho(e); + auto e_predict = (e + delta_e_predict).normalized(); + + auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() * rod2D.getDiff(); //old + new diff + + auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); + auto e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); + //apply + rod2D.setPos(pos_predictor); + rod2D.setE(e_integrated.normalized()); + */ +} + + diff --git a/src/Integratoren2d_forceless.h b/src/Integratoren2d_forceless.h new file mode 100644 index 0000000..5cf70b1 --- /dev/null +++ b/src/Integratoren2d_forceless.h @@ -0,0 +1,30 @@ +// +// Created by jholder on 21.10.21. +// + +#ifndef MYPROJECT_INTEGRATOREN2D_FORCELESS_H +#define MYPROJECT_INTEGRATOREN2D_FORCELESS_H +#include "Rod2d.hpp" +#include "Simulation.h" + +class Integratoren2d_forceless { +public: + static void Set1_Euler(Rod2d &rod2D, Simulation &sim); + + static void Set2_Heun(Rod2d &rod2D, Simulation &sim); + + static void Set3_Exact(Rod2d &rod2D, Simulation &sim); + + static void Set4_BDAS(Rod2d &rod2D, Simulation &sim); + + static const Eigen::Vector2d unitVec; + + static Eigen::Matrix2d e_2_matrix(Eigen::Vector2d m_e) ; + + static void Set5_MBD(Rod2d &rod2D, Simulation &sim); + + static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); +}; + + +#endif //MYPROJECT_INTEGRATOREN2D_FORCELESS_H diff --git a/src/LiveAgg.cpp b/src/LiveAgg.cpp new file mode 100644 index 0000000..90bf8a1 --- /dev/null +++ b/src/LiveAgg.cpp @@ -0,0 +1,33 @@ +#include "LiveAgg.hpp" +#include +LiveAgg::LiveAgg() : num_of_data(0), mean(0.0), S(0.0) +{ +} + +double LiveAgg::getSD() const noexcept +{ + return S / (num_of_data - 1); +} + +double LiveAgg::getSEM() const noexcept +{ + return S / (num_of_data - 1) / std::sqrt(num_of_data); +} + +double LiveAgg::getMean() const noexcept +{ + return mean; +} + +void LiveAgg::feed(double value) noexcept +{ + num_of_data += 1; + auto delta = value - mean; + mean += delta / num_of_data; + auto delta2 = value - mean; + S += delta * delta2; +} +int LiveAgg::getNumPoints() const noexcept +{ + return num_of_data; +} diff --git a/src/LiveAgg.hpp b/src/LiveAgg.hpp new file mode 100644 index 0000000..e15271f --- /dev/null +++ b/src/LiveAgg.hpp @@ -0,0 +1,20 @@ +#ifndef LIVEAGG_H +#define LIVEAGG_H + + +class LiveAgg +{ + public: + LiveAgg(); + [[nodiscard]] double getSD() const noexcept; + [[nodiscard]] double getSEM() const noexcept; + [[nodiscard]] double getMean() const noexcept; + [[nodiscard]] int getNumPoints() const noexcept; + void feed(double value) noexcept; + private: + int num_of_data; + double mean; + double S; +}; + +#endif // LIVEAGG_H diff --git a/src/Rod2d.cpp b/src/Rod2d.cpp new file mode 100644 index 0000000..517de77 --- /dev/null +++ b/src/Rod2d.cpp @@ -0,0 +1,76 @@ +#include "Rod2d.hpp" +#include "Eigen/Dense" + +constexpr double M_Pl = 3.141592653589793238462643383279502884; /* pi */ + +void sqrt(Eigen::Matrix2d &mat){ + const auto size = static_cast(mat.size()); + for (size_t i = 0; i < size; i++) { + mat.data()[i] = sqrt(mat.data()[i]); + } +} +Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) { + assert(L>=1.0); + if (L == 1.0) { + m_D_rot = 3.0; + m_Diff << 1, 0, + 0, 1; + } else { + const double D0 = 3.0 * M_Pl / L; + const double p = L; + auto D_para = D0 / (M_Pl * 2.0) * (log(p) - 0.207 + 0.980 / p - 0.133 / (p * p)); + auto D_ortho = D0 / (M_Pl * 2.0) * (log(p) + 0.839 + 0.185 / p + 0.233 / (p * p)); + m_D_rot = 3 * D0 / (M_Pl * L * L) * (log(p) - 0.662 + 0.917 / p - 0.050 / (p * p)); + m_Diff << D_para, 0, + 0, D_ortho; + } + m_Diff_sqrt = m_Diff; + m_D_rot_sqrt = sqrt(m_D_rot); + sqrt(m_Diff_sqrt); +} + +void Rod2d::reset() { + m_pos = {0, 0}; + m_e = {1, 0}; +} + +void Rod2d::setPos(const Eigen::Vector2d &Pos) { + m_pos = Pos; +} + +double Rod2d::getDRot() const { + return m_D_rot; +} + +const Eigen::Vector2d &Rod2d::getPos() const { + return m_pos; +} + +const Eigen::Matrix2d &Rod2d::getDiff() const { + return m_Diff; +} + +const Eigen::Matrix2d &Rod2d::getDiff_Sqrt() const { + return m_Diff_sqrt; +} + +double Rod2d::getDRot_Sqrt() const { + return m_D_rot_sqrt; +} + +const Eigen::Vector2d &Rod2d::getE() const { + return m_e; +} + +void Rod2d::setE(const Eigen::Vector2d &mE) { + m_e = mE; +} + +Eigen::Matrix2d Rod2d::getE_Base_matrix() const { + Eigen::Matrix2d mat; + mat << m_e[0], -m_e[1], + m_e[1], m_e[0]; + return mat; +} + + diff --git a/src/Rod2d.hpp b/src/Rod2d.hpp new file mode 100644 index 0000000..bc2beeb --- /dev/null +++ b/src/Rod2d.hpp @@ -0,0 +1,37 @@ +#pragma once +#include +#include +class Rod2d +{ +public: + explicit Rod2d(double L); + void reset(); + + [[nodiscard]] const Eigen::Matrix2d &getDiff() const; + + [[nodiscard]] double getDRot() const; + + [[nodiscard]] const Eigen::Vector2d &getPos() const; + + void setPos(const Eigen::Vector2d &Pos); + + [[nodiscard]] const Eigen::Matrix2d &getDiff_Sqrt() const; + + [[nodiscard]] const Eigen::Vector2d &getE() const; + + void setE(const Eigen::Vector2d &mE); + + [[nodiscard]] double getDRot_Sqrt() const; + + Eigen::Matrix2d getE_Base_matrix() const; + +private: + Eigen::Matrix2d m_Diff; + Eigen::Matrix2d m_Diff_sqrt; + double m_D_rot; + double m_D_rot_sqrt; + Eigen::Vector2d m_pos; // position + Eigen::Vector2d m_e; +}; + + diff --git a/src/Simulation.cpp b/src/Simulation.cpp new file mode 100644 index 0000000..176e752 --- /dev/null +++ b/src/Simulation.cpp @@ -0,0 +1,21 @@ +// +// Created by jholder on 21.10.21. +// + +#include "Simulation.h" + +double Simulation::getNorm(double t_norm) { + return t_norm*m_norm(m_generator); +} + +Simulation::Simulation(double t_delta_T, size_t seed): m_delta_T(t_delta_T),m_std(std::sqrt(t_delta_T*2.0)), m_generator(seed), m_norm(0, 1.0) { + +} + +double Simulation::getMDeltaT() const { + return m_delta_T; +} + +double Simulation::getSTD() const { + return m_std; +} diff --git a/src/Simulation.h b/src/Simulation.h new file mode 100644 index 0000000..8edc09a --- /dev/null +++ b/src/Simulation.h @@ -0,0 +1,26 @@ +// +// Created by jholder on 21.10.21. +// + +#ifndef MYPROJECT_SIMULATION_H +#define MYPROJECT_SIMULATION_H + + +#include + +class Simulation { +public: + explicit Simulation(double t_delta_T, size_t seed); + double getNorm(double t_norm); +private: + double m_delta_T; + double m_std; + std::mt19937_64 m_generator; + std::normal_distribution m_norm; +public: + [[nodiscard]] double getMDeltaT() const; + [[nodiscard]] double getSTD() const; +}; + + +#endif //MYPROJECT_SIMULATION_H diff --git a/src/legacy/Integratoren.cpp b/src/legacy/Integratoren.cpp new file mode 100644 index 0000000..3234c6e --- /dev/null +++ b/src/legacy/Integratoren.cpp @@ -0,0 +1,130 @@ +#include "Integratoren.hpp" +constexpr double kBT = 1.0; + +void Integratoren::integrateITOEuler_1(Rod2d &rod) +{ + //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); + double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); + double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + + double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; + double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; + + double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); + double tmpe1 = rod.e1 - rodOrtho * rod.e2; + double tmpe2 = rod.e2 + rodOrtho * rod.e1; + + rod.x1 = tmpx1; + rod.x2 = tmpx2; + rod.e1 = tmpe1; + rod.e2 = tmpe2; + rod.normalize(); +} + +void Integratoren::integrateITOHeun_2(Rod2d &rod) +{ + //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); + + double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); + double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + + double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; + double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; + + + double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); + double e1 = rod.e1 - rodOrtho * rod.e1; + double e2 = rod.e2 + rodOrtho * rod.e2; + double norm = std::sqrt(e1 * e1 + e2 * e2); + e1 /= norm; + e2 /= norm; + + double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT); + double rodOrthoRand = rod.dis(rod.generator); + double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrthoRand;//TODO Check if 0.5 faktor + double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrthoRand; + + rod.x1 = tmpx1; + rod.x2 = tmpx2; + rod.e1 = tmpe1; + rod.e2 = tmpe2; + rod.normalize(); +} + + +void Integratoren::integrateExact_3(Rod2d &rod) +{ + //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); + + double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); + double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + + double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; + double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; + + double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); + double korr = -0.5 * rod.Drot * rod.Drot / kBT / kBT * rod.deltaT; + double tmpe1 = rod.e1 - rodOrtho * rod.e2 + korr * rod.e1; + double tmpe2 = rod.e2 + rodOrtho * rod.e1 + korr * rod.e2; + + rod.x1 = tmpx1; + rod.x2 = tmpx2; + rod.e1 = tmpe1; + rod.e2 = tmpe2; + rod.normalize(); +} + +void Integratoren::integrateBDAS_4(Rod2d &rod) +{ + double deltaRTpara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); + double deltaRTortho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + double tmpx1 = rod.x1 + deltaRTpara * rod.e1 - deltaRTortho * rod.e2; + double tmpx2 = rod.x2 + deltaRTpara * rod.e2 + deltaRTortho * rod.e1; + + double deltaPhi = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); + double tmpe1 = rod.e1 * std::cos(deltaPhi) - rod.e2 * std::sin(deltaPhi); + double tmpe2 = rod.e1 * std::sin(deltaPhi) + rod.e2 * std::cos(deltaPhi); + + rod.x1 = tmpx1; + rod.x2 = tmpx2; + rod.e1 = tmpe1; + rod.e2 = tmpe2; + rod.normalize(); +} + +void Integratoren::integrateMBD_5(Rod2d &rod) +{ + /* + double movePara = std::sqrt(2.0*rod.Dpara*rod.deltaT)*rod.dis(rod.generator); // Taken from euler + double moveOrtho = std::sqrt(2.0*rod.Dortho*rod.deltaT)*rod.dis(rod.generator); // taken from euler + double x1 = rod.x1 + movePara*rod.e1- moveOrtho*rod.e2; + double x2 = rod.x2 + movePara*rod.e2+ moveOrtho*rod.e1; + */ + // unused without force + double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); + double e1 = rod.e1 - rodOrtho * rod.e2; + double e2 = rod.e2 + rodOrtho * rod.e1; + double norm = std::sqrt(e1 * e1 + e2 * e2); + e1 /= norm; + e2 /= norm; + + double movePara1 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D + double moveOrtho1 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + double movePara2 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D + double moveOrtho2 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); + double tmpx1 = rod.x1 + (movePara1 * rod.e1 - moveOrtho1 * rod.e2) + (movePara2 * e1 - moveOrtho2 * e2); + double tmpx2 = rod.x2 + (movePara1 * rod.e2 + moveOrtho1 * rod.e1) + (movePara2 * e2 + moveOrtho2 * e1); + + + double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT); + double rodOrtho2 = rod.dis(rod.generator); + + double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrtho2; + double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrtho2; + + rod.x1 = tmpx1; + rod.x2 = tmpx2; + rod.e1 = tmpe1; + rod.e2 = tmpe2; + rod.normalize(); +} diff --git a/src/legacy/Integratoren.hpp b/src/legacy/Integratoren.hpp new file mode 100644 index 0000000..bcb0168 --- /dev/null +++ b/src/legacy/Integratoren.hpp @@ -0,0 +1,19 @@ +#ifndef INTEGRATOREN_HPP +#define INTEGRATOREN_HPP +#include "Rod2d.hpp" + + +class Integratoren +{ +public: + static void integrateITOEuler_1(Rod2d &rod); + + static void integrateITOHeun_2(Rod2d &rod); + + static void integrateExact_3(Rod2d &rod); + + static void integrateBDAS_4(Rod2d &rod); + + static void integrateMBD_5(Rod2d &rod); +}; +#endif// INTEGRATOREN_HPP diff --git a/src/legacy/Simulation2d.cpp b/src/legacy/Simulation2d.cpp new file mode 100644 index 0000000..3c9b5ed --- /dev/null +++ b/src/legacy/Simulation2d.cpp @@ -0,0 +1,6 @@ +#include "Simulation2d.hpp" +#include +#include +#include "LiveAgg.hpp" +template +std::vector simulate(double deltaT, int seed, int meanSteps); diff --git a/src/legacy/Simulation2d.hpp b/src/legacy/Simulation2d.hpp new file mode 100644 index 0000000..2c66b4b --- /dev/null +++ b/src/legacy/Simulation2d.hpp @@ -0,0 +1,157 @@ +#ifndef SIMULATION2D_HPP +#define SIMULATION2D_HPP +#include "Rod2d.hpp" +#include "LiveAgg.hpp" +#include +#include +/* This is used to calculate the convergence order of the particles +*/ + +inline void printProgress(double progress) +{ + const int barWidth = 70; + std::cout << "["; + const int pos = static_cast(barWidth * progress); + for (int i = 0; i < barWidth; ++i) { + if (i < pos) { + std::cout << "="; + } else if (i == pos) { + std::cout << ">"; + } else { + std::cout << " "; + } + } + std::cout << "] " << int(progress * 100.0) << " % \r"; + std::cout.flush(); +} + +template +void simulate(double L, int seed, const std::string &filename) +{ + // we now simulate the msd after 10 steps + + + Rod2d rod(L, 1e-5, seed); + const int sampleWidth = 1000; + std::vector msd(sampleWidth, LiveAgg()); + std::vector oaf(sampleWidth, LiveAgg()); + std::vector empxx(sampleWidth, LiveAgg()); + std::vector empyy(sampleWidth, LiveAgg()); + std::vector xPos(sampleWidth, 0.0); + std::vector yPos(sampleWidth, 0.0); + std::vector e1(sampleWidth, 0.0); + std::vector e2(sampleWidth, 0.0); + std::vector targetMSD(sampleWidth, 0.0); + std::vector targetOAF(sampleWidth, 0.0); + std::vector targetEMPXX(sampleWidth, 0.0); + std::vector targetEMPYY(sampleWidth, 0.0); + std::vector samplePoints(sampleWidth, 1); + size_t temp = 1; + const double Dmean = 0.5 * (rod.Dpara + rod.Dortho); + const double deltaD = rod.Dortho - rod.Dpara; + for (size_t i = 0; i < samplePoints.size(); ++i) { + samplePoints[i] = temp; + temp += 100; + const double time = static_cast(temp) * rod.deltaT; + targetMSD[i] = 4.0 * 0.5 * (rod.Dpara + rod.Dortho) * time; + targetOAF[i] = std::exp(-rod.Drot * time); + const double u = 4.0; + targetEMPXX[i] = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time; + targetEMPYY[i] = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time; + } + unsigned long long int iteration = 0; + //const unsigned long long int maxItr = ((long int)10000)*1000000; + const unsigned long long int maxItr = 10000000; + const unsigned long long int updateInter = 10000; + while (iteration < maxItr) { + T(rod); + iteration++; + for (size_t i = 0; i < samplePoints.size(); ++i) { + const size_t sample = samplePoints[i]; + if (iteration % sample == 0) { + msd[i].feed((rod.x1 - xPos[i]) * (rod.x1 - xPos[i]) + (rod.x2 - yPos[i]) * (rod.x2 - yPos[i])); + oaf[i].feed(rod.e1 * e1[i] + rod.e2 * e2[i]); + empxx[i].feed(std::pow((rod.x1 - xPos[i]) * e1[i] + (rod.x2 - yPos[i]) * e2[i], 2)); + empyy[i].feed(std::pow((rod.x1 - xPos[i]) * e2[i] - (rod.x2 - yPos[i]) * e1[i], 2)); + xPos[i] = rod.x1; + yPos[i] = rod.x2; + e1[i] = rod.e1; + e2[i] = rod.e2; + } + } + if (iteration % updateInter == 0) { + printProgress(static_cast(iteration)/ maxItr); + } + } + + std::cout << "\n"; + + std::ofstream fileMSD(filename + "_" + std::to_string(L) + "_MSD.out"); + std::ofstream fileOAF(filename + "_" + std::to_string(L) + "_OAF.out"); + std::ofstream fileXX(filename + "_" + std::to_string(L) + "_XX.out"); + std::ofstream fileYY(filename + "_" + std::to_string(L) + "_YY.out"); + + for (size_t i = 0; i < samplePoints.size(); ++i) { + fileMSD << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetMSD[i] << "\t" << msd[i].getMean() << "\t" << msd[i].getSEM() << "\t" << msd[i].getSD() << "\t" << msd[i].getNumPoints() << "\n"; + fileOAF << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetOAF[i] << "\t" << oaf[i].getMean() << "\t" << oaf[i].getSEM() << "\t" << oaf[i].getSD() << "\t" << oaf[i].getNumPoints() << "\n"; + fileXX << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetEMPXX[i] << "\t" << empxx[i].getMean() << "\t" << empxx[i].getSEM() << "\t" << empxx[i].getSD() << "\t" << empxx[i].getNumPoints() << "\n"; + fileYY << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetEMPYY[i] << "\t" << empyy[i].getMean() << "\t" << empyy[i].getSEM() << "\t" << empyy[i].getSD() << "\t" << empyy[i].getNumPoints() << "\n"; + } + fileMSD.close(); + fileOAF.close(); + fileXX.close(); + fileYY.close(); +} + + +template +void konvergenz(double L, int seed, double deltaT, const std::string &filename) +{ + // we now simulate the msd after 10 steps + Rod2d rod(L, deltaT, seed); + const double time = 2; + const int steps = static_cast(time / deltaT); + const int minSteps = 1000; + + //const double targetMSD = 4.0*0.5*(rod.Dpara+rod.Dortho)*time; + const double targetOAF = std::exp(-rod.Drot * time); + /*const double u = 4.0; + const double Dmean = 0.5*(rod.Dpara+rod.Dortho); + const double deltaD = rod.Dortho - rod.Dpara; + const double targetEMPXX = Dmean-deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time; + const double targetEMPYY = Dmean+deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;*/ + //LiveAgg msd; + LiveAgg oaf; + //LiveAgg empxx; + //LiveAgg empyy; + int meanCount = 0; + while (meanCount < minSteps or std::abs(targetOAF - oaf.getMean()) < oaf.getSEM() * 10) { + //reset the position + rod.reset(); + for (int i = 0; i < steps; ++i) { + T(rod); + } + //msd.feed(rod.x1*rod.x1+rod.x2*rod.x2); + oaf.feed(rod.e1); + //.feed(rod.x1*rod.x1/2/time); + //empyy.feed(rod.x2*rod.x2/2/time); + if (meanCount % 100000 == 0) { + //std::cout< +#include +#include +#include "legacy/Integratoren.hpp" +#include "legacy/Simulation2d.hpp" +const int seed = 12345; +[[maybe_unused]] void overTime(int sim) +{ + const std::vector Length = { 1.0, 1.3974, 20.0 }; + const int numCases = 5; + int i = 0; + for (auto L : Length) { + if (sim == numCases * i) { + simulate(L, seed, "1"); + } + if (sim == numCases * i + 1) { + simulate(L, seed, "2"); + } + if (sim == numCases * i + 2) { + simulate(L, seed, "3"); + } + if (sim == numCases * i + 3) { + simulate(L, seed, "4"); + } + if (sim == numCases * i + 4) { + simulate(L, seed, "5"); + } + i++; + } +} + +[[maybe_unused]]void konvergenz(int sim) +{ + const std::vector Length = { 1.0, 1.3974, 20.0 }; + const std::vector DELTAT = { 2.0, 1.0, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01, 0.005, 0.0025, 0.200, 0.001, 0.0005, 0.00025 }; + int i = 0; + for (auto dt : DELTAT) { + + for (auto L : Length) { + if (sim == i++) { + konvergenz(L, seed, dt, "1"); + } + if (sim == i++) { + konvergenz(L, seed, dt, "2"); + } + if (sim == i++) { + konvergenz(L, seed, dt, "3"); + } + if (sim == i++) { + konvergenz(L, seed, dt, "4"); + } + if (sim == i++) { + konvergenz(L, seed, dt, "5"); + } + } + } +} + + +int main(int argc, char *argv[]) +{ + int sim; + auto args = std::span(argv, size_t(argc)); + if (argc > 1) { + sim = std::stoi(args[1]) - 1; + } else { + sim = 0; + } + overTime(sim); + //konvergenz(sim); + return 0; +} diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..a070fda --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,40 @@ +// +// Created by jholder on 22.10.21. +// + +#include +#include "Calculation.h" +#include "Integratoren2d_forceless.h" + +int main() { + constexpr int numStep = 1000000; + Calculation euler(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, + {Compute::Type::msd, 2}, + {Compute::Type::msd, 4}, + {Compute::Type::msd, 8}, + {Compute::Type::msd, 16}, + {Compute::Type::oaf, 1}, + {Compute::Type::oaf, 2}, + {Compute::Type::oaf, 4}, + {Compute::Type::oaf, 8}, + {Compute::Type::oaf, 16}}, 0.01, 12345); + euler.run(numStep); + for (const auto &com: euler.getComputes()) { + if(com.getType()!=Compute::Type::msd) + continue; + auto time = euler.getSim().getMDeltaT() * com.getTime(); + auto targetMSD = 4.0 * 0.5 * (euler.getRod().getDiff().trace()) * time; + std::cout << "MSD: " << com.getAgg().getMean() + << " +- " + << com.getAgg().getSEM() + << " " << targetMSD << std::endl; + } + for (const auto &com: euler.getComputes()) { + if(com.getType()!=Compute::Type::oaf) + continue; + auto time = euler.getSim().getMDeltaT() * com.getTime(); + auto targetOAF = std::exp(-euler.getRod().getDRot() * time); + std::cout <<"OAF: " << com.getAgg().getMean() + << " +- " << com.getAgg().getSEM() << " " << targetOAF << std::endl; + } +} \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt new file mode 100644 index 0000000..2d18e0c --- /dev/null +++ b/test/CMakeLists.txt @@ -0,0 +1,32 @@ +find_package(Catch2 REQUIRED) +find_package(Eigen3) + +include(CTest) +include(Catch) + +add_library(catch_main STATIC catch_main.cpp) +target_link_libraries(catch_main PUBLIC Catch2::Catch2) +target_link_libraries(catch_main PRIVATE project_options) + +set(TEST_FILES + test_LiveAgg.cpp + test_Rod2d.cpp + test_Calculation.cpp test_Compute.cpp test_Simulation.cpp) + +add_executable(tests ${TEST_FILES}) +target_link_libraries(tests PRIVATE project_warnings project_options catch_main Eigen3::Eigen3 integratoren) + +# automatically discover tests that are defined in catch based test files you can modify the unittests. Set TEST_PREFIX +# to whatever you want, or use different for different binaries +catch_discover_tests( + tests + TEST_PREFIX + "unittests." + REPORTER + xml + OUTPUT_DIR + . + OUTPUT_PREFIX + "unittests." + OUTPUT_SUFFIX + .xml) diff --git a/test/catch_main.cpp b/test/catch_main.cpp new file mode 100644 index 0000000..d8d2eca --- /dev/null +++ b/test/catch_main.cpp @@ -0,0 +1,5 @@ +#define CATCH_CONFIG_MAIN // This tells the catch header to generate a main + +#include + + diff --git a/test/test_Calculation.cpp b/test/test_Calculation.cpp new file mode 100644 index 0000000..8c24aac --- /dev/null +++ b/test/test_Calculation.cpp @@ -0,0 +1,52 @@ +// +// Created by jholder on 22.10.21. +// + +#include +#include +#include "LiveAgg.hpp" +#include "Integratoren2d_forceless.h" + +TEST_CASE("Basic run of Calculation") { + Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 10}, + {Compute::Type::oaf, 10}}, 1, 1); + Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 10}, + {Compute::Type::oaf, 10}}, 1, 1); + Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 10}, + {Compute::Type::oaf, 10}}, 1, 1); + Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 10}, + {Compute::Type::oaf, 10}}, 1, 1); + SECTION("Euler") { + euler.run(100); + CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0)); + CHECK(euler.getComputes()[0].getType() == Compute::Type::msd); + CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf); + }SECTION("Heun") { + heun.run(100); + CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0)); + }SECTION("Exact") { + exact.run(100); + CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0)); + }SECTION("Euler") { + bdas.run(100); + CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0)); + } + + SECTION("MSD") { + euler.run(100); + CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10); + CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10); + }SECTION("Deterministic") { + Calculation determ(Integratoren2d_forceless::Set0_deterministic, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 1, 1); + determ.run(10); + auto time = 1; + + auto targetMSD = pow(0.01, 2) * time; + auto targetOAF = cos(0.1); + + CHECK(determ.getComputes()[0].getAgg().getMean() == Catch::Detail::Approx(targetMSD)); + CHECK(determ.getComputes()[1].getAgg().getMean() == Catch::Detail::Approx(targetOAF)); + } +} diff --git a/test/test_Compute.cpp b/test/test_Compute.cpp new file mode 100644 index 0000000..12a951e --- /dev/null +++ b/test/test_Compute.cpp @@ -0,0 +1,19 @@ +// +// Created by jholder on 22.10.21. +// + +#include +#include "catch2/catch.hpp" +#include "Integratoren2d_forceless.h" + +TEST_CASE("Compute") { + auto rod = Rod2d(1.0); + auto com = Compute(rod, Compute::Type::msd, 10); + SECTION("Mean of same values") { + for (int i = 0; i < 100; ++i) { + com.eval(rod); + } + CHECK(com.getAgg().getNumPoints()==10); + CHECK(com.getAgg().getMean()==0.0); + } +} diff --git a/test/test_LiveAgg.cpp b/test/test_LiveAgg.cpp new file mode 100644 index 0000000..60c5021 --- /dev/null +++ b/test/test_LiveAgg.cpp @@ -0,0 +1,16 @@ +// +// Created by jholder on 21.10.21. +// + +#include +#include "LiveAgg.hpp" + +TEST_CASE("LiveAgg") { + auto ag = LiveAgg(); + SECTION("Mean of same values"){ + ag.feed(1.0); + ag.feed(1.0); + REQUIRE(ag.getMean()==1.0); + REQUIRE(ag.getNumPoints()==2); + } +} diff --git a/test/test_Rod2d.cpp b/test/test_Rod2d.cpp new file mode 100644 index 0000000..f3939d6 --- /dev/null +++ b/test/test_Rod2d.cpp @@ -0,0 +1,47 @@ +// +// Created by jholder on 21.10.21. +// + +#include "Rod2d.hpp" +#include +#include + +TEST_CASE("Rods") { + Rod2d sphere(1); + Rod2d rod(2); + SECTION("Checking translation") { + rod.setPos(rod.getPos() + Eigen::Vector2d(1, -1)); + auto newPos = rod.getPos(); + REQUIRE(newPos[0] == 1); + REQUIRE(newPos[1] == -1); + REQUIRE(sphere.getDiff_Sqrt()(0,0) == 1.0); + REQUIRE(sphere.getDiff_Sqrt()(0,1) == 0.0); + REQUIRE(sphere.getDiff_Sqrt()(1,0) == 0.0); + REQUIRE(sphere.getDiff_Sqrt()(1,1) == 1.0); + }SECTION("Checking rotation") { + SECTION("Forwards == e 90"){ + // + // ===== + // + rod.setE({0, 1}); // Lab system is 90 degree rotated + auto test = Eigen::Vector2d ({1,0}); + auto erg = (rod.getE_Base_matrix()*test).normalized(); + auto e = rod.getE(); + REQUIRE(erg.isApprox(e)); + + } + SECTION("Forwards == e 45"){ + // o + // o + // o + rod.setE(Eigen::Vector2d ({1,1}).normalized()); // Lab system is 90 degree rotated + auto test = Eigen::Vector2d ({1,0}); + auto erg = (rod.getE_Base_matrix()*test).normalized(); + auto e = rod.getE(); + REQUIRE(erg.isApprox(e)); + } + + + } + +} \ No newline at end of file diff --git a/test/test_Simulation.cpp b/test/test_Simulation.cpp new file mode 100644 index 0000000..e750a7d --- /dev/null +++ b/test/test_Simulation.cpp @@ -0,0 +1,27 @@ +// +// Created by jholder on 22.10.21. +// + +#include +#include "Simulation.h" + +TEST_CASE("Simulation") { + auto sim = Simulation(1,0); + REQUIRE(sim.getSTD()==sqrt(1*2.0)); + SECTION("Symmetrie"){ + constexpr int num = 10000000; + double sum = 0.0; + for (int i = 0; i < num; ++i) { + sum += sim.getNorm(1); + } + CHECK(sum/num == Catch::Detail::Approx(0.0).margin(0.001)); + } + SECTION("STD"){ + constexpr int num = 10000000; + double sum = 0.0; + for (int i = 0; i < num; ++i) { + sum += pow(sim.getNorm(1),2); + } + CHECK(sum/num == Catch::Detail::Approx(1.0).margin(0.001)); + } +} \ No newline at end of file