diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d5067b..ad8436e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,9 +1,18 @@ # CMakeList.txt for linear code -cmake_minimum_required (VERSION 2.8.9) - +cmake_minimum_required (VERSION 3.1) project(Timedomain1DEuler) -add_definitions(-DTRACERNAME=timedomaineuler -DARMA_USE_BLAS -DARMA_USE_LAPACK) +FIND_PACKAGE(SWIG REQUIRED) +INCLUDE(${SWIG_USE_FILE}) +FIND_PACKAGE(PythonLibs) + +INCLUDE_DIRECTORIES(${PYTHON_INCLUDE_PATH}) +INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) + +SET(CMAKE_SWIG_FLAGS "-Wall") + + +add_definitions(-DTRACERNAME=timedomaineulertracer -DARMA_USE_BLAS -DARMA_USE_LAPACK) ############################## Not so often changed @@ -15,7 +24,7 @@ add_definitions(-DTRACERNAME=timedomaineuler -DARMA_USE_BLAS -DARMA_USE_LAPACK) set (CMAKE_GCC " -Wno-unused-function -ffunction-sections -fdata-sections -Wno-unused-local-typedefs -Wno-empty-body") -set (CMAKE_GENERAL "${CMAKE_GENERAL} -std=c++11 -pipe -O2 -fPIC -Wall \ +set (CMAKE_GENERAL "${CMAKE_GENERAL} -march=native -std=c++11 -pipe -O2 -fPIC -Wall \ -Wextra -Wno-unused-parameter \ -Wno-unused-variable -Wno-unused-but-set-variable \ -Wno-return-local-addr -Wno-cpp -Wno-address") @@ -25,6 +34,7 @@ set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CMAKE_GENERAL} ${CMAKE_GCC}") # link_directories(${link_directories} /usr/local/lib) set (tdir /home/anne/wip/tmtubes/src) include_directories( + # ${PYTHON_INCLUDE} ${tdir}/nonlinear/sys ${tdir}/nonlinear/seg ${tdir}/nonlinear/geom @@ -36,9 +46,14 @@ include_directories( ) AUX_SOURCE_DIRECTORY(src src) link_directories(/home/anne/bin/lib) -add_executable(timedomaineuler timedomain.cpp ${src}) -target_link_libraries(timedomaineuler nonlin armadillo boost_program_options) +add_library(sources ${src}) +set_source_files_properties(timedomain.i PROPERTIES CPLUSPLUS ON) +set_source_files_properties(timedomain.i PROPERTIES SWIG_FLAGS "-py3") +set_source_files_properties(timedomain.i PROPERTIES SWIG_FLAGS "-Wextra") + +swig_add_module(timedomaineuler python timedomain.i) +swig_link_libraries(timedomaineuler sources nonlin armadillo ${PYTHON_LIBRARIES}) diff --git a/Plot data.ipynb b/Plot data.ipynb new file mode 100644 index 0000000..697d32f --- /dev/null +++ b/Plot data.ipynb @@ -0,0 +1,75 @@ +{ + "metadata": { + "name": "", + "signature": "sha256:cf9769ed030dba90588e4e94007d226dd2897d23c8792fd8974c73bb04381b2c" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from tmtubes import gas\n", + "air=gas('air')\n", + "L=1.\n", + "c0=air.cm(293.15)\n", + "fres=c0/(4*L)\n", + "print(\"fres:\",fres)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "fres: 85.7847876535731\n" + ] + } + ], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "plot(x,u)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 3, + "text": [ + "[]" + ] + }, + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEACAYAAABCl1qQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVNWZx/Hva7MKDEQwbN0KDqiAhEHC4kqrqIgIokZB\nExQ1kCFEHaOCOiYYYzYNMYgSY9ThGQVEQUAF2aTVRyNEYAC1UVEhLBGXSFSEh+2dP85tu23pqqK3\nW8vv8zz3oerWubfeutr11jnnnnPM3REREUnkkLgDEBGR9KdkISIiSSlZiIhIUkoWIiKSlJKFiIgk\npWQhIiJJJU0WZtbfzNaZ2TtmNraCMhOj11ebWfdkx5rZXWZWHJWfZWZNo/3tzGynma2Ktvur40OK\niEjVJEwWZpYHTAL6A52BYWbWqVyZAUAHd+8IjAQmp3DsQqCLu3cD3gZuLnPK9e7ePdpGV/UDiohI\n1SWrWfQifHlvcPc9wHRgcLkyg4ApAO6+DGhmZq0SHevui9x9f3T8MiC/Wj6NiIjUiGTJoi2wqczz\nzdG+VMq0SeFYgCuBeWWet4+aoIrM7OQk8YmISC2ok+T1VOcCscq8uZndCux296nRrq1Agbt/ambH\nA7PNrIu7f16Z84uISPVIliy2AAVlnhcQagiJyuRHZeomOtbMrgAGAGeU7HP33cDu6PFKM3sX6Ais\nLPuGZqYJrUREKsHdK/XjHnevcCMkk3eBdkA94P+ATuXKDADmRY/7AK8mO5bQ6f0G0KLcuVoAedHj\nowjJpdkB4nIJfv7zn8cdQtrQtSila1FK16JU9N2Z8Hu/oi1hzcLd95rZGGABkAc85O7FZjYqev0B\nd59nZgPMbD2wAxiR6Njo1PdGCWSRmQH81cOdT32B281sD7AfGOXu2w8i94mISA1I1gyFu88H5pfb\n90C552NSPTba37GC8jOBmcliEhGR2qUR3BmusLAw7hDShq5FKV2LUroW1cM8Axc/MjPPxLhFROJk\nZpXu4FbNQkREkkraZyGl9u2DDz+ELVtg69bw+MMP4aOP4J//hH/9C7Zvh88+gy+/hB07YOdO2LMn\nbHv3gjtYlNfz8qB+/bA1aACNGkHjxmFr2hQOOyxszZtDq1bQsmX4Nz8fWrQoPY+ISE1TM1Q5O3fC\nO+/AW2/B22/Du+/C+++HbetWaNYsfFm3aQPf/nbYDj88fKk3axa2Jk3CF3+jRtCwIdSrB3XqhK3k\nC949JI/du8O2a1dILl98AZ9/HhLPP/8Zto8/hm3bwvaPf8DmzSEZtW0L7drBv/976XbssdCxY3hP\nEZGyqtIMlbPJYv9+eO89WLkS1qyB118P25Yt0L49HHMMHH10+AJu3z5sBQWhFpAOduwISWPDhpDQ\n3nuvNMlt2ABHHAHHHQff+Q507Qr/8R9w1FGqjYjkMiWLFGzZAq++CsuWwfLlsGpVqAV07w7duoUv\n1OOOgw4dQg0gk+3eDevXh+S3Zg2sXh0+75dfwvHHw3e/C336wAknhKYtEckNShbluIdf2C+8AC+9\nFLYvv4TevUu3Hj1CX0Au2bYNVqyAv/0N/vrXkDi/9S045RTo2xdOPTXUpFT7EMlOShaEL8LnnoPF\ni+H550PtoLAwfBGeckpoUtKX4Nft3w/r1oVk+uKLIbkCnHEG9OsXttat441RRKpPTiaL/fudtWth\n1ix49tnQXt+vH5x5Zviy0y/kg+cemq8WL4YlS0LSPeIIOOccGDAgNFtlehOdSC7LyWRxzDHOrl1w\n0UUwcCCcdBLUrRt3ZNll797QVDV/PsybB3//e0gagwZB//7hFl8RyRw5mSyWLXN69lTtoTZt2gRP\nPw1z5oQ+j9NPD8n6vPPCuBARSW85mSwyMe5s8umnIXHMnAlFRSFxXHppqOU1bBh3dCJyIEoWEqvt\n2+Gpp2DatHCn1eDBcPnl4Q6rQzShjEjaULKQtPHBByFpTJkSah/Dh8OVV4ZBjSISLyULSUurV8Mj\nj8Cjj4bBj1dfDUOGaCoSkbgoWUha27ULZs+GP/8Z3nwz1DRGjgzzWolI7dEU5ZLWGjSAoUPDuI2i\nojCavkeP0LexZEkY3yEi6U01C4nFjh3w2GMwcWJ4fu218IMfhMQiIjVDzVCSsdxDjeMPf4DXXoPR\no8PWokXckYlkHzVDScYyC9OzPPMMLF0aBv4dfTT85CewcWPc0YlICSULSRudOsGDD8Ibb8Chh4bp\n1IcPh+LiuCMTESULSTutW8NvfxsWdTrmmDC475JLYO3auCMTyV1KFpK2mjWDW28NqwB+97thRuEL\nL1TSEImDkoWkvcaN4cYbQ9I48cSQNIYOVfOUSG1SspCMceih8NOfhjU3uncPzVMjRqgjXKQ2KFlI\nxmncGMaODQte5eeHjvDrroOPPoo7MpHspWQhGatpU7jjjjCFyN694W6qX/0qjBAXkeqlZCEZr2VL\nmDQpLMi0cmW4g2rKlLDGuIhUD43glqzzyiuhb2PPnjAy/JRT4o5IJD1oug+Rctxh+nQYNw569oS7\n79YstyKa7kOkHDMYNgzWrYNu3cIst+PHqz9DpLKULCSrNWwIt90Gq1aFcRmdO8OsWZoWXeRgqRlK\ncsrSpTBmDBx5ZJgevUOHuCMSqT1qhhJJ0WmnhVpGYSH06ROapnbtijsqkfSnZCE5p149uOmmkDTW\nrAl9GkuXxh2VSHpLmizMrL+ZrTOzd8xsbAVlJkavrzaz7smONbO7zKw4Kj/LzJqWee3mqPw6Mzur\nqh9QpCIFBaH/4q674Ior4PLL4eOP445KJD0lTBZmlgdMAvoDnYFhZtapXJkBQAd37wiMBCancOxC\noIu7dwPeBm6OjukMXBKV7w/cb2aq/UiNGjQorKHRvDkcdxxMnaoOcJHykn0R9wLWu/sGd98DTAcG\nlyszCJgC4O7LgGZm1irRse6+yN1LxtcuA/Kjx4OBae6+x903AOuj84jUqMaNYcIEmDsXfv1rGDgw\nrNonIkGyZNEWKPsnsznal0qZNikcC3AlMC963CYql+wYkRrRqxesWAEnnBAmKPzzn1XLEAGok+T1\nVP9MKjd83OxWYLe7Tz3YGMaPH//V48LCQgoLCysTgsg31KsH//3fMGRImAJ9xoyw3Gv79nFHJnJw\nioqKKCoqqpZzJUsWW4CCMs8L+Pov/wOVyY/K1E10rJldAQwAzkhyri0HCqxsshCpCV26hHmmJkwI\nNY4774Qf/jCMDhfJBOV/SN9+++2VPleyZqjXgI5m1s7M6hE6n+eWKzMXGA5gZn2A7e6+LdGxZtYf\nuBEY7O67yp1rqJnVM7P2QEdgeaU/nUgV1akTbrN94YVQuzjnHNhc/ueSSA5ImCzcfS8wBlgAvAk8\n7u7FZjbKzEZFZeYB75nZeuABYHSiY6NT3ws0BhaZ2Sozuz865k1gRlR+PjBaQ7UlHXTuHGoZJ50U\n+jIee0x9GZJbNN2HyEFatQouuwy6doXJk+Gww+KOSCQ1mu5DpBZ17x7umGrdGr7zHVi0KO6IRGqe\nahYiVbB4cRj9fcklYUnX+vXjjkikYqpZiMSkXz9YvRrefx969w7rgYtkIyULkSpq3hxmzgxTn/ft\nq4F8kp3UDCVSjYqLYehQOOaYkDSaNYs7IpFSaoYSSROdOsGyZdCqVegIf/XVuCMSqR6qWYjUkDlz\nYORIuPFGuP56OEQ/zSRmValZKFmI1KCNG8OdUi1awJQpoX9DJC5qhhJJU0ceCS+9BMceG0Z+L1sW\nd0QilaNkIVLD6taFu++GiRPhvPPCv6oYS6ZRM5RILXrvPbjoIujQAR56CJo0iTsiySVqhhLJEEcd\nFSYkbNo0THteXJz8GJF0oGQhUssaNAjTnd9wA5x6KjzxRNwRiSSnZiiRGK1YEZqlLroorP1dJ9ly\nZCJVoFtnRTLYJ5+EUd/uMH16uM1WpCaoz0IkgzVvDvPnQ48e0LNnWC9DJN2oZiGSRmbMgB//GO69\nN9Q2RKqTmqFEssjq1XD++SFZ/PKXkJcXd0SSLZQsRLLMxx/DxRdDw4YwdWq41VakqtRnIZJlWrSA\nBQugfXvo0wfeeSfuiCTXKVmIpKm6dWHSJLjuOjj5ZFiyJO6IJJcpWYikuVGj4PHH4bLLYPLkuKOR\nXKU+C5EM8e67MHAgnHkmTJigAXxy8NTBLZIjtm8PHd95eWEAnzq+5WCog1skRzRrBs8+Gzq+Tzop\nLK4kUhuULEQyTN26cN99cPXVcOKJsHx53BFJLlAzlEgGmzsXrroK/vQnuPDCuKORdKc+C5EctnIl\nDB4cbrG9/nqwSn0VSC5QshDJcZs2wYAB0Lcv3HOP7pSSA1OyEBH+9a+wLkaDBuFOqUaN4o5I0o3u\nhhIRmjYNd0q1aAGFhfDhh3FHJNlEyUIki9SrBw8/DOecE+6U0pxSUl3UsimSZczgF7+A/Pywxvfs\n2dC7d9xRSaZTn4VIFnvmGRgxAh55JEwVIrlNfRYickADB4aEcfXVoXlKpLKSJgsz629m68zsHTMb\nW0GZidHrq82se7Jjzex7ZvaGme0zs+PL7G9nZjvNbFW03V/VDyiS63r3hhdegDvuCCvvqVIulZGw\nz8LM8oBJQD9gC/A3M5vr7sVlygwAOrh7RzPrDUwG+iQ5di0wBHjgAG+73t27H2C/iFTSMcfAK6+E\nju9t2+CPf4RD1K4gByHZ/y69CF/eG9x9DzAdGFyuzCBgCoC7LwOamVmrRMe6+zp3f7saP4eIJNG6\ndahhrFkT1sbYvTvuiCSTJEsWbYFNZZ5vjvalUqZNCsceSPuoCarIzE5OobyIpKhp07Bc665doT/j\niy/ijkgyRbJkkWrrZnXNRrMVKIiaoa4HpppZk2o6t4gQRng/8QQccQSccQZ88kncEUkmSDbOYgtQ\nUOZ5AaGGkKhMflSmbgrHfo277wZ2R49Xmtm7QEdgZfmy48eP/+pxYWEhhYWFCT+IiJSqUwcefBDG\njQtjMRYuhLap1PsloxQVFVFUVFQt50o4zsLM6gBvAWcQfvUvB4YdoIN7jLsPMLM+wD3u3ifFY5cC\nN7j7iuh5C+BTd99nZkcBLwLHufv2cnFpnIVINfnd78La3gsXQseOcUcjNakq4ywS1izcfa+ZjQEW\nAHnAQ+5ebGajotcfcPd5ZjbAzNYDO4ARiY6NAh4CTARaAM+a2Sp3PwfoC9xuZnuA/cCo8olCRKrX\nTTdB8+Zhxtr586Fbt7gjknSkEdwiAoR+jDFj4KmnwrxSkn00gltEqux734MpU+D882HRorijkXSj\nZCEiX+nfH2bNCuMwnnoq7mgknWjWWRH5mpNPDn0X554LX34ZEoeIkoWIfEOPHrBkCZx9NuzYASNH\nxh2RxE3JQkQOqEsXKCqCM88MI72vvz7uiCROShYiUqEOHeDFF8NI75074dZb445I4qJkISIJFRSE\nCQj79QsJ4447wmp8klt0N5SIJNW6dWiSeuYZuOEGrYmRi5QsRCQlhx8OS5fCSy/BNdfA/v1xRyS1\nSclCRFL2rW+FAXsrVsB//qcSRi5RshCRg1KyJkZxMVx1FezbF3dEUhuULETkoDVpEgbubdwII0Yo\nYeQCJQsRqZRGjUKH9z/+AcOHw969cUckNUnJQkQq7dBDYe5c+Phj+P73lTCymZKFiFRJw4YwZw58\n9hlceins2RN3RFITlCxEpMoaNAiz1e7YoYSRrZQsRKRalCSML79UwshGShYiUm3q1y9NGMOGKWFk\nEyULEalWJQlj586wFoY6vbODkoWIVLv69WHmTPj8c90llS2ULESkRjRoEJZm/fRT+MEPlDAynZKF\niNSYBg1g9uwwDuPKKzXSO5MpWYhIjSoZh7F5M/zwh5p8MFMpWYhIjTv0UHj6aVi/XrPVZiolCxGp\nFY0awbPPwuuvh/UwtIBSZlGyEJFa06QJzJsHy5drxb1Mo2QhIrWqZD2M55+HW29VwsgUdeIOQERy\nT8mKe4WF4Y6pn/0s7ogkGSULEYlFixawZAn07Rs6wG+4Ie6IJBElCxGJTcuWsHgxnHpqSBijR8cd\nkVREyUJEYpWfX1rDaNgwLNMq6UfJQkRi1759qGGcdlpIGEOHxh2RlKdkISJp4eij4bnn4Mwzw5iM\n886LOyIpS7fOikja6No1jPS+6qpQ05D0oWQhImmlZ0948smweNLLL8cdjZRQshCRtHPqqfDoozBk\nCKxaFXc0AikkCzPrb2brzOwdMxtbQZmJ0eurzax7smPN7Htm9oaZ7TOz48ud6+ao/DozO6sqH05E\nMtfZZ8PkyXDuubBuXdzRSMIObjPLAyYB/YAtwN/MbK67F5cpMwDo4O4dzaw3MBnok+TYtcAQ4IFy\n79cZuAToDLQFFpvZ0e6uOSpFctCFF4bV9s46C156CY48Mu6Icleyu6F6AevdfQOAmU0HBgPFZcoM\nAqYAuPsyM2tmZq2A9hUd6+7ron3l328wMM3d9wAbzGx9FMOrlf2AIpLZrrgCPvsM+vULCaNVq7gj\nyk3JmqHaApvKPN8c7UulTJsUji2vTVTuYI4RkSx3zTVhadazzw7LtErtS5YsUp0P8htVhGqkOSlF\nhNtug9NPD30YO3bEHU3uSdYMtQUoKPO8gK//8j9QmfyoTN0Ujk32fvnRvm8YP378V48LCwspLCxM\ncmoRyWRm8PvfhzEYF1wAc+dC/fpxR5XeioqKKCoqqpZzmSeYTN7M6gBvAWcAW4HlwLADdHCPcfcB\nZtYHuMfd+6R47FLgBndfET3vDEwl9FO0BRYTOs+/FqSZld8lIjli7164+GKoUwemTYO8vLgjyhxm\nhrtXqiUoYTOUu+8FxgALgDeBx9292MxGmdmoqMw84L2oM/oBYHSiY6OAh5jZJqAP8KyZzY+OeROY\nEZWfD4xWVhCRsurUgalT4ZNP4Mc/1uJJtSVhzSJdqWYhIp9/HvowzjoL7rwz7mgyQ1VqFppIUEQy\nUpMmMH8+nHIKNG8O118fd0TZTclCRDJWixawcCGcfHJ4PHx43BFlLyULEcloBQVhavPTTgs1jHPP\njTui7KSJBEUk43XqBHPmhNHemqm2ZihZiEhW6N07zFR7wQXw+utxR5N9lCxEJGucfTb84Q9wzjmw\ncWPc0WQX9VmISFa59FL48EPo3z9MPNiiRdwRZQeNsxCRrDRuHBQVwZIlYU1vqdo4CyULEclK7jBi\nBHz0EcyeDXXrxh1R/Gpsug8RkUxlBg8+GJLGqFGaFqSqlCxEJGvVrQtPPBHujrrttrijyWzq4BaR\nrNaoETz7LJx0ErRuHSYflIOnZCEiWe/ww2HBgjAtSJs2MGRI3BFlHiULEckJ7dvD00+HW2q//e1Q\n05DUqc9CRHLG8ceHUd4XXgjFxcnLSyklCxHJKWedBb/7XRjlvXVr3NFkDjVDiUjOGT4cNm8OM9S+\n+GJYG0MS06A8EclJ7vCjH8GGDfDMM7kxaE8juEVEKmHvXjj//HC31MMPh4F82UwjuEVEKqFOHXj8\n8TBob/z4uKNJb+qzEJGc1qhRaIY64QRo1y7MJyXfpGQhIjmvZUuYNw/69oX8fDjzzLgjSj9qhhIR\nAY49Fp58Ei67DNasiTua9KNkISISOeUUuPdeGDgQtmyJO5r0omYoEZEyLrkE3n8/JAyNwSilW2dF\nRMopWQNj82aYOzfcNZUNdOusiEg1MoP77oP9++EnP9HCSaBkISJyQHXrwowZ8PLLMGFC3NHEL0sq\nVyIi1e/f/i0snHTCCXDUUbm9Dob6LEREklixIqyDMW8e9OwZdzSVpz4LEZEa1KMH/OUvYR6pjRvj\njiYeaoYSEUnB4MHw3nvhltqXXw5NVLlEzVAiIilyh9Gjw7TmTz+debfUqhlKRKQWmMHEieGW2muv\nza1bapUsREQOQskttUVFIXHkigyrRImIxK9p09Jbajt2hAED4o6o5iWtWZhZfzNbZ2bvmNnYCspM\njF5fbWbdkx1rZoeZ2SIze9vMFppZs2h/OzPbaWarou3+6viQIiLVrV07mDkTrrgC1q6NO5qalzBZ\nmFkeMAnoD3QGhplZp3JlBgAd3L0jMBKYnMKx44BF7n40sCR6XmK9u3ePttFV/YAiIjXlxBPhnnvg\nvPNg27a4o6lZyWoWvQhf3hvcfQ8wHRhcrswgYAqAuy8DmplZqyTHfnVM9O/5Vf4kIiIxuPRSuPzy\nMAZj1664o6k5yZJFW2BTmeebo32plGmT4NiW7l6Sh7cBLcuUax81QRWZ2cnJP4KISLzGj4cjj4Sr\nr87eO6SSdXCn+rFTuW/XDnQ+d3czK9m/FShw90/N7Hhgtpl1cffPyx83vszq6oWFhRQWFqYYqohI\n9TKDRx4Jy7L++tdwyy1xRxQUFRVRVFRULedKliy2AAVlnhcQagiJyuRHZeoeYH/J2lPbzKyVu39g\nZq2BDwHcfTewO3q80szeBToCK8sHVjZZiIjErWFDmDMHevcOS7RecEHcEX3zh/Ttt99e6XMla4Z6\nDegY3aVUD7gEmFuuzFxgOICZ9QG2R01MiY6dC1wePb4cmB0d3yLqGMfMjiIkivcq/elERGpR69Yh\nYfzoR7DyGz9xM1vCmoW77zWzMcACIA94yN2LzWxU9PoD7j7PzAaY2XpgBzAi0bHRqX8DzDCzq4AN\nwMXR/lOBX5jZHmA/MMrdt1fj5xURqVHdu8Of/hQ6vJctCwkkG2huKBGRGvDLX4b5o4qKQhNVOqjK\n3FBKFiIiNcA93FZ7yCHw6KOhEzxumkhQRCTNmMHDD8Pbb4c7pDKd5oYSEakhJXdI9eoFXbqENTEy\nlZqhRERq2PLlcO658Pzz0LVrfHGoGUpEJI316gV//GOoWXz0UdzRVI5qFiIiteTmm+GVV2DRIqhX\nr/bfX3dDiYhkgH37wviL/HyYPLn231/NUCIiGSAvDx57DF54IZ5kURWqWYiI1LL16+Gkk8LyrH37\n1t77qmYhIpJBOnQINYyhQ+H99+OOJjVKFiIiMejXD8aNC30YX3wRdzTJqRlKRCQm7nDVVfDZZ/DE\nEzU/JYiaoUREMpBZ6OjesgXuvDPuaBLTdB8iIjGqXx9mzYKePcPo7nSdEkTNUCIiaWD5chg4MExp\n3rlzzbyHmqFERDJcr15w112hZvHpp3FH802qWYiIpJHrroO33oJnngmD+KqTahYiIlni7rth9264\n9da4I/k6JQsRkTRSpw48/ngY3T19etzRlFIzlIhIGlq9OgzcW7wYunWrnnOqGUpEJMt06waTJsGQ\nIfDJJ3FHo5qFiEhaGzsWVqyA554LTVRVofUsRESy1L59YUnWLl3g97+v2rnUDCUikqXy8mDqVJg9\nO/wbF9UsREQywNq1cPrpsHAhdO9euXOoZiEikuW6doX77oMLLoCPP67991fNQkQkg4wdC6+9BgsW\nHHyHt2oWIiI54le/gkMOgVtuqd33VbIQEckgeXlhZPcTT4RR3rVFzVAiIhlo1So46yxYuhSOOy61\nY9QMJSKSY7p3hwkTwgjv7dtr/v1UsxARyWDXXAMbNoRxGIck+fmvmoWISI66++4wd9RvflOz76Oa\nhYhIhtuyJazh/T//E/oxKqKahYhIDmvbNkwFMnw4bNxYM++RNFmYWX8zW2dm75jZ2ArKTIxeX21m\n3ZMda2aHmdkiM3vbzBaaWbMyr90clV9nZglypIiIlCgshBtugIsvDivtVbeEycLM8oBJQH+gMzDM\nzDqVKzMA6ODuHYGRwOQUjh0HLHL3o4El0XPMrDNwSVS+P3C/man2k0BRUVHcIaQNXYtSuhalcula\n/PSn0LIljBtX/edO9kXcC1jv7hvcfQ8wHRhcrswgYAqAuy8DmplZqyTHfnVM9O/50ePBwDR33+Pu\nG4D10XmkArn0h5CMrkUpXYtSuXQtzEK/xaxZMGdO9Z47WbJoC2wq83xztC+VMm0SHNvS3bdFj7cB\nLaPHbaJyid5PREQqcNhhMG0ajBwZbqmtLsmSRaq3HKXSu24HOl90W1Oi99FtTyIiB+GEE+Cmm2Do\n0LB4UrVw9wo3oA/wXJnnNwNjy5X5EzC0zPN1hJpChcdGZVpFj1sD66LH44BxZY55Duh9gLhcmzZt\n2rQd/JboOz/RlmyC29eAjmbWDthK6HweVq7MXGAMMN3M+gDb3X2bmX2S4Ni5wOXAb6N/Z5fZP9XM\nJhCanzoCy8sHVdn7hEVEpHISJgt332tmY4AFQB7wkLsXm9mo6PUH3H2emQ0ws/XADmBEomOjU/8G\nmGFmVwEbgIujY940sxnAm8BeYLRG34mIxC8jR3CLiEjtyqgxDKkMEMxWZlZgZkvN7A0ze93Mron2\nVzjAMduZWZ6ZrTKzp6PnOXktzKyZmT1pZsVm9qaZ9c7ha/Ff0d/HWjObamb1c+VamNnDZrbNzNaW\n2VdtA6AzJlmkMkAwy+0B/svduxBuHvhx9PkPOMAxR1xLaLIsqR7n6rX4IzDP3TsB3yHcQJJz18LM\n2gI/AXq4e1dC8/dQcudaPEL4fiyr2gZAZ0yyILUBglnL3T9w9/+LHn8BFBNuAqhogGNWM7N8YADw\nF0pv3c65a2FmTYFT3P1hCH2F7v4vcvBaROoAh5pZHeBQws01OXEt3P0l4NNyu6ttAHQmJYtUBgjm\nhOgOs+7AMioe4Jjt/gDcCOwvsy8Xr0V74CMze8TMVprZg2bWiBy8Fu6+Bfg98HdCktju7ovIwWtR\nRrUNgM6kZKGeeMDMGgMzgWvd/fOyr6UwwDErmNlA4EN3X0UFA0Jz5VoQfkkfD9zv7scT7kj8WjNL\nrlwLM/sW4Zd0O8KXYWMz+37ZMrlyLQ6kqgOgMylZbAEKyjwv4OuZMeuZWV1Covhfdy8Zm7ItmosL\nM2sNfBhXfLXoRGCQmb0PTANON7P/JTevxWZgs7v/LXr+JCF5fJCD16If8L67f+Lue4FZwAnk5rUo\nUdHfRPnv0/xoX4UyKVl8NUDQzOoROmfmxhxTrTEzAx4C3nT3e8q8VDLAEb4+wDFrufst7l7g7u0J\nHZjPu/sPyM1r8QGwycyOjnb1A94AnibHrgWwEehjZg2jv5d+hBsgcvFalKjob2IuMNTM6plZeyoY\nAF1WRo3O/megAAAAsElEQVSzMLNzgHsoHeT365hDqjVmdjLwIrCG0urizYT/wDOAI4gGOLp7LSzf\nnh7MrC/wU3cfZGaHkYPXwsy6ETr66wHvEgbG5pGb12I84YfkXmAlcDXQhBy4FmY2DegLtCD0T/wM\nmEMFn93MbgGuJFyra919QcLzZ1KyEBGReGRSM5SIiMREyUJERJJSshARkaSULEREJCklCxERSUrJ\nQkREklKyEBGRpJQsREQkqf8HcogznVYp2ycAAAAASUVORK5CYII=\n", + "text": [ + "" + ] + } + ], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..6a8491a --- /dev/null +++ b/main.cpp @@ -0,0 +1,106 @@ +#include "tube.h" +#include "tracer.h" + +#include +#include +#include +#include +TRACETHIS; + +using geom::Geom; +using geom::Grid; +namespace po = boost::program_options; +using namespace td; +SPOILNAMESPACE; + +// Directory for storage +const char* dir="results"; + + +inline double min(d x,d y){return x(&f)->default_value(100)) + ; + d T=1/f; + d L=1; + us gp=500; + + d dx=L/(gp-1); // One left and right gp, so + d CFL=0.95; + // This is questionable + d dt=min(CFL*dx/gc.c0(),T/50); + cout << "dt:" << dt <<"\n"; + cout << "dx:" << dx <<"\n"; + + td::gc=tasystem::Globalconf::airSTP(0,f,1); + cout << "omg:" << td::gc.getomg() << endl; + inittrace(loglevel); + initothertrace(loglevel,nltracer); + + + // Create results dir, put all stuff in + mkdir(dir,S_IRWXU | S_IRWXG); + // Empty this directory + + td::Tube t(L,gp); + + + int nsaves_per_period=50; + + int approx_instances_per_period=round(T/dt); + if(nsaves_per_period>approx_instances_per_period){ // Not good.. + WARN("Error: to many saves for number of time instances available. Try adjusting nsaves_per_period Exiting.."); + exit(192); + } + + + d time=0; + + us loopnr=1; + us savenr=0; + us saves_per_period=20; + + d tbetweensaves=T/saves_per_period; + + d tlastsave=0; + + std::stringstream pfilename,ufilename; + pfilename << dir; ufilename << dir; + pfilename << "/p.dat"; + ufilename << "/u.dat"; + std::ofstream p(pfilename.str()); + std::ofstream u(ufilename.str()); + p << "#Pressure data vs gridpoint. gp=" << gp << "\n"; + u << "#Velocity data vs gridpoint. gp=" << gp << "\n"; + + t.getSol().save(p,data::p); + t.getSol().save(u,data::u); + const clock_t begin_time = std::clock(); + while(time<180*T){ + t.DoIntegration(dt); + time=t.getTime(); + + if(time-tlastsave>tbetweensaves){ + tlastsave=time; + savenr++; + t.getSol().save(p,data::p); + t.getSol().save(u,data::u); + cout << "Saving for t= " << time << " s.\n"; + } // ti%10 + loopnr++; + } + d comptime=std::clock()-begin_time; + p << "#Total time required to obtain solution: " << comptime/CLOCKS_PER_SEC <<" \n"; + return 0; +} + + + diff --git a/pmovie.py b/pmovie.py new file mode 100755 index 0000000..82e5474 --- /dev/null +++ b/pmovie.py @@ -0,0 +1,50 @@ +#!/usr/bin/python + +from numpy import * +from matplotlib.pyplot import * +from matplotlib import animation +import sys + + +inst=0 +print("Loading data...") +pdat=loadtxt('results/p.dat') +# save(pdat,'p.npy') + +gp=300 # Need to be taken from file + +xmax=pdat[-1,0] +x=linspace(0,1,300) + +fig = figure(figsize=(12,8)) +ax = axes(xlim=(0,1),ylim=(-180,180)) +grid('on') +line, = plot([],[],lw=2) + +ponly=pdat[:,1] +imax=int(floor(pdat.shape[0]/gp)) + +p=array_split(ponly,imax) + +def init(): + line.set_data([], []) + return line, + + +# print(imax) +def update_progress(progress): + sys.stdout.write('\r[{0}{1}] {2}%'.format('#'*int(floor(progress)),' '*int(ceil((100-progress))),progress)) + +def animate(i): + update_progress(int(round(100*i/imax))) + line.set_data(x,p[i]) + return line, +print("Creating animation...") +anim = animation.FuncAnimation(fig, animate, init_func=init, + frames=imax, interval=20, blit=True) +# print("Saving animation...") +# anim.save('p.mp4', fps=20, extra_args=['-vcodec', 'libx264']) + +print("Showing animation...") +show() + diff --git a/ptime.py b/ptime.py new file mode 100755 index 0000000..02e6920 --- /dev/null +++ b/ptime.py @@ -0,0 +1,13 @@ +#!/usr/bin/python + +from numpy import * +from matplotlib.pyplot import * + +pfile=loadtxt('results/p.dat') +gp=300 +pend=zeros(floor(pfile.shape[0]/gp)) +for i in range(pend.size): + pend[i]=pfile[(i+1)*gp-1,1] + +plot(pend) +show() diff --git a/src/solution.cpp b/src/solution.cpp deleted file mode 100644 index 502b5d0..0000000 --- a/src/solution.cpp +++ /dev/null @@ -1,54 +0,0 @@ -#include "solution.h" -#include "geom.h" -#include "tube.h" - - - -namespace td{ - SolutionAtGp::SolutionAtGp() { - rho_=gc.rho0(); - } - SolutionInstance::SolutionInstance(d time,int gp): - time(time) - { - TRACE(15,"SolutionInstance::SolutionInstance()"); - gps.resize(gp); - } - void SolutionInstance::setrho(d rho) { - TRACE(15,"SolutionInstance::setrho()"); - - for(auto gp=gps.begin();gp!=gps.end();gp++){ - gp->set(rho); - } - } // setrho() - int SolutionInstance::safe(const char* dir){ - TRACE(15,"SolutionInstance::safe()"); - std::stringstream pfilename,ufilename,rhofilename; - pfilename.str(""); ufilename.str(""); rhofilename.str(""); - pfilename << dir; - ufilename << dir; - rhofilename << dir; - - pfilename << "/p" << nr << ".dat"; - ufilename << "/u" << nr << ".dat"; - rhofilename << "/rho" << nr << ".dat"; - cout << pfilename.str() << "\n"; - - std::ofstream p(pfilename.str()); - std::ofstream u(ufilename.str()); - std::ofstream rho(rhofilename.str()); - - // Writing results for pressure to file - p << "#Pressure data vs gridpoint: t=" << time << "\n"; - u << "#Velocity data vs gridpoint: t=" << time << "\n"; - rho << "#Density data vs gridpoint: t=" << time << "\n"; - for(us j=0;j gps; - d time=0; - int nr=0; - public: - SolutionInstance(d TimeInstance,int gp); - ~SolutionInstance(){} - d getTime() const {return time;} - void setTime(d t) {time=t;} - void set(us i,SolutionAtGp& sol){ gps.at(i)=sol;} - SolutionAtGp& get(us i){ return gps.at(i); } - - void setrho(d rho); - int safe(const char* dir); - }; - -} // namespace td -#endif /* _SOLUTION_H_ */ - diff --git a/src/solutionatgp.cpp b/src/solutionatgp.cpp new file mode 100644 index 0000000..ced3dce --- /dev/null +++ b/src/solutionatgp.cpp @@ -0,0 +1,32 @@ +#include "solutionatgp.h" +#include "globalconf.h" + + +TRACETHIS + +namespace td{ + + extern tasystem::Globalconf gc; + + SolutionAtGp::SolutionAtGp() { + rho_=gc.rho0(); + } + d SolutionAtGp::u() const {return m()/rho();} + d SolutionAtGp::ekin() const {return 0.5*rho()*pow(u(),2);} + d SolutionAtGp::estat() const {return rhoE()-ekin();} + d SolutionAtGp::p() const {return estat()*(gc.gas.gamma(gc.T0)-1);} + + d SolutionAtGp::Cflux() const {return m();} + d SolutionAtGp::Mflux() const {return pow(m(),2)/rho()+p();} + d SolutionAtGp::Eflux() const {return (rhoE()+gc.p0/(gc.gas.gamma(gc.T0)-1)+p()+gc.p0)*u();} + +} // namespace td + + + + + + + + + diff --git a/src/solutionatgp.h b/src/solutionatgp.h new file mode 100644 index 0000000..2bb5e37 --- /dev/null +++ b/src/solutionatgp.h @@ -0,0 +1,36 @@ +#pragma once + +namespace td{ + + typedef double d; + class SolutionAtGp { + d rho_,m_=0,rhoE_=0; + public: + // Dependent variables + const d& rho() const {return rho_;} + const d& m() const {return m_;} // which is rho*u + const d& rhoE() const {return rhoE_;} // which is rho*E + // Derived from dependent + d u() const; + d ekin() const; + d estat() const; + d p() const; + + d Cflux() const; + d Mflux() const; + d Eflux() const; + SolutionAtGp(); + ~SolutionAtGp(){} + void setData(d rho,d m=0,d rhoE=0){ + rho_=rho; + m_=m; + rhoE_=rhoE; + } + }; +} // namespace td + + + + + + diff --git a/src/solutioninstance.cpp b/src/solutioninstance.cpp new file mode 100644 index 0000000..0213e8d --- /dev/null +++ b/src/solutioninstance.cpp @@ -0,0 +1,44 @@ +#include "solutioninstance.h" +#include "tube.h" + + + +namespace td{ + + SolutionInstance::SolutionInstance(int gp) + { + TRACE(12,"SolutionInstance::SolutionInstance()"); + gps.resize(gp); + } + void SolutionInstance::setrho(d rho) { + TRACE(15,"SolutionInstance::setrho()"); + + for(auto gp=gps.begin();gp!=gps.end();gp++){ + gp->set(rho); + } + } // setrho() + int SolutionInstance::save(std::ofstream& str,data d){ + TRACE(15,"SolutionInstance::save()"); + if(d==data::p) + for(us j=0;j gps; + d time=0; + public: + SolutionInstance(int gp); + ~SolutionInstance(){} + d getTime() const {return time;} + void setTime(d t) {time=t;} + void set(us i,SolutionAtGp& sol){ gps.at(i)=sol;} + SolutionAtGp& get(us i){ return gps.at(i); } + + void setrho(d rho); + int save(std::ofstream&,data d); + }; + +} // namespace td + diff --git a/src/tube.cpp b/src/tube.cpp index 5183536..3ddaac5 100644 --- a/src/tube.cpp +++ b/src/tube.cpp @@ -12,7 +12,7 @@ namespace td{ return pleft; } - Tube::Tube(d L,int gp): L(L),gp(gp),sol(SolutionInstance(0,gp)) { + Tube::Tube(d L,int gp): L(L),gp(gp),sol(SolutionInstance(gp)) { TRACE(15,"Tube::Tube()"); dx=L/(gp-1); VARTRACE(15,gp); @@ -26,13 +26,11 @@ namespace td{ d newt=t+dt; // Define new solutioninstance - SolutionInstance newsol=sol; - + SolutionInstance newsol(gp); + newsol.setTime(newt); SolutionAtGp newsolgp; d rho,m,rhoE; - newsol.setTime(newt); - int i=0; // Left boundary { @@ -43,16 +41,13 @@ namespace td{ d momfluxl=ip0.rho()*pow(ip0.u(),2)+pleft(t); m=ip0.m()-la*(ip1.Mflux()-momfluxl); rhoE=ip0.rhoE()-la*(ip1.Eflux()-ip0.Eflux()); - // rhoE=ip0.rhoE()-la*(ip1.Eflux()-.rhoE()*ip0.u()-ip0.u()*(gc.p0+pleft(t))); newsolgp.set(rho,m,rhoE); - // cout << "pleft: " << pleft(t) <<"\n"; - // cout << "pleft: " << newsolgp.p() <<"\n"; newsol.set(i,newsolgp); } // Leftmost node // TRACE(15,"leftmost done"); // Inner nodes - for(i=1;isol=sol;} void DoIntegration(d dt); diff --git a/timedomain.cpp b/timedomain.cpp deleted file mode 100644 index 2a0bf40..0000000 --- a/timedomain.cpp +++ /dev/null @@ -1,73 +0,0 @@ -#include "tube.h" -#include "tracer.h" - -#include -#include -#include -TRACETHIS; - -using geom::Geom; -using geom::Grid; -namespace po = boost::program_options; -using namespace td; -SPOILNAMESPACE; - -// Directory for storage -const char* dir="results"; - - - -inline double min(d x,d y){return x(&f)->default_value(100)) - ; - d T=1/f; - d L=0.01; - us gp=10; - - d dx=L/gp; - d CFL=0.5; - // This is questionable - d dt=min(CFL*dx/gc.c0(),T/50); - cout << "dt:" << dt <<"\n"; - cout << "dx:" << dx <<"\n"; - - td::gc=tasystem::Globalconf::airSTP(0,f,1); - cout << "omg:" << td::gc.getomg() << endl; - inittrace(15); - initothertrace(10,nltracer); - - - // Create results dir, put all stuff in - mkdir(dir,S_IRWXU | S_IRWXG); - // Empty this directory - system("rm results/*"); - - td::Tube t(L,gp); - - - int nsaves_per_period=50; - int nr=0; - d time=0; - - while(time<2*T){ - t.DoIntegration(dt); - time=t.getTime(); - // cout << "Iteration done. Saving data...\n"; - if(true){ - // Create filename to store data to - } // ti%10 - } - return 0; -} - - - diff --git a/timedomain.i b/timedomain.i new file mode 100644 index 0000000..6383816 --- /dev/null +++ b/timedomain.i @@ -0,0 +1,52 @@ +%module timedomaineuler +%{ + #include "solutionatgp.h" + #include "test.h" + #include "tube.h" +%} + +typedef double d; + +namespace td{ + class SolutionAtGp { + d rho_,m_=0,rhoE_=0; + public: + // Dependent variables + const d& rho() const; + const d& m() const; + const d& rhoE() const; + // Derived from dependent + d u() const; + d ekin() const; + d estat() const; + d p() const; + + d Cflux() const; + d Mflux() const; + d Eflux() const; + SolutionAtGp(); + ~SolutionAtGp(){} + void setData(d rho,d m=0,d rhoE=0); + }; + class SolutionInstance{ + public: + SolutionInstance(int gp); + ~SolutionInstance(){} + d getTime() const; + void setTime(d t); + void set(us i,SolutionAtGp& sol); + SolutionAtGp& get(us i); + void setrho(d rho); + }; + + class Tube{ + public: + Tube(double L,int gp); + ~Tube(){} + SolutionInstance& getSol() { return sol;} + void setSol(const SolutionInstance& sol) {this->sol=sol;} + void DoIntegration(d dt); + d getTime(){return t;} + }; +} +