Vectorized implementation. Does not work yet

CMakeLists.txt

@@ -11,6 +11,9 @@ INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
 
 SET(CMAKE_SWIG_FLAGS "-Wall")
 
+add_definitions(-DNDEBUG)	#This is for BOOST functions, and to disable assertions!
+# add_definitions(-DTRACER=0)
+add_definitions(-DEIGEN_NO_DEBUG -DARMA_NO_DEBUG)
 
 add_definitions(-DTRACERNAME=timedomaineulertracer -DARMA_USE_BLAS -DARMA_USE_LAPACK)
 

src/arma_numpy.cpp

@@ -1,3 +1,10 @@
+// arma_numpy.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
 #define PY_ARRAY_UNIQUE_SYMBOL npy_array_api
 #define NO_IMPORT_ARRAY
 #include "arma_numpy.h"
@@ -5,34 +12,30 @@
 // Annes conversion functions. In a later stage they have to be
 // generalized for arrays of arbitrary dimensions
 
+PyObject *npy_from_vd(const vd &in) {
+  long int size = in.size();
+  npy_intp dims[1] = {size};
 
-
-PyObject* npy_from_vd(const vd& in){
-  long int size=in.size();
-  npy_intp dims[1]={size};
   // This code should be improved massively?
-  if(size==0){
+  if (size == 0) {
     std::cout << "Size is zero!\n";
     return NULL;
   }
-  cout << "SFSF\n";    
-  PyArrayObject* array =  (PyArrayObject*) \
-          PyArray_SimpleNew(1, dims, NPY_DOUBLE); 
-  cout << "SFSF\n";    
+  PyArrayObject *array =
+      (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
   if (array == NULL) {
     return NULL;
   }
-  double *pydat= (double*) PyArray_DATA(array);
-  cout << "SFSF\n";
-  mempcpy(pydat,in.memptr(),size*sizeof(double));
+  double *pydat = (double *)PyArray_DATA(array);
+  mempcpy(pydat, in.memptr(), size * sizeof(double));
+
   return PyArray_Return(array);
 }
-vd vd_from_npy(const PyArrayObject* const in){
+vd vd_from_npy(const PyArrayObject *const in) {
 
-  double* pydata = (double*) PyArray_DATA(in);
-  vd result;
-  // long int dims=
-  // PyArra
-  //   vd res()
+  npy_intp size=PyArray_DIMS(in)[0]; // Extract first 
+  double *pydata = (double *)PyArray_DATA(in);
+  vd result(size);
+  memcpy(result.memptr(),pydata,size*sizeof(double));
   return result;
 }

src/arma_numpy.h

@@ -1,13 +1,25 @@
+// arma_numpy.h
+//
+// Author: J.A. de Jong 
+//
+// Description: Header file for two function converting between Numpy
+// arrays and armadillo matrices
+//////////////////////////////////////////////////////////////////////
 #pragma once
 #ifndef ARMA_NUMPY_H
-#define ARMA_NUMPY_H
+#define ARMA_NUMPY_H 1
+
 #include <Python.h>
 #include <numpy/arrayobject.h>
 #include "vtypes.h"
+
 SPOILNAMESPACE
 vd vd_from_npy(PyArrayObject const * const);
 PyObject* npy_from_vd(const vd& armavec);
 
 
 
+
 #endif // ARMA_NUMPY_H
+//////////////////////////////////////////////////////////////////////
+

src/solutionatgp.cpp

@@ -1,31 +0,0 @@
-#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

@@ -1,38 +0,0 @@
-#pragma once
-
-namespace td{
-
-  typedef double d;
-  typedef unsigned us;
-
-  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

@@ -2,38 +2,18 @@
 
 namespace td{
 
-  SolutionInstance::SolutionInstance(int gp)
+  SolutionInstance::SolutionInstance(int gp,d rho)
   {
     TRACE(12,"SolutionInstance::SolutionInstance()");
-    gps.resize(gp);
+    rho_=vd(gp,fillwith::zeros);
+    m_=vd(gp,fillwith::zeros);
+    rhoE_=vd(gp,fillwith::zeros);
+
   }
   void SolutionInstance::setrho(d rho)  {
     TRACE(15,"SolutionInstance::setrho()");
-
-    for(auto gp=gps.begin();gp!=gps.end();gp++){
-      gp->setData(rho);
-    }
+    rho_=ones(rho_.size())*rho;
   } // setrho()
-  vd SolutionInstance::getrho() const{
-    return gpsol_to_python(gps,&SolutionAtGp::rho);
-  }
-  vd SolutionInstance::getp() const{
-    return gpsol_to_python(gps,&SolutionAtGp::p);
-  }
-  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

@@ -1,42 +1,53 @@
+// solutioninstance.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
 #pragma once
-#include "solutionatgp.h"
-#include "vtypes.h"
+#ifndef SOLUTIONINSTANCE_H
+#define SOLUTIONINSTANCE_H 1
 
-void do_arma_to_python(const vd& armavec,double** data,int* n);
+#include "vtypes.h"
+#include "globalconf.h"
 
 namespace td{
   SPOILNAMESPACE
-
-  template<typename funptr>
-  vd gpsol_to_python(const vector<SolutionAtGp>& gps,funptr fun){
-    int size=gps.size();
-    vd result(size);
-    for(int i=0;i<size;i++)
-      result(i)=((gps[i]).*fun)();
-    return result;
-  }
-
-
-  enum data{
-    p,u
-  };
+  extern tasystem::Globalconf gc;
 
   class SolutionInstance{
-    vector<SolutionAtGp> gps;
+
+    vd rho_,m_,rhoE_;
     d time=0;
   public:
-    SolutionInstance(int gp);
+    SolutionInstance(int gp,d rho=1.2);
     ~SolutionInstance(){}
     d getTime() const {return time;}
     void setTime(d t) {time=t;}
-    void setData(us i,SolutionAtGp& sol){ gps.at(i)=sol;}
-    SolutionAtGp& get(us i){ return gps.at(i); }
 
-    vd getp() const;
-    vd getrho() const;
+    vd& rho_ref(){ return rho_;}
+    vd& m_ref() {return m_;}
+    vd& rhoE_ref() {return rhoE_;}    
+    const vd& rho() const {return rho_;}
+    const vd& m() const {return m_;}
+    const vd& rhoE() const {return rhoE_;}
+    vd p() const {return estat()*(gc.gas.gamma(gc.T0)-1);}
+    vd u() const {return m()/rho();}
+    vd ekin() const { return 0.5*rho()%pow(u(),2);}
+    vd estat() const {return rhoE()-ekin();}
+    vd Cflux() const {return m();}
+    vd Mflux() const {return pow(m(),2)/rho()+p();}
+    vd Eflux() const {return (rhoE()+gc.p0/(gc.gas.gamma(gc.T0)-1)+p()+gc.p0)%u();}
+
     void setrho(d rho);
-    int save(std::ofstream&,data d);
+
   };
 
 } // namespace td
 
+
+
+#endif // SOLUTIONINSTANCE_H
+//////////////////////////////////////////////////////////////////////
+

src/tube.cpp

@@ -1,91 +1,92 @@
 #include "tube.h"
 
-namespace td{
-  tasystem::Globalconf gc;
+TRACETHIS
 
-  d Tube::pleft(d t){
-    d pleft=std::sin(gc.getomg()*t);
-    cout << "t:"<< t<<endl;
-    cout << "omg:" << gc.getomg()<< endl;
-    // cout << "pleft:"<< pleft<<endl;
-    return pleft;
+namespace td {
+tasystem::Globalconf gc;
 
+d Tube::pleft(d t) {
+  d pleft = std::sin(gc.getomg() * t);
+  // cout << "t:" << t << endl;
+  // cout << "omg:" << gc.getomg() << endl;
+  cout << "pleft:"<< pleft<<endl;
+  return pleft;
+}
+  Tube::Tube(d L, int gp) throw(int) : L(L), gp(gp), sol(SolutionInstance(gp)) {
+  TRACE(15, "Tube::Tube()");
+  dx = L / (gp - 1);
+  if(gp<5)
+    throw 1;
+  VARTRACE(15, gp);
+  sol.setrho(gc.rho0());
+}
+
+void Tube::DoIntegration(d dt, int n) {
+  TRACE(14, "Tube::DoIntegration()");
+  int integrationnumber = 0;
+  while (integrationnumber < n) {
+    Integrate(dt);
+    integrationnumber++;
   }
-  Tube::Tube(d L,int gp):
-    L(L),
-    gp(gp),
-    sol(SolutionInstance(gp))  {
+}
+void TubeLF::Integrate(d dt) {
+  // Integrate using Lax-Friedrich method
+  d newt = t + dt;
+  // Define new solutioninstance
+  SolutionInstance newsol(gp);
+  newsol.setTime(newt);
 
-    TRACE(15,"Tube::Tube()");
-    dx=L/(gp-1);
-    VARTRACE(15,gp);
-    sol.setrho(gc.rho0());
-  }
+  const vd& oldrho=sol.rho();
+  const vd& oldm=sol.m();
+  const vd& oldrhoE=sol.rhoE();
 
-  void Tube::DoIntegration(d dt,int n){
-    TRACE(14,"Tube::DoIntegration()");
-    int integrations=0;
-    while(integrations<n){
-      Integrate(dt);
-      integrations++;
-    }
-  }
-  void Tube::Integrate(d dt){
+  vd& rho=newsol.rho_ref();
+  vd& m=newsol.m_ref();
+  vd& rhoE=newsol.rhoE_ref();
 
-    // Define new time
-    d newt=t+dt;
+  // Fluxes from previous solution
+  vd Cflux=sol.Cflux();
+  vd Mflux=sol.Mflux();
+  vd Eflux=sol.Eflux();
 
-    // Define new solutioninstance
-    SolutionInstance newsol(gp);
-    newsol.setTime(newt);
-    SolutionAtGp newsolgp;
-    d rho,m,rhoE;
+  // std::cout << "poL" << sol.p()(0) << std::endl;
+  {  // Left boundary
+    d la = dt/dx;
+    rho(0) = oldrho(0);
+    rho(0)+=-la*(Cflux(1)-Cflux(0));
+    d oldu0=oldm(0)/oldrho(0);
+    // std::cout << "oldu:"<< oldu0 << std::endl;
+    d momfluxl = oldrho(0)*pow(oldu0, 2) + pleft(t);
+    m(0) = oldm(0);
+    m(0)+=-la*(Mflux(1)-momfluxl);
+    rhoE(0) = oldrhoE(0);
+    rhoE(0)+=-la*(Eflux(1)-Eflux(0));
 
-    int i=0;
-    // Left boundary
-    {
-      SolutionAtGp& ip1=sol.get(i+1);
-      SolutionAtGp& ip0=sol.get(i);
-      d la=dt/(dx);
-      rho=ip0.rho()-la*(ip1.Cflux()-ip0.Cflux());
-      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());
-      newsolgp.setData(rho,m,rhoE);
-      newsol.setData(i,newsolgp);
+  } //  End left boundary
+  {				// Inner nodes
+    d lambda = dt/(2*dx);
+    rho.subvec(1,gp-2)=0.5*(oldrho.subvec(0,gp-3) + oldrho.subvec(2,gp-1));
+    rho.subvec(1,gp-2)+=-lambda*(Cflux.subvec(0,gp-3) -Cflux.subvec(2,gp-1));
+  
+    m.subvec(1,gp-2)=0.5*(oldm.subvec(0,gp-3) + oldm.subvec(2,gp-1));
+    m.subvec(1,gp-2)+=-lambda*(Mflux.subvec(0,gp-3) -Mflux.subvec(2,gp-1));
 
-    } // Leftmost node
-    // TRACE(15,"leftmost done");
-    // Inner nodes
-    for(i=1;i<gp-1;i++){
-      SolutionAtGp& im1=sol.get(i-1);
-      SolutionAtGp& ip1=sol.get(i+1);
-      d lambda=dt/(2*dx);
-      rho=0.5*(im1.rho()+ip1.rho())-lambda*(ip1.Cflux()-im1.Cflux());
-      m=0.5*(im1.m()+ip1.m())-lambda*(ip1.Mflux()-im1.Mflux());
-      rhoE=0.5*(im1.rhoE()+ip1.rhoE())-lambda*(ip1.Eflux()-im1.Eflux());
+    rhoE.subvec(1,gp-2)=0.5*(oldrhoE.subvec(0,gp-3) + oldrhoE.subvec(2,gp-1));
+    rhoE.subvec(1,gp-2)+=-lambda*(Eflux.subvec(0,gp-3) -Eflux.subvec(2,gp-1));
+  } // End inner nodes
+  {				// Right boundary
+    int i = gp - 1;
+    d la = dt/dx;
+    rho(i) = oldrho(i);
+    rho(i)+=-la*(0-Cflux(i-1));
+    m(i) = 0;
+    rhoE(i) = oldrhoE(i);
+    rhoE(i)+=-la*(0-Eflux(i-1));
+  } // End right boundary
 
-      newsolgp.setData(rho,m,rhoE);
-      newsol.setData(i,newsolgp);
-      
-    } // for over all gridpoints in mid
-    {
-      i=gp-1;
-      SolutionAtGp& im1=sol.get(i-1);
-      SolutionAtGp& ip0=sol.get(i);
-
-      d la=dt/(dx);
-      rho=ip0.rho()-la*(0-im1.Cflux());
-      m=0;
-      rhoE=ip0.rhoE()-la*(0-im1.Eflux());
-
-      newsolgp.setData(rho,m,rhoE);
-      newsol.setData(i,newsolgp);
-    }
-
-    // Finally, update time and solution
-    sol=newsol;
-    t=newt;
-  }
+  // Finally, update time and solution
+  sol = newsol;
+  t = newt;
+}
 
 } // namespace td

src/tube.h

@@ -1,29 +1,52 @@
-#pragma once 
-#ifndef _TUBE_H_
-#define _TUBE_H_
+// tube.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TUBE_H
+#define TUBE_H 1
 
 #include "solutioninstance.h"
 #include "globalconf.h"
 
-namespace td{
+namespace td {
 
   typedef double d;
 
-  class Tube{
-    d dx,L;                     // Grid spacing, total length
-    int gp;                     // Number of gridpoints
+  class Tube {
+  protected:
+    d dx, L;              // Grid spacing, total length
+    int gp;               // Number of gridpoints
     SolutionInstance sol; // Solutions at time instances
-    d t=0;                // Current time
+    d t = 0;              // Current time
     d pleft(d t);         // Compute pressure bc
-    void Integrate(d dt);
-  public:
-    Tube(double L,int gp);
-    ~Tube(){}
-    SolutionInstance& getSol() { return sol;}
-    void setSol(const SolutionInstance& sol) {this->sol=sol;}
-    void DoIntegration(d dt,int n=1);
-    d getTime(){return t;}
-  };
+    virtual void Integrate(d dt) = 0;
+
+  public:
+    Tube(double L, int gp) throw(int);
+    virtual ~Tube() {}
+    SolutionInstance &getSol() { return sol; }
+    void setSol(const SolutionInstance &sol) { this->sol = sol; }
+    void DoIntegration(d dt, int n = 1);
+    d getTime() { return t; }
+  };
+  class TubeLF : public Tube { // Using Lax-Friedrichs method
+  protected:
+    virtual void Integrate(d dt);
+
+  public:
+    using Tube::Tube;
+
+  };
+} // namespace td
+
+
+
+#endif // TUBE_H
+//////////////////////////////////////////////////////////////////////
+
+
 
-}                  // namespace td
-#endif /* _TUBE_H_ */

timedomain.i

@@ -35,13 +35,12 @@ class vd{
     return NULL;
   int ndims=PyArray_DIM(arr,0);
   if(ndims!=1){
-    PyErr_SetString(PyExc_TypeError,"Number of dimensions not equal to 1");
+    // PyErr_SetString(PyExc_TypeError,"Number of dimensions not equal to 1");
     return NULL;
   }
   $1= vd_from_npy(arr);
 }
 %typemap(out) vd {
-  cout << "Hoi\n";
   $result=npy_from_vd($1);
 }
 typedef double d;
@@ -95,46 +94,33 @@ namespace tasystem{
 
 namespace td{
 
-  class SolutionAtGp {
-  public:
-    const d& rho() const;
-    const d& m() const;
-    const d& rhoE() const;
-    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(int gp,d rho=1.2);
     ~SolutionInstance(){}
-    d getTime() const;
-    vd getrho() const;
-    vd getp() const;
+    vd rho() const;
+    vd p() const;
+    vd u() const;
     void setTime(d t);
-    void setData(us i,SolutionAtGp& sol);
-    SolutionAtGp& get(us i);
     void setrho(d rho);
   };
 
   tasystem::Globalconf gc;
 
   class Tube{
+    virtual void Integrate(d dt)=0;
   public:
-    Tube(double L,int gp);
-    ~Tube(){}
+    Tube(double L,int gp) throw(int);
+    virtual ~Tube();
     SolutionInstance& getSol();
     void setSol(const SolutionInstance& sol);
     void DoIntegration(d dt,int n=1);
     d getTime();
   };
+  class TubeLF:public Tube{
+    virtual void Integrate(d dt);
+  public:
+    TubeLF(d L,int gp);
+  };
 }