We give Swig a try

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})
 
 
 

main.cpp

@@ -0,0 +1,106 @@
+#include "tube.h"
+#include "tracer.h"
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <boost/program_options.hpp>
+#include <ctime>
+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<y?x:y;}
+
+int main(int argc,char *argv[]){
+
+  d f=85.785;
+  int loglevel=20;
+
+  po::options_description desc("Allowed options");
+  desc.add_options()
+    ("help", "Produce help message")
+    ("frequency",po::value<d>(&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;
+}
+
+
+

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()
+

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()

src/solution.cpp

@@ -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.size();j++){
-      SolutionAtGp& solgp=get(j);
-      p << j << "\t" <<solgp.p() << "\n";
-      u << j << "\t" <<solgp.u() << "\n";
-      rho << j << "\t" <<solgp.rho() << "\n";        
-    }
-    return 0;
-  }
-}   // namespace td
-

src/solution.h

@@ -1,61 +0,0 @@
-#pragma once
-#ifndef _SOLUTION_H_
-#define _SOLUTION_H_
-
-
-#include "vtypes.h"
-#include "geom.h"
-#include "globalconf.h"
-
-namespace td{
-  SPOILNAMESPACE
-
-  extern tasystem::Globalconf gc;
-
-  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 {return m()/rho();}
-    d ekin() const {return 0.5*rho()*pow(u(),2);}
-    d estat() const {return rhoE()-ekin();}
-    d p() const {return estat()*(gc.gas.gamma(gc.T0)-1);}
-
-    d Cflux() const {return m();}
-    d Mflux() const {return pow(m(),2)/rho()+p();}
-    d Eflux() const {return (rhoE()+gc.p0/(gc.gas.gamma(gc.T0)-1)+p()+gc.p0)*u();}
-    SolutionAtGp();
-    ~SolutionAtGp(){}
-    void set(d rho,d m=0,d rhoE=0){
-      rho_=rho;
-      m_=m;
-      rhoE_=rhoE;
-    }
-
-  };
-  
-
-  class SolutionInstance{
-    vector<SolutionAtGp> 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_ */
-

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
+
+
+
+
+
+
+
+
+

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
+
+
+
+
+
+

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.size();j++){
+        SolutionAtGp& solgp=get(j);
+        str << j << "\t" <<solgp.p() << "\n";
+      }
+    else if(d==data::u)
+      for(us j=0;j<gps.size();j++){
+        SolutionAtGp& solgp=get(j);
+        str << j << "\t" <<solgp.u() << "\n";
+      }
+    return 0;
+  }
+
+}   // namespace td
+
+
+
+
+
+
+
+
+

src/solutioninstance.h

@@ -0,0 +1,30 @@
+#pragma once
+#include "solutionatgp.h"
+#include "globalconf.h"
+
+namespace td{
+  SPOILNAMESPACE
+
+  enum data{
+    p,u
+  };
+
+  extern tasystem::Globalconf gc;
+
+  class SolutionInstance{
+    vector<SolutionAtGp> 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
+

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;i<gp-2;i++){
+    for(i=1;i<gp-1;i++){
       SolutionAtGp& im1=sol.get(i-1);
       SolutionAtGp& ip1=sol.get(i+1);
       d lambda=dt/(2*dx);
@@ -65,7 +60,7 @@ namespace td{
       
     } // for over all gridpoints in mid
     {
-      i++;
+      i=gp-1;
       SolutionAtGp& im1=sol.get(i-1);
       SolutionAtGp& ip0=sol.get(i);
 
@@ -83,7 +78,4 @@ namespace td{
     t=newt;
   }
 
-
-
-
 } // namespace td

src/tube.h

@@ -2,15 +2,14 @@
 #ifndef _TUBE_H_
 #define _TUBE_H_
 
-#include "geom.h"
-#include "globalconf.h"
-#include "solution.h"
-
+#include "solutioninstance.h"
 
 namespace td{
 
   extern tasystem::Globalconf gc;
 
+  typedef double d;
+
   class Tube{
     d dx,L;                     // Grid spacing, total length
     int gp;                     // Number of gridpoints
@@ -19,7 +18,8 @@ namespace td{
     d pleft(d t);         // Compute pressure bc
     
   public:
-    Tube(d L,int gp);
+    Tube(double L,int gp);
+    ~Tube(){}
     SolutionInstance& getSol() { return sol;}
     void setSol(const SolutionInstance& sol) {this->sol=sol;}
     void DoIntegration(d dt);

timedomain.cpp

@@ -1,73 +0,0 @@
-#include "tube.h"
-#include "tracer.h"
-
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <boost/program_options.hpp>
-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<y?x:y;}
-
-int main(int argc,char *argv[]){
-
-  d f=100;
-
-
-  po::options_description desc("Allowed options");
-  desc.add_options()
-    ("help", "Produce help message")
-    ("frequency",po::value<d>(&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;
-}
-
-
-

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;}
+  };
+}
+