From b5e8ef350a7318f910b39ed8aa01637b43ac8fa0 Mon Sep 17 00:00:00 2001 From: Jacob Date: Wed, 22 Dec 2021 13:45:26 +0100 Subject: [PATCH] test --- Auswertung/notebooks/kovergence.ipynb | 165 +++++++++++++++ Auswertung/requirements.txt | 3 +- C++/src/CMakeLists.txt | 17 +- C++/src/Calculation.cpp | 94 --------- C++/src/Calculation.h | 57 ------ C++/src/Compute.cpp | 45 +++-- C++/src/Integratoren2d.h | 47 +++-- C++/src/LiveAgg.cpp | 6 +- C++/src/force_lambdas.h | 2 +- C++/src/konvergenz_runner.cpp | 51 +++++ C++/src/konvergenz_runner.h | 18 ++ C++/src/legacy/Integratoren2d_force.cpp | 184 ----------------- C++/src/legacy/Integratoren2d_force.h | 52 ----- C++/src/legacy/Integratoren2d_forceless.cpp | 123 ------------ C++/src/legacy/Integratoren2d_forceless.h | 25 --- C++/src/logging.cpp | 20 ++ C++/src/logging.hpp | 15 ++ C++/src/main.cpp | 210 +++----------------- C++/test/CMakeLists.txt | 4 +- C++/test/catch_main.cpp | 10 +- C++/test/test_Compute.cpp | 36 +++- C++/test/test_Integratoren.cpp | 62 ------ C++/test/test_Rod2d.cpp | 2 + C++/test/test_Simulation.cpp | 28 --- 24 files changed, 413 insertions(+), 863 deletions(-) create mode 100644 Auswertung/notebooks/kovergence.ipynb delete mode 100644 C++/src/Calculation.cpp delete mode 100644 C++/src/Calculation.h create mode 100644 C++/src/konvergenz_runner.cpp create mode 100644 C++/src/konvergenz_runner.h delete mode 100644 C++/src/legacy/Integratoren2d_force.cpp delete mode 100644 C++/src/legacy/Integratoren2d_force.h delete mode 100644 C++/src/legacy/Integratoren2d_forceless.cpp delete mode 100644 C++/src/legacy/Integratoren2d_forceless.h create mode 100644 C++/src/logging.cpp create mode 100644 C++/src/logging.hpp delete mode 100644 C++/test/test_Integratoren.cpp delete mode 100644 C++/test/test_Simulation.cpp diff --git a/Auswertung/notebooks/kovergence.ipynb b/Auswertung/notebooks/kovergence.ipynb new file mode 100644 index 0000000..af6414c --- /dev/null +++ b/Auswertung/notebooks/kovergence.ipynb @@ -0,0 +1,165 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 36, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 37, + "outputs": [], + "source": [ + "data1 = np.loadtxt(\"../data/1Euler_msd_harmonic.dat\")\n", + "data2 = np.loadtxt(\"../data/2Heun_msd_harmonic.dat\")\n", + "data3 = np.loadtxt(\"../data/3Exact_msd_harmonic.dat\")\n", + "data4 = np.loadtxt(\"../data/4BDAS_msd_harmonic.dat\")\n", + "data5 = np.loadtxt(\"../data/5MBD_msd_harmonic.dat\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 38, + "outputs": [ + { + "data": { + "text/plain": "[]" + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxJUlEQVR4nO3dZ1RU59rG8f+eoUqVIkrviKKigoDdWBNjbDG2lJOenJPkJCfV9Gb6m95PetFEjVFj1GjsBVHslS7SBASkw8DMfj+gOTGCsQAzw9y/tfzA8Dx77llr1sV2P01RVRUhhBAdn8bYBQghhGgfEvhCCGEhJPCFEMJCSOALIYSFkMAXQggLIYEvhBAWwsrYBTRHUZQJwAQnJ6fbw8PDjV2OEEKYjV27dp1UVdWzud8ppjwPPyYmRk1OTjZ2GUIIYTYURdmlqmpMc7+TRzpCCGEhJPCFEMJCSOALIYSFkMAXQggLYZKBryjKBEVRPi0vLzd2KUII0WGYZOCrqvqLqqp3uLi4GLsUIYToMEwy8C/Xio/nkLzmB2OXIYQQJsUkF15djqryUly+WILOsIQd+kYGjLve2CUJIYRJ6HB3+I4ubugf/RdaA+ifmMv2FV8ZuyQhhDAJHS7wAYZOvQfd0/ejUYEnXyVx2afGLkkIIYyuQwY+wJBJd2J44SFQQPPsW2z5+UNjlySEEEZlkoHfWtMyE8bfCnMfw6ABm+ffY9PC91qpQiGEMD8mGfitOS0zftxNaF9+mkYrsJ/7IRvm/18rVCiEEObHJAO/tQ0YPRObV55HZwOOr3zG+u9fNXZJQgjR7iwi8AFiRk6j0+tzqbUD59e+Yu3Xc41dkhBCtCuLCXyAvsOm4PTmq9TYQef/+47VXzxn7JKEEKLdWFTgA0QPugaXd96kqhN4vP0Dv336hLFLEkKIdtEhA3/D90fZ/Vs2Br2h2d/3jr8S93ffo9IBvN5bzKqPHmvnCoUQov11uMDX6w3UVjaQ+HMGP722i5K8qmbb9Rwwii4ffESZE3T9YCkr3n+wnSsVQoj2ZZKBfznz8LVaDePujGLMbT2pLK1jwUs72flrFvrGc+/2I/sNp9tH/6XUFXw+XsHyd/7dCtULIYRp6tCHmNdW6ti8II20nYW4+zhyxY3d6RLgfE67tAOJHP/nLXiUwvFbRzLhP+9fTtlCCGE0FnuIub2TDWNu7cmVd/WitkrHold3kbgkg8YG/VntwnolEPjJNxS7QcDna1n62l1GqlgIIdpOhw78M4KjPZn5dBwR8V3ZvSqbBXN3ciLz7MdFIT1iCftsPoUeCsFfbWTJK7cZqVohhGgbFhH4AHYO1oy8MZIJ9/ahoV7PT6/vYsvCNBp0/7vbD4iIJuLzHynoohDyzVaWzP2H8QoWQohWZjGBf4Z/T3dmPh1H1BAf9q3N4YcXdpCXUva/34f2IvKLReR5KYR+l8Ti528wYrVCCNF6LC7wAWzsrRg2K4JJD/QFVWXJW3vYOC8FXV0jAH7BPej99RJyvBXC5yez+OmZRq5YCCEun0UG/hk+EZ2Z8VQcfUb6cXBzHvOfT+L4oRIAvP3Dif7qF457K0Qs3MtPT15n5GqFEOLyWHTgA1jbahk8LYypD/fH2kbLL+/tY+03R6irbqCbXwgx363gmK+G7osOsGjOVGOXK4QQl8ziA/+MrsEuXPdELP3HBZCy/QTzn08ia18xXboFEv/9b2T5a4j8+TCLHplo7FKFEOKSSOD/iZW1lvhJIUx7LAZ7RxtWfHSA1Z8folOnLgz8/ncyAzT0XJbKwoeuNnapQghx0Uwy8FvriMNL5envxLQ5MQyYEETG7iLmP5dESbaGwfPWkhqkJWp5BgvvH2eU2oQQ4lJ16K0VWkNJXhXrvjlCUXYlwdGe9L3Kg+13jSEis5EDY/y47t3VRq1PCCH+zGK3VmgN7j6OTH2kPwmTQ8g+WMLyt9MIeHAxR0Ot6LU6hwX3jDR2iUIIcUEk8C+ARquh39gApj8ZS+euDmxbmIMy8ksOdfeg1+/5LLh7hLFLFEKIvyWBfxE6d3Vg8kP9GHxdGIVZVZQGvEBy9FCi1p9g4Z1DwYQfjwkhhAT+RdJoFPpc4ceMp+LoEuhMhet0tib8m5AkAwtuHyKhL4QwWRL4l8jF056J9/dl+OwIVOdIEuOewDUrih9vGYxBr//7CwghRDuTwL8MiqLQc4gPM5+Jw7eHF2lh19FQ9w9+vO0aCX0hhMmRwG8FTm52XHNfX4ZdH0aFUzdOWd3LN7c/SkN9g7FLE0KIP0jgtxJFUYga7MesucPR6lOptrmKL+6ZR1H2KWOXJoQQgAR+q3P1dOTmj+/GoeZ7FL0rC1/aQdIv6c0eoi6EEO1JAr8NWNtYc/3nn+Cgexuvoj0k/3qchS/toCi7wtilCSEsmAR+G7GytmLGV2vQOf5MrwOfUJZbwKJXk5s9RF0IIdqDBH4b0lpZMf2rreSE5TBo2wvYVSe1eIi6EEK0NQn8NqbVapnx+WYODOhEQtI3uJ36GF1dY7OHqAshRFuSwG8HWq2WmZ9tYtdgT6L3HsA55QEiEzybPURdCCHaikkGvrH3w28LGo2G2f/dyK6hXkQeqaJi0Syuuius2UPUhRCiLch++O1MVVW+/+co+q/PJzXMihFfrWf/ujL2rcvBsbMtI67vjn8Pd2OXKYQwU7IfvglRFIXrP1rLrpG+hKc1suGm4fQb6/a/Q9Tf3ce6b45QXyOrdIUQrUsC30iu/2ANyWMCCEvXs3bGMJw713HdE7H0GxfA0e0nmPdc0yHqQgjRWiTwjeiGd1ex68oQQjP1rJt1BXUVhSRMCuHaR/ufdYh6baXO2KUKIToACXwju/6t5SSPDyMky8DGWaOpKM6lS4DzWYeoz3s2iaOJBZjyeIsQwvRJ4JuAG/5vGcnXRBKYbWDL7DGUFWWjtdIQOz6I656IxdWrE2u/PsLSt/dyqqjG2OUKIcyUBL6JuOG1xeydHEXAcZXts6+kND8TAHdvR6Y81I9hsyIozq7ghxd2kLzymGzGJoS4aBL4JmT2ywvZO7UP/rkqO264mpK8dAAUjULUUB9mPRtPYJQ7SUszWfCSbM8ghLg4EvgmZvbcH9g7rT+++SrJN1zDyZyUP37n4GrLuDt7cdXdvdDVNm3PsHF+CvW1smBLCPH3JPBN0Oznv2P/jDi8T6jsvnEyBalnLz4L6uPJzGfi6D3Cl4Ob8pj/7HYy98gUTiHE+Ungm6hZz3zFwVkD6VaocmzWDax840740ywdGzsrhlwXzrWPxGDnZMPKTw6w4qP9VJXVGbFqIYQpk60VTNzK+W9j+/4ndCuBI71sGPzG93QJiDqrjV5vYN/vOexcnoWiVYifGELUMB80GsVIVQshjOV8WytI4JuBktJifrt/In12lFHmAtU3jGbsPe+e0668uJaN81PIOVyKV5Azw2d3x8PX0QgVCyGMRfbSMXPubp7M+mYbWY/eiEEB//fXsGBWX0rzU89q5+Jpz4R7+zDq5h5UnKxl4Us7Sfw5g0bZc18Igdzhm53iwjzW3T+V3nvKKewMupuvZswdr5/Trq6qga2L0zm6rQBnT3uGz4rAL9LNCBULIdqT3OF3IJ5ePkyfv53Uf1+HdSP4vLWcH27sz6nCrLPa2TlaM/LGSCY+0BdFgWXv7OX3Lw/LvjxCWDC5wzdjBbnZbL3/WnoerCLfA7h9KiNvevGcdo0NenatzGb3qmxs7K0YNC2UiLiuKIoM6grR0cgdfgfVzTeAaxftJOXOiTjUQpdXf2LerbFUluad1c7KWkvcNcGn9+WxZ+1XR1j2juzLI4SlkTv8DiInM5XkB2bSPaWG413B9u7ZDJ/+5DntVIPKoS35JC5OR69XiR0fSPRof7Ra+dsvREcg0zItyOKX/43fgtXY6ODAEBcmv7oEB5eu57SrPlXP5h9TydhTjJu3AyOu707XYBcjVCyEaE0S+BYm6+hBDj54A6EZdWT6gPO9tzBk0sPNt91XzKYfUqk6VU+voT7ETwrBxt6qnSsWQrQWCXwLZDAYWPrCP/FfvBGNAQ6McGPaS0uxd/Q4p62urpGkpZns35CLg7MNQ2dEENzX0whVCyEul0kEvqIoDsCHgA7YoKrq93/XRwL/8qXvTyb1oVsJOq4jNQC6/PufJFx1b7NtC7MqWP/9UUpyqwjq48HQGeE4drZr54qFEJejzWbpKIryhaIoRYqiHPzL6+MURUlRFCVdUZTHTr88BVikqurtwDWX877iwoX2jmHcqj0cnhyPfz7YzPmQr/8zlPraU+e09QpqOloxYXIIOYdLmfdcEgc25GIwmO7/AoUQF+5yp2Z8BYz78wuKomiBD4ArgR7ATEVRegC+QM7pZrLWvx1pNBqmvvwldl98RomHNQNWFPPLtQnsWPvpOW21Wg39xgYw4+k4ugY5s+mHVBa/vouTuVVGqFwI0ZouK/BVVd0ElP7l5QFAuqqqmaqq6oAfgIlALk2hf973VRTlDkVRkhVFSS4ulj3eW1Nk7CBG/babI1f1JfwY8NBbfPXYCBrqKs5p6+Jpz4T7ohl1cw/Ki2VfHiE6graYfO3D/+7koSnofYDFwFRFUT4Cfmmps6qqn6qqGqOqaoynpwwctjaNlRVT3pyH5uMPqHTWErfkBIunx5O8+etz2iqKQkRcV2Y/G094fFd2/5bN/Bd2kHPkr3/jhRDmoC0Cv7n1+qqqqtWqqt6squrdFzJgK9pW1JArGL5mD0dG9qRHmoru/lf48unRNNSf++jmj3157o9GQfblEcJctUXg5wJ+f/rZF8hvg/cRl0lrY82UDxahf+cNdPYa4hfk8uOsAezbvqDZ9r7d3Zjx1AD6XxlA2s5C5j2bxNHtBZjy1F4hxP9c9rRMRVECgeWqqkad/tkKSAVGAnnATmCWqqqHLuKaE4AJoaGht6elpV1WfeLC6GrrWHHfdMK2pFLqBBkTgrjh0Z/Q2tg3274kr4oN3x/lRGYFvt07M2xWBK5dOrVz1UKIv2qzefiKoswHhgMeQCHwjKqqnyuKchXwNqAFvlBVde6lXF/m4be/XcsXUTP3GTzKDCRFK8Q+/CI9+09ptq1qUDm0OY/EnzNkXx4hTIRJLLy6FBL4xlFfWcWqe6YRnnSME50he1IYNz64EI2VbbPtZV8eIUyHbI8sLoqtkyMTv15J5QuPY6tXiP0qjS//0Y8jB1Y2297B1ZZxd/biqrt7oatt5KfXd7Fpfgq62sZ2rlwIcT4meYcvz/BNR+2pCn6/ewqhe/LI9YD8KT254b55aKxsmm1/1r48LrYMnRFOcLRMrxWivcgjHXHZEr/9L5q336JTrcq2OC1XPPIOYZEjW2xfmFXB+u+OUpIn+/II0Z4k8EWrqDpZwsa7phJ8sJBjXnByWh+uv/s7FG3z2ynr9Qb2/Z7DzuVZKFqF+IkhRA3zQaORoxWFaCsS+KJVbfn0XWw/+ghbHWweqOWqRz8kKHRoi+3Li2vZOO8oOUfKcPd1ZOj0cLzDXNuvYCEsiAS+aHUV+QVsvXsagSklpPlA5fQBzLr1CxStttn2qqqSsbuYrYvSqCqrJyzWi4FTQuQxjxCtTAJftAlVVdn83ms4fv4VigG2DLFm0sOf4BeU0GKfBp2e3b9ls+e34yga6H9lINGj/LCybv4PhRDi4phd4MssHfNyKvs4O+6ejl/mKQ4HgG7GIGbc+EmLd/sAFSdr2fpTOpl7inH2sGPQtWEE9fFAUeT5vhCXw+wC/wy5wzcfqqqy8bXncf3uB/QKbB5hw7SHv8DHt/95++UcLWXzj2mUFVTj18ONwdPCcOvm0E5VC9HxSOCLdnMyLYN9/5yJd04l+0JAM3M402Z9gKJpeY2fXm/g4MY8dvySRWO9nl5X+BI7PghbOUxdiIsmgS/alWowsOGFx3FfsJQ6a9g80pbrH/wGL+/e5+1XW6lj+9JMDm/Nx97RmoTJIXSP74Yi0ziFuGAS+MIoCg8c5vB9N9C1oIbkCLCfNYap17553mf7AEXZFWz+MZUTmRV0CXBiyIxwugbJ3jxCXAizC3wZtO041MZG1j/1EJ5Lf6PKDrYOsWL8zU8RGX3d+fupKqk7Ctm2OJ2ach3dE7oSPykEB5fmN3ATQjQxu8A/Q+7wO478XXtIfegOvAqqSPOGY1e4cfPtn+Lk1fO8/XR1jexaeYy9v+egtdYQOz6I3iN80VrJvn9CNEcCX5gE1WBg5ycfwX8/xKHGQGIvcBrbi+tmf4pi73revqcKa9iyKI3sAyW4enVi8HVhBPR0b5/ChTAjEvjCpDRUVLLu8QfwXreVemvYkgBDJk+n/6gnoYV9ec44duAkWxamUV5US2BvDwZPC8XFU07aEuIMCXxhkk4eSWH3Q3fil1FIjgekDNUyc9bzeEQ1f8LWGfoGA/vW5ZC84hh6vYG+o/zpNy4AGzuZximEBL4wWaqqcmDRT1S+8QJu5Tp2hQPD3Jkx6yOsuvU6b9/q8noSf84gZfsJHFxsGDg1lLBYL1mtKyyaBL4wefq6Ota++CxdlixFAbbEqPQe2Y8hk98Bx/MfoHIis5xNP6RSfLySbqEuDLkuHE9/p/YpXAgTY3aBL9MyLVfl8Vw2PfwvgvelUuwM+wbrmXTV9fgNewSsW95ZUzWoHEksYPuSDGqrGug52Ju4icHYOzZ/MpcQHZXZBf4ZcodvudJ+X0fuc4/StbiKQwFQNVBh5uRnse01Dc7zyKa+poGdy4+xf0MuNnZaBkwIJmqoNxqtTOMUlkECX5gltbGRTe+8g8M3X2CnM7C1r0pwXGdGT30PfJv9Pv+hNL+azQtSyT1ahpu3A0Omh+Mb0bmdKhfCeCTwhVmrLyll9aP3E7xlJ5WdIDmhkbGDYgi76jVw9Wuxn6qqZO09yZZFaVSW1BHSrwsDp4bg7G7fjtUL0b4k8EWHkLtzF4cefwD/nGIyukLhQB3Th8/CcdgjYNvyIG2jTs+eNcfZvSobgH7jAug72h8rGzl0RXQ8Eviiw1BVlZ1ffon+w7dxrWpge0/wjDEwfvQclH43gKblEK8srWPbT+mk7yrCyc2OQdeGEtzXU6Zxig5FAl90OPqqKlY9/Th+v62hQQuJcXqG93Eh6srXIGTEefvmpZSxeUEqJXnV+ER0Zsj0MNy9HdupciHalgS+6LBKU9PZ9vC9hKQcI88NMgfVc22faNzGvgKe4S32M+gNHNqcT9KyTHR1enoN8yH26iDsHKzbsXohWp8EvujwDiz+mdLXX6RLWQ17QsE2ppqJsdOwGvEEdHJrsV9dVQNJyzI5tDkPWwdr4icGEznIG40cuiLMlNkFviy8EpfCUF/P2pfn4v7TIrQGlW0xKjGR9cQOfRAG3AFWLS/CKs6pZPOPqRSkl+Pp78SQ6eF0C5FDV4T5MbvAP0Pu8MWlqMrLZ90j9xG26xAnneDgwAYmBnfCe8xc6H51iwu3VFUlPbmIrT+lU32qnvA4LwZODsXBVQ5dEeZDAl9YpIz1G8l87jF8T5ziqB/UxtUwNSQS27Evg3ffFvs11OvZteoYe9YcR6PVEHtVIH2u8ENrLat1hemTwBcWS9Xr2fz+e9h+9TmOdY1s6wPhPU4xrNckGPkUOHu32Le8uIati9LJ2ncSF097Bl8XRmAvj/YrXohLIIEvLJ6u7BSrHvsPwZsSqbGFnQl6rvKvJij+Hhh0H9g4tNj3+KESNi9I41RhDQFR7gyeFoarlxy6IkyTBL4Qp+Xv2cfeOfcTdOwEx7rAiYG1XOdph8PIp6H3DNA0/9hG32jgwIZcdizPQt9goPcVfvQfFyDTOIXJkcAX4k9UVWXnt9+ge+8t3Cvr2RGp4NXnFGO7haCMfRkCB7fYt6ZCx/YlGRxJLMDW3oqYqwLpNcxXnu8LkyGBL0Qz9NXV/Pbsk/isWIVBA9viFK4IKCIybCyMfh7cQ1rsezK3isSf0zl+qBQndzvirgkmPNYLRebvCyOTwBfiPErTM9n6yL2EHs7khCukDWpkWudTdI69A4Y+BPYtb6ucc6SUbYvTOZlThYefIwOnhuLXveWFXkK0NQl8IS7AoWXLKHr1ebqWVLMvWMGmfyUTHbVoh8+BmFtA2/zzetWgkrqzkKSlmVSW1uHf042EyaF4+Mr+PKL9SeALcYHUhgZ+f/Ul3BcswLrRwOb+WuLCiujf2R/GvAjhY1tcuNXYoOfgxjySVxyjvraR7vFdGTAhGCe3lo9mFKK1mV3gy9YKwtiqCwr5/dH7CN+xn1JH2DtIw1T3AroGDIaxc6Frrxb71lU3sHtVNvvX54ICfa7wpd/YAGw7yYwe0fbMLvDPkDt8YWyZm7eQ/swj+OWXkeKjUD2wgWu1xdj2nQ3D54CLb4t9K0pq2bEsi5QdJ7DtZEXMlTKjR7Q9CXwhLoNqMLDlow+w/vxTnGoa2drbioioUwxV6lFib4Mh/wGHllfgFudUkrg4nZwjZTi52xE/KZiw/jKjR7QNCXwhWoGuvIJVc/5D0Iat1FvDlng7xvnm0V1rDQn/goR7wM65xf7HD5ewbXEGJblVePo7MXBqqBysLlqdBL4QrejEgUMkP34/IWm5nHCF/UMdmemYThdb16a7/djbwLr5g9JVg0rqjhNsX5pJVVk9AVHuJEwOwd1HZvSI1iGBL0QbOLh0KUWvv0i3k1Uc9tNQOsyJ6/VHsXfsBsMegb7XtziVs7FBz4H1eexadXpGT0I34iYE4dhZZvSIyyOBL0QbURsb2fTu29h9+zWOtY1sjbLGM8GeiZVH0bgFw4gnoOeUFvfoqatuYNfKY+zfkIuiKPQZ6dc0o8feqp0/iegoJPCFaGMN5RWsfnoOfr+vQ6+BtQMcGdwH4ktTwatX01bMYWNanMNfcbKWpGWZpO4oxM7BmpjxgUQN9UFrJTN6xMWRwBeinZzKyGTznAcI3Z9KiRNsGtaFGX5lBJdlg188jHwaAge12L/4eCXbFqeTe7QMZw874ieFENq/C0oLfyiE+CsJfCHaWcb6jWS88Dh++aVkdFU4NCaIO+0ycK0ogNBRTcHfrU+zfVVVJedw0x49JXnVdAlomtHjEy4zesTfk8AXwghUg4GdX3xJwyfv4VZZT1K4FfVX9uCWiiSsa8ug5+SmZ/weYc32NxhUUpNOkLSsaUZPYC934ieH4O4tM3pEyyTwhTAifU0Na+e+iOeypWgNBtb0dyBsTBQTslejNNZB9CwY/liLq3YbdXr2r89l16psGuoa6T6wGwOuDsaxsxyuLs4lgS+ECajOK2D94w8SlLSHKntYOdSLawaH0e/okqYGsbfBkAdbXLVbV9VA8spjHNiQi0aj0GeUH/3GBGAjM3rEn0jgC2FCCpJ3s++phwjIKiDXHTZd2YO7I13xPrgYrDv97ard8uKmGT1pOwuxc7QmdnwQPYd4y4weAUjgC2FyVFXlwMKfOPXWK3iWVbM3SEPm5KHc71KBw5HlTYeuDP4PDLi9xVW7RdkVbFucTl7KKVw87YmfFEJIP0+Z0WPhzC7wZXtkYSlUnY5Nb76Fw/xvsdfpWdfHFrsZE7mlYj/azHXg5H3eVbuqqnL8UNOMntL8aryCnBk4JQTvMJnRY6nMLvDPkDt8YSl0pWX8/uRj+G3YhM4Kfh3UmbgZMxib9gvk7oC/WbVrMKikbC8gaVkW1afqCeztQcLkENy6ORjh0whjksAXwkyUHk1l+5z/EHQkgyIXWDk6hJuumULkrs+h8ODfrtpt0OnZvy6H3auyaajXEznImwETgnBwkRk9lkICXwgzk75qNcdfeoZuRac46quwY9IQHkyIx2PLO1CW9berdmurdCSvOMbBjXlotArRo/zpO8YfGzuZ0dPRSeALYYZUvZ4dH3+K+vnHuNTo2NzTirJZM7jfzxPbzW9C5d+v2i0vrmH70kzSk4uwd2qa0dNjiDdarczo6agk8IUwY41VVax/7jm6rPgVFJUVAxzwu/1+ZutPomx9G2rLoMckuOLJFlftFmY1zejJTzuFSxd7EiaFENxXZvR0RBL4QnQAldk5bJ7zEEG791PmAL+M6Mb4O55mcO5WSPwA/mbVrqqqZB8oYdvPGZQVnJ7RMzUU71DX9v8wos1I4AvRgeRtS+Lw04/im1tIlhesuzqae2Y/RuC++ZD8eVOj86zaNegNHN3etEdPTbmOwN4exE8MllO3OggJfCE6GFVVOfDdfCrfewO3ilp2hmlIu248D4/9B86JH8DeeX9atfsvsHM55xoN9Xr2rc1hz5rj6OoaCY/1YsCEIFw8OxnhE4nWIoEvRAdlqKtj86uv47ToB6z1Blb3s8X25ru4K2oYVhtfhsNL/3bVbl11A3tWZ7N/XS4GvUrkYG9irwrEwVWmcpojCXwhOri6wiI2PDEHv63bqLaFZUM6E3f3s1zt4oGy7kXIWAtOZ87avaHZVbvV5fUkrzjG4c35KFqF3sN96Tc2ADvH5s/lFaZJAl8IC3Fy30F2P/EQfunZ5LvBirEh/OOu1+lVWwK/P3dBq3bLi2vZuTyLlB0nsLHVEj3anz4j/WQOv5mQwBfCgqiqSvrSFeS//jxdSirYH6CQPHkID1//Al3ydsO6F06v2o2CK56C8LHNrtotyasiaVkmWftOYu9kTf9xgfQc6o2VtdYIn0pcKAl8ISyQ2tDAjvc+RPPN53Sqa2B9HysqZs/iP+P+jd3RX2H93KZVu/4DYcyL4Nu/2eucyCpn+5JM8lLKcOxsS+zVQXSP74pGFm+ZJAl8ISxYQ9kpNj37LF3WrKbBSuWXBAcCb3uQG/pOQbPnW9jwClQXQ9S1Tat2Owc0e52co6VsX5JJ0bEKXL06EXdNMCF9PVE0snjLlEjgCyGoSMskcc7D+B88TLEzLB3ZjYl3vcRQrx6w9R1IfB9UA8Td1TSH3971nGuoqkrWvpNsX5pJWUE1nv5OxE0Mxr+Hm6zaNRES+EKIP+Ss20T680/S9UQxqd6w7soobr3pOXrbusC6F2Hf/KapnMMehZhbwMrmnGsYDCppO06Q9EsWlSV1eIe5Ej8xmG6yatfoJPCFEGdR9Xr2f/kttR+9g0t1HfsDFXaMieaOm54nQlcNq5+ErE1NM3pGPQeRE5od2NU3Gji8JZ/kFceoqdAR0Mud+InBePg6GeFTCZDAF0K0wFBby473PoYfvsalpp79gQq7xsVy903PE1ScAmueguKj4J9wemC32RyhoV7P/vU57Fl9nPqaRsJiujBgQjCuXrJqt71J4AshzktfU0PiO+9jvfB7nGt07AtU2Dd+IP+64Wl809fC+pegugiipp4e2A1s9jr1NQ3sWX2cfety0DeqRA7sRuz4QBw727XvB7JgEvhCiAuir65my1vvYvfTfJxrG9gbpOHI1UP518xH8No7D7a9B6oe4u48PbDb/Nm5NRU6klce49CmPBRFIWq4D/3HBWDveO54gGhdEvhCiIuir65mwxtv4rhkIc61DewJ1pB2zSjunXIXbtveh73fN83iGfYoxNza7MAuQMXJWnb+mkXK9hNY2WqJHuVP9Eg/bOxl1W5bkcAXQlySxsoq1r3xBq5Lf8KprpHdIVqyJ43jnrFTcV7/MmRthM5BMPo5iLym2YFdgNL8anb8kknGnmLsHKzpf2UAUUN9sLKRVbutzSQCX1GUYOAJwEVV1WsvpI8EvhCmoaGyijWvvozH8qU41enZHWJFwdQJ3J2QgMO6uVB8pOmc3bFzWxzYBSjKrmD70kxyDpfi4GpL7PhAug/sJkcutqLLDnxFUb4ArgaKVFWN+tPr44B3AC3wmaqqr1zAtRZJ4AthnnQVlfz28ot4rfwVpzo9u0KsKZk2iTt6+GG/8bWmgd2eU5oGdt2CWrxOXkoZiUsyKMyqwKWLPXETggnt30VW7baC1gj8oUAV8M2ZwFcURQukAqOBXGAnMJOm8H/5L5e4RVXVotP9JPCFMHP15RWsfOl5fFatxLHeQHKoDVXTp3CbjwbrbR+AobFpYHfoQy0O7KqqyrEDJSQtzaAkrxp3X0fiJwYTEOUuq3YvQ6s80lEUJRBY/qfATwCeVVV17Omf5wCoqvrXsP/rdc4b+Iqi3AHcAeDv798/Ozv7guoTQrS/2rJyfn3pGQJWr2kK/jBbGqZP4kbHIqz2zms6aWvYo01HLrYwsKsaVNKSC0n6JYuK4lq6hboQPzEE7zDX9v0wHURbBf61wDhVVW87/fMNQJyqqve00N8dmEvT/wg++7s/DCB3+EKYi+qSMn598UmC123Aod5Acpgd2mlXM0PdjyZrQ9PA7qhnocfEFgd29XoDR7YWkPxrFtXlOvx7uhE/MQRPf1m1ezHaKvCnAWP/EvgDVFW9t1WqRgJfCHNTVVzKLy/OIWzDZhzqVXaGd8Jpykgm1WxCKToCvgOaBnb9BrR4jUadngMb8tj12zHqqxsJ6deFuGuC6NzVoR0/ifkyqUc6F0MCXwjzVF50kl9eeJTIjYl00qnsjHDAc0Ic406tQakqhB6TYNQzTXv1tKC+tpG9a46zd20O+gYD3RO6Ejs+CCc3WbV7Pm0V+FY0DdqOBPJoGrSdparqoVYoeAIwITQ09Pa0tLTLvZwQwkhKThSy/IVHidq8g046leTujviNjWR46RrQN8CAO5oGdju5tXiNmgodu1dlc2BTLgC9hvrSb1wAnZxl1W5zWmOWznxgOOABFALPqKr6uaIoVwFv0zQz5wtVVee2VtEgd/hCdBRFeQX8+uLD9N66Gzudyp5IJ0KH+xBftr5pYHfowzDgdrCybfEalaV17Pw1i6PbCrCy0dJnpB/Ro/2xlVW7ZzGJhVeXQgJfiI4l73gOq+Y+TN9t+7BtgH2RzkQmONC3Zhe4BjQN7Pac3OLALkDZiWqSlmWRsbsIWwcr+o0NoPdwX1m1e5oEvhDCpGQfy2b1Sw8Rk3gQmwY40MOJ3n319NCng28sjJkL/nHnvUbx8Uq2L83g+KFSHFxsiBkfRI9B3Sz+rF0JfCGEScrISGfty48Qm3QEmwY43MORvlEVhGpONE3hHPXseQd2AfLTyti+JJOCjHI8/BwZPqs7XkHO7fMBTJDZBb4M2gphWY6kHmHjK48RtzMVm0Y4GtmJmIgiAuxrm57tD334vAO7qqqSuaeYzT+mUl2ho9dQH+ImhVjk832zC/wz5A5fCMuy/8gBtr72GPHJmdg0QnqkHTGhefi42Z0e2L3jvAO7utpGti/L5MCGXDo52zDkunBC+nla1FYNEvhCCLOSvH8PO958nIHJx7DWQ1Z3a2JD8/Hy8W6av99zynkHdouyK9jwfQrFxysJiHJn6IxwnD3s2/ETGI8EvhDCLG3bk8Set55i0O4crPSQ211DbGgh7uF9mlbs+se32NegN3BgQx5JyzJRDSqxVwfRZ5Rfh9+KWQJfCGHW1u/czKF3n2XInny0ejgRoRITVoxrzFVNA7vuIS32rSytY/OPqWTtO4m7jwPDZ3ena7BLu9Xe3swu8GXQVgjRnN8S15L2wQsM2VOIRoWT4Y30iyjHZfgtMOyR8w7sZu5tGtStOlVPz8HexE8Kwc7Buh2rbx9mF/hnyB2+EKI5yzavJPvjlxm6txiNCmXhDfTtqcNp/PkHdnV1jexYnsX+tTnYOdkweFooYTFeHWpQVwJfCNEhLVq/hPzPXmfY3lK0KlT2qKP/ACdsrn4RIie0OLBbfLySDd8fpSi7Er8ebgybGY6LZ6d2rr5tSOALITq0735fSMWnrzNsfyW19mDbp4qoIb1QrnwJvKOb7WMwqBzcmMf2pRkY9CoxVwXSd7Q/WivzHtSVwBdCdHh6g573vn0Dv2+/o0duIyUeKv7Rp/AfObnpjF3nbs32qyqrZ8vCVDJ2F9O5mwPDZ0WY9WlbEvhCCItRXFXOx//3AEN/TaRLBRQF64nuU43L2Ptg4L1g0/yjm2MHTrJpfiqVpXVEDurGwCmhZjmoa3aBL7N0hBCXa1fmIda+8gCjE3PQGqAmqo6YvnZYXfkM9LoONOc+ummo17Pz1yz2/p6DnYMVg6aGEh7X1awGdc0u8M+QO3whxOX6cd1iKt5/icGHq6lwAMfoCnoMCEUZ9woEJDTb52RuFRu+P0phVgU+EZ0ZPisCVy/zGNSVwBdCWLSGxkbe/XwuYT8sJKxAT7EXhEWX0G3QOBj1HLgFndNHNagc2pJP4s8ZNDbo6T8ukP5jA9Bam/agrgS+EEIABeUlfPnqvxm2ehduVVAUpicmqhyHEaePWrQ7dwVudXk9WxemkZZchKtXJ4bPisAnorMRqr8wEvhCCPEn21P2sO3VBxmZVICqQH3vWvpHabEa/Tj0uwm0526rfPxQCRvnp1Bxso7u8V0ZeG0o9o6md66uBL4QQjTjuxXzaPz4DeJSaznlDJ2jy4no7Ycybi6EjjqnfYNOT/KKY+xdfRxrey2DpobSPaGbSQ3qSuALIUQL6hsbeO+jZ+i5YCmBxQYKvSGyz0m8+g9tOmqxS/dz+pTkV7Hx+xQKMsrxDnNl2KwI3Lo5GKH6c5ld4Mu0TCFEe8spLeT7l+5lxNoDONZCcaSeAd1L6DT4HzB8Djh4nNVeNagc2VbAtsXpNNTr6Tc2gP5XBmBlbdzD1M0u8M+QO3whRHvbfCCJPa89zIhdxTRYQWN0LTHhejQjHoa4O8/ZmK2mQsfWn9JITSrExdOeYbMi8ItsedfOtiaBL4QQF0FVVb5e8gU2/32Pvpn1lHRW8IouIzTCE2XMCxB5zTkbs+UcKWXjvBTKi2sJH+DFoGvD6OTc/oO6EvhCCHEJanX1fPDuHKIXr8KnVCU/QKFPVBEeUQOaTtzy7ntW+8YGPbtWZbN7VTbWtloSJofQY5A3iqb9BnUl8IUQ4jJkFuex6MV7uGLDUex0UBxlYGBYIbaxM05vzOZ9VvvSgmo2zkshP+0U3UJcGDYrAncfx3apVQJfCCFawdpdG0l543GG7S2lxg6U6Fr6BdejGfLv0xuz/W+mjqqqHE08wbaf0tHVNhI92p+Y8YFY27TtoK4EvhBCtBKDwcAXCz7E9fP/0jNHR5GHBr/oEoKDOjfd7feeftbGbLVVOrb9lM7RxBM4e9gxdGYEAT3d26w+CXwhhGhllXW1fPJ/DxG7bD1dylVyQzTE9DhB59AoGPsSBA46q31eShkb5qVwqrCG0JguDJ4WhoNL80cxXg6zC3yZhy+EMBcpeVksf/FertiSgdYAJb1VBgefwLr31TD6eXAL/qOtvsHA7tXZ7FqZjdZKIWFyCD2H+LTqoK7ZBf4ZcocvhDAXq7b9Rs6bTzP4YAXlDgo2/Wrp61eFEn8nDHkI7F3/aHuqsIYN81LISynDK8iZ4bMj8PB1apU6JPCFEKIdGAwG/vvtm3T76mvCChop8NISEl1MgI9D02rd/jf/sTGbqqqk7ihk66I06qob6TPSjwFXB2Fte3mDuhL4QgjRjk7VVPLZqw+Q8OtW3KogO8KKgd3zcPYLbZq/Hzb6j7Z11Q0kLk7n8NYCHN1sGTYjgsDeHue5+vlJ4AshhBEcyjrKmrn/ZkTicVQNlEVrGBKQh1XEFU3B3yXyj7b5aafYMC+FsoJqgvt6MmJ2d+wcL/5M3fMFvmkf3SKEEGasZ1B37v/sN3I+eIW9wY503Wkg6Tcf9m7bh/rhQFj+AFSfBMA7zJXpT8QSNzGY8qIarGxaP57lDl8IIdqB3qDnk89eJui7Hwgs0pPjY01k7yL8ulg3nbYVd9cfG7MZ9AY02ksLfLnDF0III9NqtPzzjicZ8OsWVkyOxbmsgYqVnVmZ7E7VL8/B+7FwaAmo6iWH/d+RwBdCiHbk7uTKgy9/g90PC1gf543vkXpSV3ixfp+K/oeb4MsrIW93m7y3BL4QQhhBdHgv7vl6LWlvPsNRX3u6JjayZY0f+/dmoH42GiryW/09JfCFEMKIpoybwbRlO1j/z6noNCrWq61YnhhEdnF1q7+XSQa+oigTFEX5tLy83NilCCFEm7O2suKf971I7xUbWTW+D/oGHc6enq3+PjJLRwghTExNfR2dbO0uqa/M0hFCCDNyqWH/dyTwhRDCQkjgCyGEhZDAF0IICyGBL4QQFkICXwghLIQEvhBCWAgJfCGEsBAmvfBKUZRiILuFX7sA51uK6wGcbPWijOPvPqu5vGdrXPNSrnExfS607YW0O18b+X6a5vte7jUvtX9rfkcDVFVtfpmuqqpm+Q/49G9+n2zsGtvrs5rLe7bGNS/lGhfT50LbXki787WR76dpvu/lXvNS+7fFd7S5f+b8SOcXYxfQjozxWdviPVvjmpdyjYvpc6FtL6SdpXxHjfU5TfE7eqn92+I7eg6TfqRzORRFSVZb2E9CCGOT76cwBnO+w/87nxq7ACHOQ76fot112Dt8IYQQZ+vId/hCCCH+RAJfCCEshAS+EEJYCIsIfEVRHBRF+VpRlP8qijLb2PUI8VeKogQrivK5oiiLjF2L6LjMNvAVRflCUZQiRVEO/uX1cYqipCiKkq4oymOnX54CLFJV9XbgmnYvVliki/mOqqqaqarqrcapVFgKsw184Ctg3J9fUBRFC3wAXAn0AGYqitID8AVyTjfTt2ONwrJ9xYV/R4Voc2Yb+KqqbgJK//LyACD99N2SDvgBmAjk0hT6YMafWZiXi/yOCtHmOlr4+fC/O3loCnofYDEwVVGUj7Cc5e7CNDX7HVUUxV1RlI+BvoqizDFOaaKjszJ2Aa1MaeY1VVXVauDm9i5GiGa09B0tAe5q72KEZelod/i5gN+ffvYF8o1UixDNke+oMJqOFvg7gTBFUYIURbEBZgDLjFyTEH8m31FhNGYb+IqizAcSgQhFUXIVRblVVdVG4B7gN+AIsEBV1UPGrFNYLvmOClMjm6cJIYSFMNs7fCGEEBdHAl8IISyEBL4QQlgICXwhhLAQEvhCCGEhJPCFEMJCSOALIYSFkMAXQggLIYEvhBAW4v8BNM4ZNaQqDaEAAAAASUVORK5CYII=\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.loglog(data1[:,0], data1[:,1])\n", + "plt.loglog(data2[:,0], data2[:,1])\n", + "plt.loglog(data3[:,0], data3[:,1])\n", + "plt.loglog(data4[:,0], data4[:,1])\n", + "plt.loglog(data5[:,0], data5[:,1])\n", + "\n", + "del data1,data2,data3,data4,data5" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 39, + "outputs": [], + "source": [ + "data1 = np.loadtxt(\"../data/1Euler_oaf_harmonic.dat\")\n", + "data2 = np.loadtxt(\"../data/2Heun_oaf_harmonic.dat\")\n", + "#data3 = np.loadtxt(\"../data/3Exact_oaf_harmonic.dat\")\n", + "#data4 = np.loadtxt(\"../data/4BDAS_oaf_harmonic.dat\")\n", + "data5 = np.loadtxt(\"../data/5MBD_oaf_harmonic.dat\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 40, + "outputs": [ + { + "data": { + "text/plain": "" + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5BklEQVR4nO3dd1AV1/vH8ffSe1NsoIBiQ0EUbNhiL7HFXiK2SKwxlhhNYmKqPbHFrrEbwRo1Yo29gr2jCIgVRUGatP39cb/Jzxgb9V7gec04E7hnzzk7c+eT5ezZZxVVVRFCCJH/6Wl7AkIIIXKHBL4QQhQQEvhCCFFASOALIUQBIYEvhBAFhAS+EEIUEAbansCbFC5cWHV2dtb2NIQQIs8IDg5+pKqq/as+0+nAd3Z2JigoSNvTEEKIPENRlPDXfSZLOkIIUUBI4AshRAEhgS+EEAWETq/hCyHEu0hJSSEyMpKkpCRtTyXXmJiY4OjoiKGh4TsfI4EvhMjzIiMjsbS0xNnZGUVRtD2dHKeqKo8fPyYyMhIXF5d3Pk6WdIQQeV5SUhKFChUqEGEPoCgKhQoVyvBfNPky8C+f/Y3Y++e0PQ0hRC4qKGH/t8ycb75b0klPS+XT4GlE60ED1ZwPXFtTu/pw9E2ttT01IUQ+pq+vj7u7+z8/d+vWjbFjx762/bJlywgKCmLOnDm5MT0gHwa+ih4uesOxig3gmOU9dt0KwD7kd943c6J9lf6UqdgRCtiVgBAi55mamnL27Nkc6z81NRUDg6xFtk4u6SiK0kZRlIUxMTEZPlZfX48Fff2Y4buFdsXX4hrVFJskS1Y+v037U9/SbakHa7d9REzU1RyYuRBC/JuzszOPHj0CICgoiPfee+8/baKioujYsSPVq1enevXqHDlyBIAJEybg5+dHs2bN8PX1zfJcdPIKX1XVrcBWb2/vAZntw9HWjM9bupPSbBq7Lz9gzdFTJMWuItb6Gj89PsHU7Z1ooFjSvkxbfGp8gqGReTaegRBCW77deonLd2OztU+3ElZ806bSG9skJibi6en5z8/jxo2ja9eu79T/8OHDGTFiBHXr1iUiIoLmzZtz5coVAIKDgzl8+DCmpqaZnv/fdDLws5Ohvh6t3IvTyr0tN6MasfZEBAfOb6ewWSBBVg/Yc3MNdtdX8b6FC+2q+FG+XBtZ8hFCZFhWlnT27NnD5cuX//k5NjaWZ8+eAdC2bdtsCXsoAIH/ojL2FnzV2o3Rzcuz/Xx31hwPwfbpesxsTrFWL4yVx7+kwtHxtCtam1a1RmNn56rtKQshMuhtV+K5zcDAgPT0dIDXbqNMT0/n2LFjrwx2c/PsW33QyTX8nGZiqE9HL0c2DGnIpL4TKVNiCcZhY3B/UI6UFJXJUUdo/Ed7hq2sw95jU0lJTtT2lIUQeZSzszPBwcEAbNiw4ZVtmjVr9q/dOjl187dABv6L3EpY8eMH7uwZ14NW9X8mKWEO9qHdqBpjz8WUp3x6fQWNVlfnp/XtuHRjB6qqanvKQggd9Pca/t///t6S+c033zB8+HDq1auHvr7+K4+dNWsWQUFBeHh44Obmxvz583NkjoouB5i3t7ea2/XwVVXlzO2nrDoezs7zYbiZBGJVKJizpokkKwquqgHtivvwfo1R2NuWztW5CSFe7cqVK1SsWFHb08h1rzpvRVGCVVX1flX7ArWG/y4URaFaKVuqlbLlaWs31gd7sOZEBMrtm9QutJ0YqxCm3z/IL1sO4GNgS7tyHWhYdSDGhtlzU0UIIXKKBP4b2JgZ8VG90vSv68Kxm5VZfaIqwZfuUc4wiApFD3De+D6fXVmK5aWltLQqS9uqA/FwaVrgHvEWQuQNEvjvQFEUfFwL4+NamIexSfgHlWPtyXo8fhpNA7s9GNsE80fsNfwPjcL5kCHtiteldc2RFLN21vbUhRDiH7KGn0lp6Sr7rz1k9YkI/rr2EAe9SBo57CbM8CpnjBQUVaWWYSHalu9E4yr9MTU00/aUhci3ZA3//8kafg7Q11NoXLEojSsWJfJJAmtPRrDuVBkexyXS2PYspez+4mjafcZdWoj5xYU0ty5P26oDqebUWJZ8hBBaIYGfDRxtzfiseQWGNy7H7ssPWHXcnqU3q2Gjn0DPkgd5anCMwKdX2HhgBCUxoo1DfdpW/xQHaydtT10IUYBI4GcjIwM93vcozvsexbnxMI61JyNYG2xNTGIL6tvdp6b9Ts6kXmRe5G7m3tlDdcNCtC3fhWYefTCTJR8h8jQLCwvi4uL++Vkb5Y/fRgI/h7gWsWB8azc+a16ebefvsfqEDVOvFcPcUGWAy2X09ALZmXCP8Rfn8dOFeTS1qUBbz4+p7tQIPaXAPw8nhMgBEvg5zMRQn05ejnTycuTinRjWnIxg4RkDEpIrUbt4Ol8V2ceV+APsfHKZPw6MoLhiRBuH92jnNZRSNu/+rkohhO6Kiopi4MCBREREADBjxgzq1KnDhAkTsLCwYPTo0QBUrlyZbdu2AdCyZUvq1q3L0aNHcXBwYMuWLVkuoqaTga8oShugjatr/ipeVtnBmp8+cGdcywpsPnuX1cfD+fxcEyyMm+NX7ikl9DaxJ/Y8i2/vZGHkLqoaFaZtha60quwrSz5CvKsdY+H+hezts5g7tJz0xiYvl0eOjo6mbdu2wJvLH79OSEgIa9euZdGiRXTp0oUNGzbw4YcfZuk0dDLws6Mevi6zNDGkVy0nPqxZitMRT1l9PJw5F1SSU32pXcqcKXZnuf1kM1vj7vHt+V+ZfX4BfSr0oGu1IRL8Quiol8sj/72GD28uf/w6Li4u//wPxMvLi7CwsCzPUScDv6BQFAUvJ1u8nGwZ39qNDacjWX0igo+PlsXW7Av6uhsxigDWPNzDz1dX8Nu1tfhW7EV3Tz/MDeWFLUK80luuxLXhdeWPXyydDP8un2xsbPzPf+vr65OYmPWqvXJ3UEfYmmvKOOwd2YDVH9WkdplCzDqVgO+JVtiZzWWhmTeV458x8/JSmq+tz4LgmTxLfvMVghBCN7yu/LGzszOnT58G4PTp09y6dStH5yGBr2P09BTquBZmbk8vjo5txNCGrmy6noLv2c6UL/wry40qUjXuKXMuLqb5uveYGzyTmOcZf/evECL3vK78cceOHYmOjsbT05N58+ZRrly5HJ2HlFbIAx7GJjFzbwi/n7qNiYEeX1RLpUrMXBbHX2OvuRkWekb0cPOlV6Xe2JjYaHu6QuQ6Ka3w/95UWkGu8POAIlYm/PiBO7tH1KdBeXu+PK7Q5/ZIWjlNJiDZFp/YJyy6sJjmAY2ZETyD6KRobU9ZCKGDJPDzkNL2Fszt6cWmwT6UtjfH75ApA59MoEfpb9gQb0iD2CcsvbCEFgHN+DnoZx4lPtL2lIUQOkQCPw+qWsqWdX61+K1PdYwNDeh4oDBj0n/h4zKj2fw0lcax0Sy/9Bst1zdnyqkpRCVEaXvKQggdIIGfRymKQsMKRfhzeD2mda7Cg7hUGu93YqLlQj4t/TF/PHxG85gnrLm8khYbmjPxxEQexD/Q9rSFEFokgZ/H6espdPJyZN/o9/iyVUVORiZS62BlFjqsZIxTd7befUzrmBj8r66l5caW/HD8B+7F3dP2tIUQWiCBn0+YGOozoH5pDn7WkI/rl2H95Ti8j9dlY/k1fFGiJdsi79Hu2TM2XAug1aZWfHvsW+7E3dH2tIUQuUgCP5+xNjNkbMsK/DX6PdpXLcHMU/F4n2vDfs81fGVbkz8jbtMxLpEtIRtpvfF9vjn6Dbef3db2tIXI8xRFoVevXv/8nJqair29Pa1btwY0pRbs7e3x9PSkUqVKdOrUiYSEBAAmTJiAg4MDnp6elC1blg4dOvyrFEN2kcDPp0rYmDKlUxUCP61PzdKF+PJQMj63+nCh5kq+MHVlR/htuiSmse3GFtpsasOXh78kPDZc29MWIs8yNzfn4sWL/5RA2L17Nw4ODv9q07VrV86ePculS5cwMjJi3bp1/3w2YsQIzp49S0hICF27dqVRo0ZERWXvhgsJ/HyuXFFLFvf2xv/j2jjYmOK3D5o9Hk2kzwLGqtYEhofTI1mfXbd20HZzW8YdGkdoTKi2py1EntSyZUu2b98OwNq1a+nevfsr26WmphIfH4+tre0rP+/atSvNmjVjzZo12To/KZ5WQNRwsWPDIB92XX7AlMCrdNxthnepSUypc4MxF36h390IVpSqyO9hu9geup0Wzi3w8/DD1TZ/lagW+d/kk5O5Gn01W/usYFeBz2t8/tZ23bp147vvvqN169acP3+efv36cejQoX8+X7duHYcPH+bevXuUK1eONm3avLavatWqcfVq9p6HXOEXIIqi0LxSMXZ+Wp+JHdyJeJJEo132fGw9n/Q63zEy6iGBt0Lpp1+YA7f/osMfHRi1fxTXoq9pe+pC5AkeHh6EhYWxdu1aWrVq9Z/P/17SuX//Pu7u7kydOvW1feVE2Ru5wi+ADPT16F6jFO09HVh65Bbz99+k1rXS9PRcxRibvXwaPI8+aYmsLFuLNXcOsyt8F41LNWZglYFUsKug7ekL8UbvciWek9q2bcvo0aPZv38/jx8/fmUbRVFo06YNs2fPZuzYsa9sc+bMGby9X1kSJ9PkCr8AMzXSZ0hDVw6OaUjfOi6sO/8E70PVmFk5ANOq/RkWcpLAW2EMsijPyXsn6Ly1M8P2DePSo0vanroQOqtfv358/fXXuLu7v7Hd4cOHKVOmzCs/27BhA7t27XrtPYDM0snAVxSljaIoC2NipOxvbrA1N2J8azf2jmrA++7FmXEsmhrBTVhbYxMWFdsz+MIedkbcYahtVU4/OE237d0YvGcw56POa3vqQugcR0dHhg8f/srP1q1bh6enJx4eHpw5c4bx48f/89kvv/zyz7bMVatWsW/fPuzt7bN1blIeWfzH5buxTA68yoHrUZSwNmFCLYWmd+ejhOwkzqo4v1dsyPKYizx9/hSfEj4MqjIIzyKe2p62KMCkPPL/k/LIIkPcSlixvF8N1nxUk8KWxvjtTKRl1FBON16LuXUpPjqxhp1R8Yx0bM7V6Cv02tGLj3Z9RNB9+Z+zELpMAl+8lo9rYbYMqcOcHlVJTEmjw3aVbikTuNlkEWZ6RvQ9tIgdMXqMdvmAG09u0HdnX/oG9uXkvZM5ssNACJE1EvjijRRFobVHCfaMbMD37Spx81E8jbeZM8R6Fg8b/YxZXBS9980k8Lk1Y8t/SERsBP139adPYB+O3j0qwS+EDpHAF+/EUF+PXrWdOfBZQz5tUpa/rkdTO7A43zgt51mDbzG5d46egT/xJ458UdmPO3F3+Hj3x/Ta0YvDdw5L8IscV9C+Y5k5X7lpKzIl6tlzZu8LYc2JCAz19RhS254BBtswPjUf0pJJrubL5pKVWRyyjnvx96hcqDKDPAdRz6EeiqJoe/oin7l16xaWlpYUKlSoQHy/VFXl8ePHPHv2DBcXl3999qabthL4IkvCHsUzbdc1tp2/h525EWPqWNM5fi36Z1aAvhEptQaytbgrC6+s5E7cHWoVr8XYGmMpY/Pq/cdCZEZKSgqRkZEkJSVpeyq5xsTEBEdHRwwNDf/1ewl8kePORz5l0o6rHL35mJJ2pnzjY0Lje4tQLm0EUztS6o0gwMqGXy8sID4lnu4VujPIcxBWRlbanroQ+YoEvsgVqqpyMOQRk3Zc5cq9WCqVsOKHGqlUvT4TQv8C65I8afg5cxJvEXB9PTbGNnxS7RM+cP0AfT19bU9fiHxBAl/kqvR0lT/O3WXarmtEPkmkXtnCfOf+CJfTE+H+BXCqw1WfgUwM3cDph6epaFeRcTXHUbVIVW1PXYg8TwJfaMXz1DRWHY9gzr4QniSk0L5KUb4vdQbLIxMh8QlqVV8Cy9Vh2oUFPEx4SOvSrRnhNYIiZkW0PXUh8iwJfKFVsUkpLDhwk0WHbmFpbMDk90vR5OEyOLkQjMxJqD+axcbpLL+8En09ffw8/PB188VI30jbUxciz5HAFzrh+oNnjPQ/y8U7sXSo6sC3PoZY7v8Kbu6DwuW5/d4opj08yr7b+yhpWZIx1cfQwLFBgdhmJ0R2kVo6QieUK2rJpsF1GN64LFvO3aXZqvscqrkAuv8OacmUXO/HzAcPWVBzAoZ6hgzbN4xBewdxK+aWtqcuRL4gV/hCK85HPmWk/zluPIyjVy0nxjVzwez0Ijg4FdKSSak5kN+LOTH34hKSUpP40O1DPvb4GAsjC21PXQidJks6QiclpaQxbec1lhy5hZOdGdO7VMHLLhn2fgdnV4NFUR43+IxZybfZdGMzdiZ2fOr1KW3LtEVPkT9OhXgVCXyh046HPmZ0wDnuPk1kQP3SjGxaDuP7Z2HHGLgTBA5eXKoziIlhWzgXdQ73wu6MqzEOd/s3v1FIiIJI1vCFTqtVuhCBn9ana/WSLDgQStvZR7ik5wr9d8MHCyDmDpX8P2JFkjk/eY/hfvx9evzZg68Of8WjxEfanr4QeYZc4Qud8tfVh3y+4TzR8ckMb1yWQe+VwSA1AQ5Nh2NzQN+I+DqfsNBMnxVXV2Osb8xAj4H0rNgTQ33Dtw8gRD4nSzoiT3makMz4LZfYeu4uVUraML1zFVyLWEB0KOwaD1e3ga0L4Q1GMOXRSQ7eOYizlTOf1/icug51tT19IbRKlnREnmJjZsTs7lWZ06Mq4Y/jeX/WIZYevkW6jQt0Ww29NoOBMU6bP+HXu3f41WscKiqD9gxi2N5hRMRGaPsUhNBJcoUvdNrD2CTGbrzAvqsPqV26EFM7e+BoawZpqRC0BP76EZ7HkVK9P6uKuzD/8nJS0lPwdfPFz8MPM0MzbZ+CELlKlnREnqaqKgFBkXy79RKKovB1azc6eztqnsCNf6wJ/eDfwMSGqHojmJF2nz9Ct1LEtAifen1K69Kt5WldUWBI4It84XZ0Ap+tP8fx0GgaVyjCxI7uFLE00Xx4/wLsGAvhh6GoO+fq+DEx4k8uPb6Ep70nY2uOpVKhSto9ASFygQS+yDfS01WWHQ1jcuBVTI30+aF9ZVp7lNB8qKpweYvmxm5MBOkV27GlQgNmXF3Bk6QndCjbgU+qfYKdiZ12T0KIHCSBL/KdGw/jGOV/lnORMbSpUoLv21XCxux/1TVTEuHobDj0M6DyrPYg5psbs+a6P6YGpgz2HEzXCl0x1JNtnCL/kcAX+VJqWjrz9t9k5t4Q7MyNmNzJg4blX6ilHxMJu7+Bi+vBypHQ+sOYHH2Go/eOUsa6DJ/X+JzaJWpr7wSEyAF5LvAVRWkDtHF1dR0QEhKi7ekIHXfxTgyj/M9x7cEzutcoyZfvu2FhbPD/DcKPaco03D+PWqo2+726MOVGAJFxkTQu1ZjR3qNxtHTU3gkIkY3yXOD/Ta7wxbt6nprGz7uvs/BgKI62pkzrVIWapQv9f4P0NDizSlOYLeExz6v2ZKVjORZeXU1aehp9K/elX+V+so1T5HkS+KLACAqLZlTAOSKiE+hfx4XRzctjYvjCC9ITn2pKMJ+YD4bm3K8zhF/Ux/wZtoOiZkUZ7T2a5s7NZRunyLMk8EWBEv88lYk7rrDqeASuRSz4uUsVPBxt/t0o6jrsHAc39kDhcpyuPYBJd/dwJfoKXkW9GFdjHOXtymtl/kJkhQS+KJAOXo9izPrzRMU9Z2hDV4Y2csVQ/6VqItd3QuA4iL5JWtnmbKz4HrOuryU2OZbO5Toz1HMoNiY2Wpm/EJkhgS8KrJiEFL7deomNZ+7g7mDNz12qULao5b8bpSZrlngOTIHUJGJq9GeelTm/39iIuaE5w6oOo1O5ThjoGbx6ECF0iAS+KPACL97ji00XiXueymfNytOvrgv6ei+t0z97APu+gzOrwdyekLqDmRx7kRP3T1LWtizjaoyjerHq2jkBId6RBL4QQNSz53yx6QK7Lz+ghrMd0zpXoVShV+zKuXMadnwOkSdRS1Rlr1c3pt7axN34uzRzasZo79EUtyie+ycgxDuQwBfif1RVZePpO0z44xJpqsqX71ekR41S/92Vo6pwIQB2fw3P7pHk3onfSlZg6XV/APpV7kfvSr1lG6fQORL4QrzkztNExqw/x5Ebj2lQzp7JHT0oZm3y34bP4+DwL5pSDXoG3Kvtx3Qlhp0ReyhiWoRh1YbJS9WFTpHAF+IV0tNVVp0I56c/r2BsoM937SrRtkqJV+/Bj74Fu8fDla1g68wZn4+Z+vAwFx5doIJdBUZ7j6Zm8Zq5fxJCvEQCX4g3uPUonlH+Zzkd8ZRW7sX4ob07duZGr24cul+zjfPhZdLdPiDQvSUzLi3hXvw93nN8j5HeI3GxdsnV+QvxIgl8Id4iLV1lwcGb/LL7OtamRkzq4E4Tt6KvaZwCR2bCgclgbElS8x9ZpcSz+OJinqc+p3P5zgyqMghbE9vcPQkhkMAX4p1duRfLSP9zXLkXS2cvR75u44alyWvKKD+8CluGwJ0gKP8+j5t8xbybG1l/fT1mBmb4efjRo2IPjPRf89eCEDlAAl+IDEhOTWfm3uvM23+T4tamTO3sgU+Zwq9unJ4Gx+fBvu/BwBiaT+Smc02mB//MoTuHcLBwYITXCJo5NZP6PCJXSOALkQmnI54w2v8coY/i6ePjzOctKmBqpP/qxo9vwh/DIPwIlGkMbWZwNCGSaUHTCHkSgqe9J6Orj6aKfZXcPQlR4EjgC5FJiclpTA68yrKjYZQubM70LlWoWuo1a/Pp6RC0BPZM0Pzc9FvSqvZmy62tzD4zm0eJj2jh3IJPvT7FwcIh185BFCwS+EJk0ZEbj/gs4Bz3Y5MY/J4rnzQui5HBa/beP42ArcPh5j5wqgttZ5FgVZylF5ey/NJy0tV0err1ZID7ACyNLF/dhxCZJIEvRDaITUrh+62XCQiOxMPRmkW+3hS1esXDWqB5Uvfsagj8AtKSofF4qDmQ+4lRzD4zm603t2JjbMNgz8FSmE1kKwl8IbJR4MX7jPI/i6WJIUv6eFOphPXrG8feg20j4PoOcKwO7X4F+/JcfnyZaUHTOHX/FC7WLoz2Hk09h3pyY1dk2ZsCX54HFyKDWlQuRsBAHxQFOs8/xr6rD17f2Ko4dF8LHZdobuzOrwsHp+FmU5YlzZYwq+EsVFVlyN4hDNg9gGvR13LvRESBI1f4QmTSg9gk+i8/xeW7sYxv7UbfOm95wjYuCnZ8Bpc2QTEPzdV+cQ9S0lPwv+bP/HPziXkeQ3vX9gytOpQiZkVy50REviJLOkLkkITkVIb/fpbdlx/Qu7YT41u7YfDyW7VedmUrbBsJidFQdyTUHw0GxsQmx7Lo/CJWX1mNgZ4BfSv3pbebVOQUGSOBL0QOSktXmRx4lYUHQ3mvvD2zu1d9/dO5f0uIhp1fwrk1YF9Rc7Xv6AXA7We3mRE8g13hu6Qip8gwCXwhcsGaExGM33KRskUsWNKnOg42pm8/KGS3Zgvns3tQewg0/BIMNcedeXiGaaemcf7ReanIKd6ZBL4QueRQSBSDV53GxEifJb298XC0eftBSbGaF60E/wZ2ZaDdHHDyATQvbAkMC2RG8Azuxt+VipzirSTwhchFIQ+e0XfZKR7FPWdG16q0qFzs3Q4MPaApz/A0HGr4QeNvwNgCgOdpz1l1eRWLLywmKTVJKnKK15LAFyKXRT17jt/KIM7efsrYFhXwq1/63fbYJ8fD3u/hxHywKQltZkGZhv98HJ0Uzdyzc6Uip3gtCXwhtCApJY1RAefYfv4e3aqX5Pv2lTF82w6ev0Uchy1D4XEIVPOFZj+Ayf8/4BX6NJTpwdM5GHlQKnKKf5HAF0JL0tNVft59nTl/3aCOayHm9vTC2vQtO3j+lpIEByZpXrZiURRaz4DyLf7V5NjdY0wLmsb1J9elIqcAJPCF0Lr1wZGM23ieUnZm/NanBqUKZWBv/Z3Tmqv9h5fAvQu0nAxmdv98nJaexh83/2DWmVlSkVNI4AuhC47dfMzAVcEY6Cks9PXCy8nu7Qf9LTUZDv8MB6eCqS20mgaV2v+rSUJKAr9d+o1lF5dJRc4CTAJfCB0RGhVHv2WnuBuTxLTOVWhbpUTGOrh/UfNaxXtnoWJbTfBb/vvduw/iHzD7zGz+uPmHVOQsgCTwhdAh0fHJDFwZzMmwaEY1LcfQRq4Zu9malgrHZsNfE8HIDFpMBo8u8FIfVx5fYVrQNE7ePykVOQsQCXwhdMzz1DTGbbjAxjN36FDVgYkd3TE2eM3rE18n6jr8MRRun4CyzTQ3da3/vW6vqir7b+/n5+CfCYsNo2bxmnzm/Rnl7cpn27kI3SKBL4QOUlWVOftuMH33dWq42LHgQy9szTO4nz49DU4ugr3fgp4BNPseqvX+z9V+SnoKAdcCmHdunlTkzOck8IXQYVvO3uGz9ecpYW3C0j7VKW1vkfFOom9pntINOwQuDaDtLLB1/k+z/1TkrNSXPpX7YGrwDnV/RJ4ggS+EjgsOj2bAimDSVZX5H3pRq3ShjHeSng6nl8Ou8aCmQZMJUH0A6P33Ya8XK3K6WLswvcF0ytqWzfqJCK2TN14JoeO8nOzYPLgOhcyN6LXkBOuDIzPeiZ4eePeFIcfBqQ7sGAO/tYRHIf9pWtKyJNPfm86iZot4lvyMHtt7sPnG5qyfiNBpEvhC6IhShczYOLgONVzsGB1wjmk7r5Genom/wK0doWcAtJ8PUVdgXh04PEOzu+cltYrXIqBNAFXsqzD+yHi+PPwlCSkJWT8ZoZMk8IXQIdamhizrW4Ou3iWZ89cNPvn9DEkpaRnvSFHAszsMOQllm8Keb2BJE3hw6T9NC5sWZkHTBQyqMoitN7fSY3sPbj69mQ1nI3SNBL4QOsZQX49JHd0Z17IC287fo/ui4zyKe565ziyLQddV0Ok3eHobFjSA/ZM0T+6+QF9Pn8Geg1nQdAFPnj+h+/bubL25NRvORugSCXwhdJCiKHzcoAzzP6zGlXuxtP/1CCEPnmW2M6jcAYacALd2sH8iLGoId8/8p2ntErVZ32Y9lQtX5ovDX/D1ka9JTE3M4tkIXZFrga8oSmlFUZYoirI+t8YUIq9rUbk46/xq8zw1nQ5zj3IoJCrznZkXhk5LoNtaiH8EixrDngmaqpwvsDezZ2HThfh5+LH5xmZ6bO9BaExo1k5E6IR3CnxFUZYqivJQUZSLL/2+haIo1xRFuaEoytg39aGqaqiqqv2zMlkhCqIqJW3YPKQODram9PntFGtPRmStwwqtNFf7nt3h8C+w7kPNA1wvMNAzYFjVYcxvMp/HiY/ptq0b20O3Z21coXXveoW/DPhXIW5FUfSBX4GWgBvQXVEUN0VR3BVF2fbSP3mcT4gscLAxJWBgbeqVLcy4jRf46c8rmdvB8zdTG2j3q6b42o3dsO/7VzbzcfAhoE0AFe0qMvbQWCYcnUBSatIr2wrd906Br6rqQSD6pV/XAG7878o9GfgdaKeq6gVVVVu/9O9hNs9biALH0sSQxb7e+NZ2YuHBUAauCiYh+b9bLTOkxgDw6qO50r+44ZVNipoXZUnzJXzk/hEbQjbw4Z8fEhYTlrVxhVZkZQ3fAbj9ws+R//vdKymKUkhRlPlAVUVRxr2hnZ+iKEGKogRFRWVhvVKIfMhAX4/v2lXmmzZu7LnygK4LjvMgNotX3C2nQslasHkI3Dv/6nH1DBhebThzG8/lQcIDum7rSuCtwKyNK3JdVgL/VTVWX/s3pqqqj1VVHaiqahlVVSe+od1CVVW9VVX1tre3z8L0hMi/+tZxYZGvNzej4mj/6xEu343NfGcGRtBlhebFKr/30NzQfY16jvUIaBNAebvyfHbwM74/9j3P0zK5ZVTkuqwEfiRQ8oWfHYG7WZuOEOJdNa5YlICBtVFV6Dz/KPuuPsh8Z5ZFodtqiI8C/96QlvLapsXMi7Gk+RL6Vu6L/3V/ev3Zi4jYLN5IFrkiK4F/CiirKIqLoihGQDfgj+yZlhDiXVQqYc2WoXVwsTfno+VBLDtyK/OdOVSDNrMg/DDs/OKNTQ31DBnpNZI5jeZwN/4uXbZ1YWfYzsyPLXLFu27LXAscA8orihKpKEp/VVVTgaHATuAK4K+q6n+f2xZC5KiiVib4f1ybRhWKMmHrZb7ZcpHUtPTMdValK9QeCicXwumVb23eoGQDAloH4GrjyugDo/nx+I8kpyW/9TihHVIeWYh8Ii1dZeKfV1h8+BYNy9szu0c1LIwz8R7btFRY3QnCj0CfP6Fk9bcekpKewszgmSy/vBy3Qm5MazCNkpYl33qcyH5SHlmIAkBfT+Gr1m780L4yB0Me0WneUe4+zURZBH0D6LQUrEpoHsqKvffWQwz1DBldfTSzGs7i9rPbdN3alT3hezJxFiIn6WTgK4rSRlGUhTExMdqeihB5zoe1nPitT3XuPEmk3a9HOB/5NOOdmNlpSjA8f6YJ/ZR32/rZsFRDAtoE4GztzIj9I5h0chIpb7gBLHKXTga+qqpbVVX1s7a21vZUhMiT6pezZ8NgH4z09eiy4Bg7L93PeCdF3eCD+XAnCLaPhHdc/nWwcGB5i+V8WPFDVl9Zje8OX+7E3cn4+CLb6WTgCyGyrlxRSzYPqUOFYlYMXBXMwoM3yfA9O7e20OBzOLsaTix458MM9Q35vMbnzHhvBuGx4XTe2pl9EfsyeAYiu0ngC5GP2Vsa87tfLVpVLs5Pf17li00XScnoDp4GY6H8+5qtmqEHMnRoY6fG+Lfxp5RlKYb/NZwpp6bIEo8WSeALkc+ZGOozu3tVhjQsw9qTEfT97RQxiRkIXT09zdJOIVcI6ANPwjI0vqOlIytarqBHhR6svLySPoF9uBsnz2hqgwS+EAWAnp7CZ80rMLWTByduPabjvKPcjs7Au2tNrKD7WlDT4PeekByfofGN9I0YV3Mc0xtMJzQmlM5bO3Pgdsb+WhBZJ4EvRAHS2bskK/rVJOrZc9r/eoTg8CfvfnChMprtmg8vw+bB73wT90XNnJvh39ofBwsHhu4byvSg6aSkyxJPbpHAF6KAqV2mEBsH+2BhYkD3RceZu/8GTxPe8elY1ybQZAJc3gyHpmdq/JJWJVnZaiVdy3dl2aVl9Avsx/34TOwiEhmmk0/aKorSBmjj6uo6ICQkRNvTESJfio5PZqT/WfZfi8LYQI/2ng709nHGrYTVmw9UVdjwkaZ+fo91UK55pucQeCuQCccmYKhnyI91f6S+Y/1M9yU03vSkrU4G/t+ktIIQOe/KvVhWHAtj05k7JKWkU8PZDl8fJ5pXKoah/msWAZITYGlzzQ3cj/aCfblMjx8eG87oA6O5Gn2VfpX7MazqMAz0MlESQgAS+EKIdxCTkIJ/0G1WHA/jdnQixaxM6FmzFN1qlMLe0vi/Bzy9DQvf09TRH7AXTDL/oOTztOdMPjmZgOsBVCtSjSn1p1DUvGjmT6YAk8AXQryztHSV/dcesuxoGIdCHmGkr8f7HsXp7eOMZ0mbfzcOOwIr2kKZxppdPHr6WRp7e+h2vj32LSb6JkysN5E6DnWy1F9BJIEvhMiUm1FxrDwWzvrgSOKep1LF0ZrePs6871EcY4P/hfupxbB9FNQbBY2/zvKYt2JuMerAKEKehDDAfQCDPQfLEk8GSOALIbLkWVIKG0/fYfmxMEKj4ilkbkT3GqXoWasUxa1MYOtwOL0cOv0GlTtkebzE1EQmnZzExpCNeBf1ZnL9yRQxK5INZ5L/SeALIbKFqqocvvGI5UfD2Hv1IXqKQotKxehTswTeB3xR7l+A/rugmHu2jLf15la+P/49pgamTKw3EZ8SPtnSb34mgS+EyHa3oxNYeTycdaduE5OYgk/RNBY/H4WJsTF6fgfAvFC2jBP6NJRRB0Zx8+lN/Dz8GFRlEPpZvFeQn0ngCyFyTGJyGlvO3mHZ0TCMHpwlwPg77ll6oO+7iZL22VPiPCElgZ9O/MSWm1uoUawGk+tPprBp4WzpO7/Jc4EvD14JkfeoqsrJW9Fc3bmQ3g8m8Vtqc46UHUMfH2fquBZCUZQsj7H5xmZ+PP4j5obmTK4/mZrFa2bDzPOXPBf4f5MrfCHyprg/xmBxegETlMEsS6xLGXtzevs406GaY+bes/uCG09uMOrAKG7F3GJQlUH4efjJEs8LJPCFELkrLRVWd0QNP8r+2suYcdWac5ExWBgb0MnLEd/aTpS2t8h09wkpCfxw/Ae2hm6lVvFaTKw3UZZ4/kcCXwiR+xKiNU/ipj4Hv/2ceWrCimPhbDt/l5Q0lfrl7Old24mG5Yugp5fx5R5VVdl0YxM/nfgJSyNLptSfQvVi1bP/PPIYCXwhhHY8uASLm2rej9tnOxgYE/XsOWtPRrD6RDgPYp9Tys4M39pOdPYqibWZYYaHuP7kOqP2jyLiWQRf1fqKzuU658CJ5B0S+EII7bm8Bfx9oeqH0HYO/O/mbUpaOjsv3Wf50TBOhT3B1FCf9lUd6O3jRIVib6nY+ZL4lHjGHBzDwciDfFrtU/q798+JM8kTJPCFENq170c4OAVaToWafv/5+NLdGFYcDWfz2Ts8T02nVmk7etd2pqlbUQxeV7HzJSnpKXx56Et2hO2gf+X+DK82PFt2BuU1EvhCCO1KT4d1PeH6TvDdAi71XtnsSXyypmLnsXDuPE2khLUJPWs50a16SQpZvKJi50vS0tP48cSPBFwPoEu5LnxZ60v0lIL1nicJfCGE9iXFwuImEB8FfvvB1um1TdPSVfZeecCKY+EcvvEIIwM92niUoI+PM+6Ob36YS1VVZp6eyZKLS2jp0pIf6/6IoV7G7w3kVRL4Qgjd8OgGLGoENqWg/04wMn/rITcePmP50XA2nI4kITmNqqVs6OPjTMvKxTEyeP3V+5ILS5hxegb1HeszvcF0TAxMsvNMdJYEvhBCd4TsgTWdwa295qXo77jOHpuUwobgSFYcC+fWo3gKWxjTo2YpetYsRVGrV4e5/zV/fjj+A15FvZjdaDYWRpnf+59X5LnAl9IKQuRzh2fAnm+g8TdQb2SGDk1PVzkYEsWKY+H8de0h+opCS/fi9PFxwsvJ7j/t/wz9ky8Pf0k5u3LMbzIfWxPbbDoJ3ZTnAv9vcoUvRD71rxeh+0O5ZpnqJuxRPCuPh+MfdJtnSal8UNWB79pVwtLk32v2ByMPMnL/SBwsHFjQdAHFzItlx1noJAl8IYTu+edF6OGad+IWLpvprhKSU1lwIJTZ+0JwsDVlRlfP/1ztB90PYui+oVgbWbOo2SJKWZXK6hnopDcFfsHarySE0B1GZtBtNegbwNrukBST6a7MjAwY0bQcAQNrA9B5/jF+2X2d1LT0f9p4F/NmSfMlJKQm4LvDl2vR17J8CnmNBL4QQntsSkGXFfDkFmz00+zXzwIvJzv+/KQe7T0dmLk3hC4LjhHxOOGfzysVqsTyFsvR19On786+nH14NosnkLdI4AshtMu5LrSYBNcD4a8fs9ydpYkhP3f1ZGY3T0IextFq1iE2no7k7+Xr0jalWdFyBbbGtvjt9uPo3aNZHjOvkMAXQmhf9Y+gmi8cmgaXNmdLl+08HdgxvB5uxa0Y6X+O4b+fJSYxBQAHCweWt1xOScuSDN07lL3he7NlTF0ngS+E0D5FgVbTwLEGbB4E9y9mS7eOtmas9avF6Gbl2H7hHq1mHuLkrWgACpsWZmnzpbgVcmPkgZFsvrE5W8bUZRL4QgjdYGAMXVeCiTX83h3iH2dLt/p6CkMblWX9wNoY6Ct0W3iM6buukZKWjrWxNQubLqRmsZqMPzKeVZdXZcuYukoCXwihOyyLQdfV8OwBrO+jeXNWNqlaypbtn9SjYzVHZu+7Qaf5xwh7FI+ZoRlzGs+hSakmTD41mXln56HL29WzQgJfCKFbHL2gzUy4dRB2fZWtXVsYGzC1cxV+7VGNW1FxvD/rEAFBtzHUM2Rqg6m0K9OOuefmMuXUFNLVrO0Y0kVZe5uwEELkBM/ucP88HJ8LxT3As0e2dv++R3GqlrJhpP9ZPlt/nv3XovjpA3e+q/MdlkaWrLqyimfJz5jgMwEDvfwTkzp5Ji/U0tH2VIQQ2tL0e80rErd+CoXLa678s1EJG1NWf1SLBQdv8vOu65yOeMLPXTwZU30MVsZWzD07l7iUOKbUn4KRvlG2jq0tUlpBCKG7/n4Relqypoa+Zc7UwDkf+ZThv58l7HE8AxuUYUSTcgSErGXSyUnUKl6LmQ1nYmZoliNjZzcprSCEyJvM7KDbGk3ZhXW9IPV5jgzj4WjD9k/q0tW7JPP236TT/KPULtyOH+r8wMn7JxmwewAxzzNf+kFXSOALIXRbscrwwXyIPAnbR2kqbeYAMyMDJnX0YP6H1YiITuD9WYdJjK7K9AbTufL4Cn139uVR4qMcGTu3SOALIXSfWzuo/xmcWQmnFufoUC0qFydweH2qOdkwduMF1h+0ZXKdmUQ+i8R3hy934u7k6Pg5SQJfCJE3vPcFlGsJgWMh7HCODlXM2oSV/WryRasK7L36gK9+f84wt2k8ff4U3x2+hD4NzdHxc4oEvhAib9DTgw4Lwa40+PvC04gcHk7Br34ZNg2ug4WxAV/7P6Oe+TekpafTO7A3lx5dytHxc4IEvhAi7zCxgm5rNU/g/t5T8xKVHFbZwZptw+rRo0Yp1h1NxShqKEZ6pvTf1Z9T90/l+PjZSQJfCJG3FHaFTkvg/gX4Y2iO3cR9kamRPj9+4M4iX28ePbHizuV+mCh2DNoziIORB3N8/OwigS+EyHvKNoUm32jeiXtkRq4N29StKIHD61G9pAvhF/tgkFacT/Z9wp+hf+baHLJCAl8IkTfV+RQqd4Q930LI7lwbtoiVCcv71uCrltWJvtGP9CRnxh4ai/81/1ybQ2ZJ4Ash8iZFgbZzNPv01/eHRzdybWg9PYX+dV3YPKgxReKGkPKsPN8f/575Zxfl2hwyQwJfCJF3GZlpnsTVN9DU0M/hnTsvcythxbZhjejk+CUpMZ78em4W4w9O1tnyyhL4Qoi87e8XoT+NgFnVNE/jxt7NteFNDPX5vr0nc5pOQS/Oh823VtFryxhSs7GWf3aRwBdC5H3OdWHYaajWC4KXw0xP2PG55kUquaRJxeLs7jWTYmorzsUE0milH/di43Jt/Hehk4GvKEobRVEWxsTk/WJFQohcYu0ArX+BYcHg0QVOLoKZVTQvUYnPnRo4RaxM2NV7EvUK9eaJcooWa/qy83LuLjO9iZRHFkLkT9GhcGAqnP8dDEyh5sfgM0xTgTMXzDm1kgWXppKa6ESHEl8zvlU1TAz1c3zcN5VHlsAXQuRvj0Jg/yTNnn0jC6g1CGoPAVObHB/6jxt/8tWRL0hNLEaJpGHM6VaPCsWscnRMqYcvhCi4CpfVPJk7+Bi4NoKDU2CGBxyYAkmxOTp0W9dWzGk8CxOzKB5a/Ezb+dtZevgW6enaudCWK3whRMFy/wL8NRGubQdTW/D5BGr4gbFFjg0Z/CCYIXuGkJJiQvTNvtRzqci0zh4UsTTJ9rHkCl8IIf5WzB26r9G8MtGxBuz9VnNz9+jsHCvG5lXUi6UtlmJukk7R8os5GXmBFjMOsedy7u0iAgl8IURBVaIq9PSH/nuguIdmN88sTzg+H1KSsn04t0JuLGu5DAtjY2zKLMHG9g4frQjiq80XSExOy/bxXkUCXwhRsJWsDr02Qd8dULgcBH4Os6pq3qyVmpytQ5W2Ls2KlisobGbHM5tfaV0zllXHI2gz5zCX7ub8NnQJfCGEAHDygT7boPdWzdO720fB7GqaB7nSUrJtmBIWJVjWYhlOVk4ciZvCyPZJxCam8MGvR1l8KDRHb+hK4AshxItc6kO/QPhwI1gUha2fwBxvOLtG8+KVbFDYtDBLWyylcqHKLLn+HZ+0f0qD8vb8sP0KvktP8iA2+5eUQAJfCCH+S1HAtTF8tAd6+IOJNWweBHNrwvkASM/6mruVkRULmi6gVvFaTA76lnpel/npA3eCwqNpMeMg92ISs+FE/k0CXwghXkdRoFxz8DsAXVeDgQls/Ajm+cClTZCenqXuzQzNmN1oNk2dmjI1aCpPjLeydWhdfGs7U8wq+7dsSuALIcTbKApUbA0fH4LOyzS/C+gDC+rBlW1Zes2ikb4RU+pP4QPXD5h/bj7rw+YwvIkriqJky9RfJIEvhBDvSk8PKn0Ag45Ch8WQmgTresLCBnB9Z6aD30DPgG99vqWXWy/WXF3D+CPjSU3P/vLKEvhCCJFRevrg0RkGn4D28yApBtZ0gcVN4MbeTAW/oih85v0ZQzyHcPL+SZ4kPcn2aUtpBSGEyKq0FM0unoNTIeY2lKwFDb+A0g0y1V1scixWRpkrsialFYQQIifpG4JXb81LWN7/WfP2rRVtYVlrCD+W4e4yG/ZvI4EvhBDZxcAIqveHT85Ayynw6Dr81gJWtIfbp7Q9Owl8IYTIdoYmmheufHIWmv2gqdC5pAms7gx3TmttWhL4QgiRU4zMNG/ZGn4OGn8DkadgUUNY20PzP4FcppOBL++0FULkK8YWUG8kDD8PDb+C8MMwvy74+8LDK7k2DdmlI4QQuS3xKRyfC8fmQnIcVO4A743TvJ0ri2SXjhBC6BJTG822zU/PQ90RcC0Qfq0BmwbC45s5NqwEvhBCaIuZHTT5RrPGX2uwpj7PnOqwZSgkRGf7cBL4QgihbRb20PxHTfDX8IOww2Bolu3DSOALIYSusCwGLSfB0FOarZ3ZTAJfCCF0jb5hjnQrgS+EEAWEBL4QQhQQEvhCCFFASOALIUQBIYEvhBAFhAS+EEIUEBL4QghRQOh08TRFUaKA8Nd8bA28qZxmYeBRtk9KO952rnllzOzoMzN9ZOSYd237Lu3e1Ea+n7o5blb7zOzx2fkddVJV1f6Vn6iqmif/AQvf8nmQtueYW+eaV8bMjj4z00dGjnnXtu/S7k1t5Pupm+Nmtc/MHp8T39FX/cvLSzpbtT2BXKSNc82JMbOjz8z0kZFj3rXtu7QrKN9RbZ2nLn5HM3t8TnxH/0Onl3SyQlGUIPU1NaGF0Db5fgptyMtX+G+zUNsTEOIN5Pspcl2+vcIXQgjxb/n5Cl8IIcQLJPCFEKKAkMAXQogCokAEvqIo5oqiLFcUZZGiKD21PR8hXqYoSmlFUZYoirJe23MR+VeeDXxFUZYqivJQUZSLL/2+haIo1xRFuaEoytj//boDsF5V1QFA21yfrCiQMvIdVVU1VFXV/tqZqSgo8mzgA8uAFi/+QlEUfeBXoCXgBnRXFMUNcARu/69ZWi7OURRsy3j376gQOS7PBr6qqgeB6Jd+XQO48b+rpWTgd6AdEIkm9CEPn7PIWzL4HRUix+W38HPg/6/kQRP0DsBGoKOiKPMoOI+7C930yu+ooiiFFEWZD1RVFGWcdqYm8jsDbU8gmymv+J2qqmo80De3JyPEK7zuO/oYGJjbkxEFS367wo8ESr7wsyNwV0tzEeJV5DsqtCa/Bf4poKyiKC6KohgB3YA/tDwnIV4k31GhNXk28BVFWQscA8orihKpKEp/VVVTgaHATuAK4K+q6iVtzlMUXPIdFbpGiqcJIUQBkWev8IUQQmSMBL4QQhQQEvhCCFFASOALIUQBIYEvhBAFhAS+EEIUEBL4QghRQEjgCyFEASGBL4QQBcT/AdDna8HF1aPOAAAAAElFTkSuQmCC\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.loglog(data1[:,0], data1[:,1],label=\"Euler\")\n", + "plt.loglog(data2[:,0], data2[:,1],label=\"Heun\")\n", + "#plt.loglog(data3[:,0], data3[:,1])\n", + "#plt.loglog(data4[:,0], data4[:,1])\n", + "plt.loglog(data5[:,0], data5[:,1],label=\"MBD\")\n", + "plt.legend()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 40, + "outputs": [], + "source": [], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/Auswertung/requirements.txt b/Auswertung/requirements.txt index 9ee03f4..972860b 100644 --- a/Auswertung/requirements.txt +++ b/Auswertung/requirements.txt @@ -1,4 +1,5 @@ numpy~=1.19.5 seaborn~=0.11.2 matplotlib~=3.3.4 -pandas~=1.1.5 \ No newline at end of file +pandas~=1.1.5 +jupyter~=6.4.6 \ No newline at end of file diff --git a/C++/src/CMakeLists.txt b/C++/src/CMakeLists.txt index 9275e10..37250a3 100644 --- a/C++/src/CMakeLists.txt +++ b/C++/src/CMakeLists.txt @@ -3,11 +3,11 @@ find_package(spdlog) find_package(Eigen3) +add_library(logging logging.cpp) -# Generic test that uses conan libs add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp Simulation.h Integratoren2d.cpp - Calculation.cpp Compute.cpp) + Compute.cpp konvergenz_runner.cpp konvergenz_runner.h) target_include_directories(integratoren PUBLIC .) target_link_libraries( integratoren @@ -15,7 +15,17 @@ target_link_libraries( project_warnings fmt::fmt spdlog::spdlog - Eigen3::Eigen3) + Eigen3::Eigen3 + logging + ) + +target_include_directories(logging PUBLIC .) +target_link_libraries( + logging + PRIVATE project_options + project_warnings + fmt::fmt + spdlog::spdlog) add_executable(main main.cpp) find_package(OpenMP) @@ -28,4 +38,5 @@ target_link_libraries( project_warnings integratoren Eigen3::Eigen3 + fmt::fmt ) \ No newline at end of file diff --git a/C++/src/Calculation.cpp b/C++/src/Calculation.cpp deleted file mode 100644 index 4f47d41..0000000 --- a/C++/src/Calculation.cpp +++ /dev/null @@ -1,94 +0,0 @@ -// -// Created by jholder on 22.10.21. -// - -#include "Calculation.h" - -#include -#include - -void Calculation::run(size_t steps) { - std::ofstream file ("dump.xyz"); - for (size_t step = 0; step < steps; ++step) { - m_integrator(rod, sim); - for (auto &comp : computes) { - comp.eval(rod); - } - file< &Calculation::getComputes() const { - return computes; -} - -const Simulation &Calculation::getSim() const { return sim; } - -Calculation::Calculation( - std::function t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, double length) - : rod(length), - start_rod(rod), - sim(deltaT, seed), - m_integrator(std::move(t_integrator)) { - for (const auto &pair : t_computes) { - computes.emplace_back( - Compute(rod, pair.first, Compute::Force::zero_F, pair.second, sim)); - } -} - -Calculation::Calculation( - inte_force_type t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, force_type force, Compute::Force type_force, - torque_type torque, double length) - : rod(length), start_rod(rod), sim(deltaT, seed) { - m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { - t_integrator(t_rod, t_sim, force, torque); - }; - for (const auto &pair : t_computes) { - computes.emplace_back( - Compute(rod, pair.first, type_force, pair.second, sim)); - } -} - - -Calculation::Calculation( - std::function t_integrator, - const std::vector> &t_computes, - double deltaT, size_t seed, double length) - : rod(length), - start_rod(rod), - sim(deltaT, seed), - m_integrator(std::move(t_integrator)) { - for (const auto &pair : t_computes) { - computes.emplace_back( - Compute(rod, pair.first, Compute::Force::zero_F, pair.second, sim)); - } -} - -Calculation::Calculation( - Calculation::inte_force_type t_integrator, - const std::vector>& t_computes, double deltaT, - size_t seed, Calculation::force_type force, Compute::Force type_force, - Calculation::torque_type torque, double length) - : rod(length), start_rod(rod), sim(deltaT, seed) { - m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { - t_integrator(t_rod, t_sim, force, torque); - }; - for (const auto &pair : t_computes) { - computes.emplace_back( - Compute(rod, pair.first, type_force, pair.second, sim)); - } -} - -void Calculation::reset() { - rod = start_rod; - for (auto &c : computes) { - c.reset(start_rod); - } -} diff --git a/C++/src/Calculation.h b/C++/src/Calculation.h deleted file mode 100644 index a745b6d..0000000 --- a/C++/src/Calculation.h +++ /dev/null @@ -1,57 +0,0 @@ -// -// Created by jholder on 22.10.21. -// - -#pragma once - -#include -#include "Compute.h" -#include "Rod2d.hpp" -#include "Simulation.h" - -class Calculation { - private: - Rod2d rod; - const Rod2d start_rod; - Simulation sim; - std::function m_integrator; - std::vector computes; - - public: - void reset(); - - [[nodiscard]] const Simulation &getSim() const; - - using force_type = - std::function; - using torque_type = std::function; - using inte_force_type = - std::function; - - Calculation( - inte_force_type t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, force_type force, Compute::Force force_type, - torque_type torque, double length = 1.0); - - Calculation( - std::function t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, double length = 1.0); - - Calculation(std::function t_integrator, - const std::vector> &t_computes, - double deltaT, size_t seed, double length = 1.0); - - Calculation(inte_force_type t_integrator, - const std::vector>& t_computes, - double deltaT, size_t seed, force_type force, - Compute::Force force_type, torque_type torque, - double length = 1.0); - - [[nodiscard]] const std::vector &getComputes() const; - - void run(size_t steps); - - [[nodiscard]] const Rod2d &getRod() const; -}; diff --git a/C++/src/Compute.cpp b/C++/src/Compute.cpp index ca0159a..636b56c 100644 --- a/C++/src/Compute.cpp +++ b/C++/src/Compute.cpp @@ -2,26 +2,34 @@ // Created by jholder on 22.10.21. // +//#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE + #include "Compute.h" #include #include +#include #include "Rod2d.hpp" #include "Simulation.h" + void Compute::evalMSD(const Rod2d &rod2D) { + SPDLOG_TRACE("Evaluating MSD"); const auto &new_pos = rod2D.getPos(); auto old_pos = start_rod.getPos(); auto diff = new_pos - old_pos; auto msd = diff.dot(diff); + SPDLOG_TRACE("Feeding {}", msd); agg.feed(msd); } void Compute::evalX(const Rod2d &rod2D) { + SPDLOG_TRACE("Evaluating X"); const auto &new_pos = rod2D.getPos(); auto msd = (new_pos[0]); agg.feed(msd); + SPDLOG_TRACE("Feeding {}", msd); } void Compute::evalX_squared(const Rod2d &rod2D) { @@ -56,7 +64,8 @@ void Compute::eval_empYY(const Rod2d &rod2D) { void Compute::eval(const Rod2d &rod2D) { time_step++; - if (time_step % every == 0) { + SPDLOG_TRACE("resetting: {} time_step {}", resetting, time_step); + if (time_step % every == 0){ if (resetting or time_step == every) { switch (type) { case Type::msd: @@ -78,10 +87,7 @@ void Compute::eval(const Rod2d &rod2D) { evalX_squared(rod2D); break; } - if (resetting) { - start_rod = rod2D; - } - return; + start_rod = rod2D; } } } @@ -94,7 +100,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, time(sim.getMDeltaT() * static_cast(every)) { switch (type) { case Type::msd: { - resetting = false; + resetting = true; type_str = "msd"; switch (force) { case Force::zero_F: @@ -102,22 +108,24 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, break; case Force::const_F: target = 0.0; - break; + break; //TODO missing analytic case Force::harmonic_F: - target = 2 * (1 - 2 * std::exp(- time) + std::exp(-2 * time)) + 2 * (1 - std::exp(-2 * time)); + target = + 2 * (1 - 2 * std::exp(-time) + std::exp(-2 * time)) + + 2 * (1 - std::exp(-2 * time)); break; } break; } case Type::oaf: { - resetting = false; + resetting = true; type_str = "oaf"; target = std::exp(-rod.getDRot() * time); break; } case Type::empxx: { - resetting = false; + resetting = true; type_str = "empxx"; const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; @@ -130,12 +138,12 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, break; case Force::const_F: case Force::harmonic_F: - target = 0.0; + target = 0.0; //TODO missing analytic break; } } break; case Type::empyy: { - resetting = false; + resetting = true; type_str = "empyy"; const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; @@ -148,7 +156,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, break; case Force::const_F: case Force::harmonic_F: - target = 0.0; + target = 0.0; //TODO missing analytic break; } } break; @@ -160,7 +168,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, target = rod.getPos().sum(); break; case Force::const_F: - target = 0.0; + target = 0.0; //TODO missing analytic break; case Force::harmonic_F: target = rod.getPos()[0] * exp(-1 * time); @@ -174,10 +182,8 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, switch (force) { case Force::zero_F: - target = 0.0; - break; case Force::const_F: - target = 0.0; + target = 0.0; //TODO missing analytic break; case Force::harmonic_F: target = pow(rod.getPos()[0], 2) * exp(-2 * time) + @@ -185,10 +191,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, break; } } break; - - } - resetting = false; } const LiveAgg &Compute::getAgg() const { return agg; } @@ -207,6 +210,8 @@ std::ostream &operator<<(std::ostream &os, const Compute &compute) { } void Compute::reset(const Rod2d &rod) { + SPDLOG_TRACE("Resettign Compute for new run"); + time_step = 0; start_rod = rod; } diff --git a/C++/src/Integratoren2d.h b/C++/src/Integratoren2d.h index 4142807..71e96d7 100644 --- a/C++/src/Integratoren2d.h +++ b/C++/src/Integratoren2d.h @@ -9,36 +9,49 @@ class Integratoren2d { using Vector = Eigen::Vector2d; -public: - static void Set1_Euler_f(Rod2d &rod2D, Simulation &sim, const std::function &force, - const std::function &torque); + public: + using integrator_t = std::function + &, + const std::function &)>; + using force_t = const std::function; + static void Set1_Euler_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); static void Set1_Euler(Rod2d &rod2D, Simulation &sim); - static void Set2_Heun_f(Rod2d &rod2D, Simulation &sim, const std::function &force, - const std::function &torque); + static void Set2_Heun_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); static void Set2_Heun(Rod2d &rod2D, Simulation &sim); - static void - Set3_Exact_f(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque); + static void Set3_Exact_f( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque); static void Set3_Exact(Rod2d &rod2D, Simulation &sim); - static void - Set4_BDAS_f(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque); + static void Set4_BDAS_f( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque); static void Set4_BDAS(Rod2d &rod2D, Simulation &sim); - static void - Set5_MBD_f(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque); + static void Set5_MBD_f( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque); static void Set5_MBD(Rod2d &rod2D, Simulation &sim); }; diff --git a/C++/src/LiveAgg.cpp b/C++/src/LiveAgg.cpp index 8b4fc34..a3b6da4 100644 --- a/C++/src/LiveAgg.cpp +++ b/C++/src/LiveAgg.cpp @@ -12,6 +12,9 @@ double LiveAgg::getSEM() const noexcept { double LiveAgg::getMean() const noexcept { return mean; } void LiveAgg::feed(double value) noexcept { + if (value != value) { + return; + } num_of_data += 1; auto delta = value - mean; mean += delta / num_of_data; @@ -21,6 +24,7 @@ void LiveAgg::feed(double value) noexcept { int LiveAgg::getNumPoints() const noexcept { return num_of_data; } std::ostream &operator<<(std::ostream &os, const LiveAgg &agg) { - os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean << " S: " << agg.S; + os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean + << " S: " << agg.S; return os; } diff --git a/C++/src/force_lambdas.h b/C++/src/force_lambdas.h index d805705..a801b74 100644 --- a/C++/src/force_lambdas.h +++ b/C++/src/force_lambdas.h @@ -1,7 +1,7 @@ // // Created by jholder on 27.10.21. // - +#include "Eigen/Dense" #pragma once [[maybe_unused]] const static auto harmonic_Force = [](const Eigen::Vector2d &pos, diff --git a/C++/src/konvergenz_runner.cpp b/C++/src/konvergenz_runner.cpp new file mode 100644 index 0000000..5f68be1 --- /dev/null +++ b/C++/src/konvergenz_runner.cpp @@ -0,0 +1,51 @@ +// +// Created by jholder on 16.12.21. +// +#include "konvergenz_runner.h" +#include +#include +#include "force_lambdas.h" + +void Konvergenz::konvergenz(Compute::Force force, + const Integratoren2d::integrator_t& integrator, + size_t every, std::ofstream& file) { + Integratoren2d::force_t force_lambda = harmonic_Force; + switch (force) { + case Compute::Force::zero_F: + break; + case Compute::Force::const_F: + break; + case Compute::Force::harmonic_F: + break; + } + double delta_t = 1.0 / static_cast(every); + Rod2d rod(1.0); + Simulation sim(delta_t, 1234); + Compute com(rod, Compute::Type::oaf, force, every, sim); + + while (true) { + for (int i = 0; i < 100000; ++i) { + integrator(rod, sim, force_lambda, zero_Torque); + com.eval(rod); + } + std::cout << every << " " << com.getDifference() << " " + << com.getAgg().getSEM() << " \r"; + if (com.getDifference() >= com.getAgg().getSEM() * 10 and + com.getAgg().getSEM() < 0.01) { + std::cout << "\n"; + file << every << " " << com.getDifference() << " " + << com.getAgg().getSEM() << '\n'; + file.flush(); + break; + } + } +} +void Konvergenz::runner(Compute::Force force, + const Integratoren2d::integrator_t& integrator, + std::string filename) { + std::ofstream file; + file.open(filename); + for (size_t i = 1; i < 100; i *= 2) { + konvergenz(force, integrator, i, file); + } +} diff --git a/C++/src/konvergenz_runner.h b/C++/src/konvergenz_runner.h new file mode 100644 index 0000000..287721f --- /dev/null +++ b/C++/src/konvergenz_runner.h @@ -0,0 +1,18 @@ +// +// Created by jholder on 16.12.21. +// + +#ifndef MYPROJECT_KONVERGENZ_RUNNER_H +#define MYPROJECT_KONVERGENZ_RUNNER_H +#include "Compute.h" +#include "Integratoren2d.h" +class Konvergenz { + public: + static void runner(Compute::Force force, + const Integratoren2d::integrator_t& integrator, + std::string filename); + static void konvergenz(Compute::Force force, + const Integratoren2d::integrator_t& integrator, + size_t every, std::ofstream& file); +}; +#endif // MYPROJECT_KONVERGENZ_RUNNER_H diff --git a/C++/src/legacy/Integratoren2d_force.cpp b/C++/src/legacy/Integratoren2d_force.cpp deleted file mode 100644 index 142af06..0000000 --- a/C++/src/legacy/Integratoren2d_force.cpp +++ /dev/null @@ -1,184 +0,0 @@ -// -// Created by jholder on 21.10.21. -// - -#include "Integratoren2d_force.h" -#include "vec_trafos.h" - -using Vector = Eigen::Vector2d; -using Matrix = Eigen::Matrix2d; - -const Vector Integratoren2d_force::unitVec = {1, 0}; - -void Integratoren2d_force::Set1_Euler( - Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - auto e = rod2D.getE(); - // translation - Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise - Vector trans_part = rod2D.getDiff_Sqrt() * rand; - Vector trans_lab = rotation_Matrix(e) * trans_part; - - // Force - Vector F_lab = force(rod2D.getPos(), e); - Vector F_part = rotation_Matrix(e).inverse() * F_lab; - Vector F_trans = rotation_Matrix(e)*rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT(); - - // rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - double M_rot = torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT(); - Vector new_e = e + (rot + M_rot) * ortho(e); - new_e.normalize(); - // apply - rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); - rod2D.setE(new_e); -} - -Vector Integratoren2d_force::Heun_predictor_pos( - const Rod2d &rod2D, Simulation &sim, - const std::function &force) { - auto standard_dev = sim.getSTD(); - Vector rand_pred = {sim.getNorm(standard_dev), - sim.getNorm(standard_dev)}; // gaussian noise - Vector trans_pred_part = - rod2D.getDiff_Sqrt() * - rand_pred; // translation vector in Particle System - Vector trans_pred_lab = rotation_Matrix(rod2D.getE()) * trans_pred_part; - - Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_pred_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab; - Vector F_pred_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_pred_part; - F_pred_trans *= sim.getMDeltaT() ; - - return F_pred_trans + trans_pred_lab; -} - -Vector Integratoren2d_force::Heun_predictor_rot( - const Rod2d &rod2D, Simulation &sim, - const std::function &torque) { - auto std = sim.getSTD(); - double rot_predict = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - const Vector& e = rod2D.getE(); - double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * - rod2D.getDRot() * sim.getMDeltaT(); - Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e); - - return e_change_predict; -} - -void Integratoren2d_force::Set2_Heun( - Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force); - Vector pos_predictor = rod2D.getPos() + delta_pos_predictor; - - Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque); - Vector e_predictor = rod2D.getE() + delta_e_predictor; - e_predictor.normalize(); - - Rod2d pred_rod = rod2D; - pred_rod.setPos(pos_predictor); - pred_rod.setE(e_predictor); - - Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque); - - Vector e_integrated = - rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); - // apply - rod2D.setPos(pos_predictor); - rod2D.setE(e_integrated.normalized()); -} - -void Integratoren2d_force::Set3_Exact( - Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - // translation - Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise - Vector trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - - Vector F_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; - - Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT(); - Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; - - // rotation - double rot = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - Vector e = rod2D.getE(); - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT(); - auto correction = - -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT(); - Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e; - new_e.normalize(); - // apply - rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); - rod2D.setE(new_e); -} - -void Integratoren2d_force::Set4_BDAS( - Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - Vector e = rod2D.getE(); - // translation - Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise - Vector trans_part = rod2D.getDiff_Sqrt() * rand; - Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; - - // Force - Vector F_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; - Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT(); - - // rotation - double rot = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT(); - Vector new_e = Eigen::Rotation2D(rot + M_rot) * e; - new_e.normalize(); - // apply - rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); - rod2D.setE(new_e); -} - -void Integratoren2d_force::Set5_MBD(Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque) { - - Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force); - [[maybe_unused]] Vector pos_predictor = rod2D.getPos() + delta_pos_predictor; - - Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque); - Vector e_predictor = rod2D.getE() + delta_e_predictor; - e_predictor.normalize(); - - - auto std = sim.getSTD(); - Eigen::Vector2d e = rod2D.getE(); - - - Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predictor)*rod2D.getDiff()); - Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL(); - - auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand; - - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predictor)); - //apply - rod2D.setPos(pos_integrated); - rod2D.setE(e_integrated.normalized()); -} diff --git a/C++/src/legacy/Integratoren2d_force.h b/C++/src/legacy/Integratoren2d_force.h deleted file mode 100644 index 5816954..0000000 --- a/C++/src/legacy/Integratoren2d_force.h +++ /dev/null @@ -1,52 +0,0 @@ -// -// Created by jholder on 21.10.21. -// - -#pragma once - -#include "Rod2d.hpp" -#include "Simulation.h" - -class Integratoren2d_force { - public: - static void Set1_Euler( - Rod2d &rod2D, Simulation &sim, - const std::function - &force, - const std::function &torque); - - static void Set2_Heun( - Rod2d &rod2D, Simulation &sim, - const std::function - &force, - const std::function &torque); - - static void Set3_Exact( - Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque); - - static void Set4_BDAS( - Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque); - - static const Eigen::Vector2d unitVec; - - static void Set5_MBD(Rod2d &rod2D, Simulation &sim, - const std::function& force, - const std::function& torque); - -private: - static Eigen::Vector2d Heun_predictor_pos( - const Rod2d &rod2D, Simulation &sim, - const std::function - &force); - - static Eigen::Vector2d Heun_predictor_rot( - const Rod2d &rod2D, Simulation &sim, - const std::function &torque); - - void Set1_Euler(Rod2d &rod2D, Simulation &sim, const std::function &force, - const std::function &torque); -}; diff --git a/C++/src/legacy/Integratoren2d_forceless.cpp b/C++/src/legacy/Integratoren2d_forceless.cpp deleted file mode 100644 index 836220a..0000000 --- a/C++/src/legacy/Integratoren2d_forceless.cpp +++ /dev/null @@ -1,123 +0,0 @@ -// -// Created by jholder on 21.10.21. -// - -#include "Integratoren2d_forceless.h" -#include "vec_trafos.h" - -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) - Eigen::Vector2d e = rod2D.getE(); - // translation - Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; - Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; - Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; - - // rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - Eigen::Vector2d new_e = e + rot * ortho(e); - - // apply - rod2D.setPos(rod2D.getPos() + trans_lab); - rod2D.setE(new_e.normalized()); -} - -void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { - // Predict with Euler - auto std = sim.getSTD(); - Eigen::Vector2d e = rod2D.getE(); - - Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - Eigen::Vector2d trans_part = - rod2D.getDiff_Sqrt() * rand; - Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; - Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab; - - // rotation - double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - 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); // pos with Euler - rod2D.setE(e_integrated.normalized()); -} - -void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { - auto std = sim.getSTD(); - Eigen::Vector2d e = rod2D.getE(); - // translation - Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; - Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; - Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; - - // rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - auto correction = - -0.5 * pow(rod2D.getDRot(), 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()); -} - -// this is a slow implementation to keep the same data structure for performance -// analysis this must be altered -// Normalisation should not be necessary if a proper angular representation -// is used. But with vector e it is done in case of numerical errors -void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { - auto std = sim.getSTD(); - Eigen::Vector2d e = rod2D.getE(); - // translation - auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; - - auto rot = - Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); - Eigen::Rotation2Dd rotation(rot); - auto e_new = (rotation.toRotationMatrix() * e).normalized(); - - rod2D.setPos(rod2D.getPos() + trans_lab); - rod2D.setE(e_new); -} - -void Integratoren2d_forceless::Set5_MBD(Rod2d & rod2D, - Simulation & sim) { - - auto std = sim.getSTD(); - Eigen::Vector2d e = rod2D.getE(); - - //rotation - double rot_predict = sim.getNorm(std) * - rod2D.getDRot_Sqrt(); - auto delta_e_predict = rot_predict * ortho(e); - auto e_predict = (e + delta_e_predict).normalized(); - - - Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predict)*rod2D.getDiff()); - Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL(); - - auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand; - - 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_integrated); - rod2D.setE(e_integrated.normalized()); - -} diff --git a/C++/src/legacy/Integratoren2d_forceless.h b/C++/src/legacy/Integratoren2d_forceless.h deleted file mode 100644 index efeb594..0000000 --- a/C++/src/legacy/Integratoren2d_forceless.h +++ /dev/null @@ -1,25 +0,0 @@ -// -// Created by jholder on 21.10.21. -// - -#pragma once - -#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 void Set5_MBD(Rod2d &rod2D, Simulation &sim); - - static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); -}; diff --git a/C++/src/logging.cpp b/C++/src/logging.cpp new file mode 100644 index 0000000..2bf80ef --- /dev/null +++ b/C++/src/logging.cpp @@ -0,0 +1,20 @@ +// +// Created by jholder on 16.12.21. +// + +#include "logging.hpp" +#include + +void LOGGER::setLogging(LOG_LEVEL level) { + switch (level) { + case LOG_LEVEL::WARN: + spdlog::set_level(spdlog::level::warn); + break; + case LOG_LEVEL::DEBUG: + spdlog::set_level(spdlog::level::debug); + break; + case LOG_LEVEL::TRACE: + spdlog::set_level(spdlog::level::trace); + break; + } +} \ No newline at end of file diff --git a/C++/src/logging.hpp b/C++/src/logging.hpp new file mode 100644 index 0000000..f1d9b42 --- /dev/null +++ b/C++/src/logging.hpp @@ -0,0 +1,15 @@ +// +// Created by jholder on 16.12.21. +// + +#ifndef MYPROJECT_LOGGING_H +#define MYPROJECT_LOGGING_H + +enum LOG_LEVEL{WARN,DEBUG, TRACE}; +class LOGGER{ + public: + static void setLogging(LOG_LEVEL level); +}; + + +#endif // MYPROJECT_LOGGING_H diff --git a/C++/src/main.cpp b/C++/src/main.cpp index 85f0d51..0e0268e 100644 --- a/C++/src/main.cpp +++ b/C++/src/main.cpp @@ -1,189 +1,29 @@ -// -// Created by jholder on 22.10.21. -// - -#include -#include - -#include -#include -#include "Calculation.h" -#include "Integratoren2d.h" -#include "force_lambdas.h" -#include -constexpr size_t SEED = 1364; -constexpr double stepSize = 0.01; -constexpr size_t n_computes = 10; -constexpr size_t num_integratoren = 5; -constexpr size_t num_forces = 3; -constexpr size_t num_length = 1; -constexpr size_t delta_compute = 10; -constexpr size_t numStep = 50000; -inline const std::string header = {"time,val,target,std,numPoints\n"}; - -void run(size_t integrator_index, size_t force_index, size_t length_index) { - using type = std::pair; - std::vector vec({type(Integratoren2d::Set1_Euler_f, "euler"), - type(Integratoren2d::Set2_Heun_f, "heun"), - type(Integratoren2d::Set3_Exact_f, "exact"), - type(Integratoren2d::Set4_BDAS_f, "bdas"), - type(Integratoren2d::Set5_MBD_f, "mbd")}); - auto [integrator, name] = vec[integrator_index]; - - auto force = std::vector< - std::function>( - {zero_Force, const_Force, harmonic_Force})[force_index]; - auto force_name = std::vector( - {"zero_force_", "const_force_", "harmonic_force_"})[force_index]; - auto force_type = std::vector( - {Compute::Force::zero_F, Compute::Force::const_F, - Compute::Force::harmonic_F})[force_index]; - - namespace fs = std::filesystem; - - std::string folder = - "out/" + force_name + name + "_L" + std::to_string(length_index) + "_"; - fs::create_directories("out"); - // double length = std::vector({1.0, 1.5, 2.0})[length_index]; - double length = 1.0; - { - std::vector> computes; - for (size_t i = 1; i < n_computes; ++i) { - computes.emplace_back(Compute::Type::msd, i * delta_compute); - computes.emplace_back(Compute::Type::oaf, i * delta_compute); - computes.emplace_back(Compute::Type::empxx, i * delta_compute); - computes.emplace_back(Compute::Type::empyy, i * delta_compute); - computes.emplace_back(Compute::Type::x, i * delta_compute); - computes.emplace_back(Compute::Type::x_squared, i * delta_compute); - } - - Calculation euler(integrator, computes, stepSize, SEED, force, - force_type, zero_Torque, length); - - for (size_t i = 0; i < numStep; ++i) { - euler.run(n_computes * delta_compute); - euler.reset(); - } - - std::ofstream msdFile(folder + "msd.dat"); - std::ofstream oafFile(folder + "oaf.dat"); - std::ofstream empxxFile(folder + "empxx.dat"); - std::ofstream empyyFile(folder + "empyy.dat"); - std::ofstream xFile(folder + "x.dat"); - std::ofstream xsqFile(folder + "x_squared.dat"); - - xFile << header; - xsqFile << header; - msdFile << header; - oafFile << header; - empxxFile << header; - empyyFile << header; - - for (const auto &com : euler.getComputes()) { - if (com.getType() == Compute::Type::msd) { - msdFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::oaf) { - oafFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::empxx) { - empxxFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::empyy) { - empyyFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::msd) { - msdFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::x) { - xFile << com.getTime() << ", " << com.getAgg().getMean() << ", " - << com.getTarget() << ", " << com.getAgg().getSEM() - << ", " << com.getAgg().getNumPoints() << std::endl; - } - if (com.getType() == Compute::Type::x_squared) { - xsqFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << ", " - << com.getAgg().getNumPoints() << std::endl; - } - } - } - std::cout << "Finished run " << folder << std::endl; +#include "konvergenz_runner.h" +void msd(){ + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set1_Euler_f, + "1Euler_msd_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set2_Heun_f, + "2Heun_msd_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set3_Exact_f, + "3Exact_msd_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set4_BDAS_f, + "4BDAS_msd_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set5_MBD_f, + "5MBD_msd_harmonic.dat"); } - - -void konv(int delta, std::ofstream &file, - Calculation::inte_force_type integrator, size_t seed) { - std::vector> computes; - double step_Size = 1.0 / delta; - computes.emplace_back(Compute::Type::msd, delta); - Calculation euler(integrator, computes, step_Size, seed, harmonic_Force, - Compute::Force::harmonic_F, zero_Torque, 1.0); - - for (int j = 0; j < 1000000; ++j) { - for (int i = 0; i < 1000; ++i) { - euler.run(static_cast(delta)); - euler.reset(); - } - std::cout << euler.getComputes()[0].getAgg().getSEM() << " " - << euler.getComputes()[0].getDifference() << std::endl; - if (euler.getComputes()[0].getAgg().getSEM() * 10 < - euler.getComputes()[0].getDifference()) { - file << delta << " " << euler.getComputes()[0].getAgg().getMean() - << " " << euler.getComputes()[0].getAgg().getSEM() << " " - << euler.getComputes()[0].getTarget() << " " - << euler.getComputes()[0].getDifference() << std::endl; - return; - } - } - file << "timeout: " << delta << " " - << euler.getComputes()[0].getAgg().getMean() << " " - << euler.getComputes()[0].getAgg().getSEM() << " " - << euler.getComputes()[0].getTarget() << " " - << euler.getComputes()[0].getDifference() << std::endl; +void oaf(){ + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set1_Euler_f, + "1Euler_oaf_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set2_Heun_f, + "2Heun_oaf_harmonic.dat"); + // Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set3_Exact_f, + // "3Exact_oaf_harmonic.dat"); + //Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set4_BDAS_f, + // "4BDAS_oaf_harmonic.dat"); + Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set5_MBD_f, + "5MBD_oaf_harmonic.dat"); } -Calculation::inte_force_type integrator(int index) { - if (index % 5 == 0) return Integratoren2d::Set1_Euler_f; - if (index % 5 == 1) return Integratoren2d::Set2_Heun_f; - if (index % 5 == 2) return Integratoren2d::Set3_Exact_f; - if (index % 5 == 3) return Integratoren2d::Set4_BDAS_f; - if (index % 5 == 4) return Integratoren2d::Set5_MBD_f; - return Integratoren2d::Set1_Euler_f; -} -void konv_runner(int index) { - std::ofstream file; - file.open(fmt::format("out/{}_konv.out", index)); - for (int i = 1; i < 256; i *= 2) { - konv(i, file, integrator(index), static_cast(index / 5)); - } - file.close(); -} - -int main(int argc, char *argv[]) { - int index = 0; - if (argc >1) index = std::stoi(argv[1]); - konv_runner(index); - +int main(){ + oaf(); return EXIT_SUCCESS; -#pragma omp parallel for default(none) - for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) { - run(j % (num_integratoren), - (j % (num_integratoren * num_forces)) / num_integratoren, - j / (num_integratoren * num_forces)); - } } \ No newline at end of file diff --git a/C++/test/CMakeLists.txt b/C++/test/CMakeLists.txt index fd13677..ba9a7a4 100644 --- a/C++/test/CMakeLists.txt +++ b/C++/test/CMakeLists.txt @@ -5,13 +5,13 @@ 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 PUBLIC Catch2::Catch2 logging) target_link_libraries(catch_main PRIVATE project_options) set(TEST_FILES test_LiveAgg.cpp test_Rod2d.cpp - test_Compute.cpp test_Simulation.cpp test_Integratoren.cpp) + test_Compute.cpp) add_executable(tests ${TEST_FILES}) target_link_libraries(tests PRIVATE project_warnings project_options catch_main Eigen3::Eigen3 integratoren) diff --git a/C++/test/catch_main.cpp b/C++/test/catch_main.cpp index db3af94..e089537 100644 --- a/C++/test/catch_main.cpp +++ b/C++/test/catch_main.cpp @@ -1,3 +1,11 @@ -#define CATCH_CONFIG_MAIN // This tells the catch header to generate a main +#define CATCH_CONFIG_RUNNER #include + +#include "logging.hpp" + +int main( int argc, char* argv[] ) { + LOGGER::setLogging(LOG_LEVEL::TRACE); + auto result = Catch::Session().run( argc, argv ); + return result; +} \ No newline at end of file diff --git a/C++/test/test_Compute.cpp b/C++/test/test_Compute.cpp index 5437766..3151667 100644 --- a/C++/test/test_Compute.cpp +++ b/C++/test/test_Compute.cpp @@ -2,19 +2,41 @@ // Created by jholder on 22.10.21. // -#include "Compute.h" - #include +#include "Compute.h" TEST_CASE("Compute") { Rod2d rod(1.0); + Rod2d rod_old(1.0); Simulation sim(0.1, 1); - auto com = Compute(rod, Compute::Type::msd, 10, sim); - SECTION("Mean of same values") { - for (int i = 0; i < 100; ++i) { + SECTION("MSD") { + Compute com(rod, Compute::Type::msd, Compute::Force::zero_F, 2, sim); + for (int i = 0; i < 11; ++i) { + rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1}); com.eval(rod); } - CHECK(com.getAgg().getNumPoints() == 10); - CHECK(com.getAgg().getMean() == 0.0); + CHECK(com.getAgg().getNumPoints() == 5); + CHECK(com.getAgg().getMean() == 8); + } + SECTION("X") { + Compute com(rod, Compute::Type::x, Compute::Force::zero_F, 2, sim); + + for (int i = 0; i < 11; ++i) { + rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1}); + com.eval(rod); + } + CHECK(com.getAgg().getNumPoints() == 1); + CHECK(com.getAgg().getMean() == 3); + + for (int i = 0; i < 11; ++i) { + rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1}); + com.eval(rod); + if ((i + 1) % 2 == 0) { + rod = rod_old; + com.reset(rod_old); + } + } + CHECK(com.getAgg().getNumPoints() == 5); + CHECK(com.getAgg().getMean() == 3); } } diff --git a/C++/test/test_Integratoren.cpp b/C++/test/test_Integratoren.cpp deleted file mode 100644 index ef59e70..0000000 --- a/C++/test/test_Integratoren.cpp +++ /dev/null @@ -1,62 +0,0 @@ - -// -// Created by jholder on 24.10.21. -// -#include "Calculation.h" -#include "Compute.h" -#include "Integratoren2d.h" -#include "force_lambdas.h" - -#include -#include - -TEST_CASE("Integrator Check") { - const size_t SEED = Catch::rngSeed(); - const double length = GENERATE(1.0, 1.4, 2.0); - constexpr int steps = 10000; - constexpr double delta = 0.1; - - using type = std::tuple, Calculation::inte_force_type, std::string>; - auto[integrator, integrator_force, name] = GENERATE( - type(Integratoren2d::Set1_Euler, Integratoren2d::Set1_Euler_f, "Euler"), - type(Integratoren2d::Set2_Heun, Integratoren2d::Set2_Heun_f, "Heun"), - type(Integratoren2d::Set3_Exact, Integratoren2d::Set3_Exact_f, "Exact"), - type(Integratoren2d::Set4_BDAS, Integratoren2d::Set4_BDAS_f, "BDAS"), - type(Integratoren2d::Set5_MBD, Integratoren2d::Set5_MBD_f, "MBD")); - - Calculation calc(integrator, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED, length); - calc.run(steps); - SECTION("ForceLess") { - CAPTURE(name); - CAPTURE(length); - size_t i = 0; - for (const auto &c: calc.getComputes()) { - CAPTURE(c); - CHECK(c.getDifference() <= delta * c.getAgg().getMean()); - ++i; - } - } - - Calculation calc_f(integrator_force, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED, zero_Force, zero_Torque, length); - calc_f.run(steps); - - SECTION("ForceLess") { - CAPTURE(length); - CAPTURE(name); - for (const auto &c: calc_f.getComputes()) { - CAPTURE(c); - CHECK(c.getDifference() <= delta * c.getAgg().getMean()); - } - } - -} diff --git a/C++/test/test_Rod2d.cpp b/C++/test/test_Rod2d.cpp index 1a77b1a..2a84f95 100644 --- a/C++/test/test_Rod2d.cpp +++ b/C++/test/test_Rod2d.cpp @@ -7,10 +7,12 @@ #include "Rod2d.hpp" #include "vec_trafos.h" + TEST_CASE("Rods") { Rod2d sphere(1); Rod2d rod(2); SECTION("Checking translation") { + rod.setPos({0,0}); rod.setPos(rod.getPos() + Eigen::Vector2d(1, -1)); auto newPos = rod.getPos(); REQUIRE(newPos[0] == 1); diff --git a/C++/test/test_Simulation.cpp b/C++/test/test_Simulation.cpp deleted file mode 100644 index 97f19a7..0000000 --- a/C++/test/test_Simulation.cpp +++ /dev/null @@ -1,28 +0,0 @@ -// -// 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