Improved interface

CMakeLists.txt

@@ -52,6 +52,7 @@ link_directories(/home/anne/bin/lib)
 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 "-includeall")
 swig_add_module(timedomaineuler python timedomain.i)
 swig_link_libraries(timedomaineuler sources nonlin armadillo ${PYTHON_LIBRARIES})
 

gc.i

@@ -0,0 +1,48 @@
+%module gc
+%{
+  #include "globalconf.h"
+%}
+namespace tasystem{
+
+  class Globalconf{
+  public:
+    d T0,p0;			/* Reference temperature and pressure (used to initialize a lot of variables. */
+				// finite volume size, speed of sound,
+				// deltax of volume
+    d kappa;			// Artificial viscosity tuning factor,
+    // typically between 0.25 and 0.75
+
+    gases::Gas gas;
+
+    Globalconf(us Nf=0,d freq=100,const string& gasstring="air",d T0=293.15,d p0=101325.0,d kappa=1.0,bool driven=true);
+    const us& Nf() const;
+    const us& Ns() const;
+    static Globalconf airSTP(us Nf,d freq,d kappa=1.0);
+    ~Globalconf(){TRACE(-5,"~Globalconf()");}
+    bool isDriven() const;
+    void setDriven(bool d);
+    d getomg() const;
+    d getfreq() const;
+    d c0() const;
+    d rho0() const;
+    d deltanu0() const;
+    vd omgvec;    
+    void setNf(us);
+    /* %rename set setNffreq */
+    /* void set(us Nf,d freq);	// Set data for new frequency and */
+    // number of samples
+    void setomg(d omg);
+    void setfreq(d freq);
+
+    void setGas(const string& mat);
+    const string& getGas() const;
+
+    void setMass(d mass);
+    d getMass() const;
+
+    void show() const;
+    //    void setgas(string g){ gas(g);}
+
+  }; /* Class Globalconf */
+
+} /* namespace tasystem */

run.py

@@ -1,35 +1,56 @@
 #!/usr/bin/python
 
+if(!os.path.isfile('timedomaineuler.py')):
+    os.system('cmake . && make')
 from timedomaineuler import *
+
 from numpy import *
+import os
 import sys
-#Globalconf accessable with cvar.gc
-f=85.785                                  #Frequency of oscillation
+from argparse import ArgumentParser
+import argparse
 
+parser=ArgumentParser(description="Compute response to tube\
+with sinusoidal pressure boundary condition",\
+ formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+parser.add_argument('scheme',help="Scheme to use, either LF (Lax-Friedrichs) or\
+ (MCM) MacCormack",choices=["LF","MCM"])
+parser.add_argument('-f','--freq',type=float,dest="f",default=85.785,\
+                  help="Frequency of oscillation")
+parser.add_argument('-L','--length',type=float,dest="L",default=1.0,\
+                  help="Tube length")
+parser.add_argument('-g','--gridpoints',type=int,dest="gp",default=300,\
+                  help="Number of gridpoints")
+parser.add_argument('-c','--CFL',type=float,dest="CFL",default=0.5,\
+                  help="CFL Nummber")
+parser.add_argument('-n','--saves_per_period',type=int,dest="nr_p_period",
+                  default=50,help="Number of saves per period")
+parser.add_argument('-p','--periods',type=float,dest="periods",
+                  default=10,help="Periods to compute")
+args=parser.parse_args()
 
-L=1.                                      #Length of tube
-gp=300                                    #Number of gridpoints
+globals().update(vars(args))    #Put options in global namespace
+print("Frequency:",f)
 
-nr_p_period=50                  #Number of save per oscillation period
-CFL=0.5;                        # CFL number
-periods=100                      #Number of periods to compute
 
 T=1/f
 gc=cvar.gc #Reference!
-dx=L/(gp-1);                    # One left and right gp, so
-
-gc.setfreq(f)
-tube=TubeLF(L,gp)
+dx=L/(gp-1);                    # Grid spacing
+gc.setfreq(f)                   # Set frequency in gc (for c++)
+x=linspace(0,L,gp)              # Grid. Not really needed
+if scheme=="LF":
+    tube=TubeLF(L,gp)
+else:    
+    tube=TubeMCM(L,gp)
 dt=min(CFL*dx/gc.c0(),T/nr_p_period)
 
-intsteps=int(floor(1./(gc.getfreq()*dt)/nr_p_period))
-x=linspace(0,L,gp)
+intsteps=int(floor(1./(gc.getfreq()*dt)/nr_p_period))   #Number of
+                                                        #steps per save
 
 # Create tube instance
 tube=TubeLF(L,gp)
 
 # To create a nice progress bar
-
 def update_progress(progress):
         nbars=60
         text1="\r[{0}{1}] {2}%".format('#'*int(floor(nbars*progress)),' '*int(nbars-floor(nbars*progress)),"%0.0f" %(100*progress))
@@ -37,7 +58,6 @@ def update_progress(progress):
         sys.stdout.write(text1)
         sys.stdout.flush()
 if(__name__=="__main__"):
-
     # Create data storage
     uarr=zeros((nr_p_period*periods+1,gp),dtype=float)
     parr=zeros((nr_p_period*periods+1,gp),dtype=float)

src/solutioninstance.h

@@ -17,7 +17,6 @@ namespace td{
   extern tasystem::Globalconf gc;
 
   class SolutionInstance{
-
     vd rho_,m_,rhoE_;
     d time=0;
   public:
@@ -29,9 +28,11 @@ namespace td{
     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);}

src/tube.cpp

@@ -21,17 +21,18 @@ d Tube::pleft(d t) {
   sol.setrho(gc.rho0());
 }
 
-void Tube::DoIntegration(d dt, int n) {
+void Tube::DoIntegration(d dt, int maxnr) {
   TRACE(14, "Tube::DoIntegration()");
-  int integrationnumber = 0;
-  while (integrationnumber < n) {
-    Integrate(dt);
-    integrationnumber++;
+  int intnr = 0;
+  for (intnr=0; intnr < maxnr; intnr++) {
+    sol=Integrate(dt);		// Update solution  
   }
 }
-void TubeLF::Integrate(d dt) {
+SolutionInstance TubeLF::Integrate(d dt) {
   // Integrate using Lax-Friedrich method
+  d t=sol.getTime();
   d newt = t + dt;
+
   // Define new solutioninstance
   SolutionInstance newsol(gp);
   newsol.setTime(newt);
@@ -84,9 +85,57 @@ void TubeLF::Integrate(d dt) {
     rhoE(i)+=-la*(0-Eflux(i-1));
   } // End right boundary
 
-  // Finally, update time and solution
-  sol = newsol;
-  t = newt;
+  return newsol;
 }
+  SolutionInstance TubeMCM::Integrate(d dt){ // Integrate using
+					     // MacCormack method
+    d t=sol.getTime();
+    d newt=t+dt;
 
+    d lambda=dt/dx;
+
+    SolutionInstance predictor(gp);
+    {				// Predictor step
+      const vd& oldrho=sol.rho();
+      const vd& oldm=sol.m();
+      const vd& oldrhoE=sol.rhoE();
+    
+      vd& prho=predictor.rho_ref();
+      vd& pm=predictor.m_ref();
+      vd& prhoE=predictor.rhoE_ref();
+      // Fluxes from previous solution
+      vd Cflux=sol.Cflux();
+      vd Mflux=sol.Mflux();
+      vd Eflux=sol.Eflux();
+
+      // Density
+      prho.tail(gp-1)=oldrho.tail(gp-1)
+	-lambda*(Cflux.tail(gp-1)-Cflux.head(gp-1));
+      // The first element
+      prho(0)=oldrho(0)-lambda*(Cflux(1)-Cflux(0));
+
+      // Momentum
+      d oldu0=oldm(0)/oldrho(0);
+      d momfluxl = pow(oldm(0), 2)/oldrho(0) + pleft(t);
+
+      pm.tail(gp-1)=oldm.tail(gp-1)
+	-lambda*(Mflux.tail(gp-1)-Mflux.head(gp-1));
+      // The first element
+      pm(0)=oldm(0)-lambda*(Mflux(1)-momfluxl);
+      pm(gp-1)=0;
+
+      // Energy
+      prhoE.tail(gp-1)=oldrhoE.tail(gp-1)
+	-lambda*(Eflux.tail(gp-1)-Eflux.head(gp-1));
+      // The first element
+      prhoE(0)=oldrhoE(0)-lambda*(Eflux(1)-Eflux(0));
+
+    }
+    SolutionInstance corrector(gp);
+    // Temp test hack:
+    corrector=predictor;
+    corrector.setTime(newt);
+
+    return corrector;
+  }					     // Integrate(dt)
 } // namespace td

src/tube.h

@@ -21,9 +21,8 @@ namespace td {
     d dx, L;              // Grid spacing, total length
     int gp;               // Number of gridpoints
     SolutionInstance sol; // Solutions at time instances
-    d t = 0;              // Current time
     d pleft(d t);         // Compute pressure bc
-    virtual void Integrate(d dt) = 0;
+    virtual SolutionInstance Integrate(d dt) = 0;
 
   public:
     Tube(double L, int gp) throw(int);
@@ -31,12 +30,18 @@ namespace td {
     SolutionInstance &getSol() { return sol; }
     void setSol(const SolutionInstance &sol) { this->sol = sol; }
     void DoIntegration(d dt, int n = 1);
-    d getTime() { return t; }
+    d getTime() { return sol.getTime(); }
   };
   class TubeLF : public Tube { // Using Lax-Friedrichs method
   protected:
-    virtual void Integrate(d dt);
+    virtual SolutionInstance Integrate(d dt);
+  public:
+    using Tube::Tube;
 
+  };
+  class TubeMCM : public Tube { // Using Lax-Friedrichs method
+  protected:
+    virtual SolutionInstance Integrate(d dt);
   public:
     using Tube::Tube;
 
@@ -44,7 +49,6 @@ namespace td {
 } // namespace td
 
 
-
 #endif // TUBE_H
 //////////////////////////////////////////////////////////////////////
 

timedomain.i

@@ -9,7 +9,7 @@
   #include "arma_numpy.h"
   #include "tube.h"
 
-  #include "globalconf.h"
+
   SPOILNAMESPACE
     namespace td{
     extern tasystem::Globalconf gc;
@@ -50,54 +50,12 @@ class vd{
 typedef double d;
 typedef unsigned us;
 
-namespace tasystem{
-
-  class Globalconf{
-  public:
-    d T0,p0;			/* Reference temperature and pressure (used to initialize a lot of variables. */
-				// finite volume size, speed of sound,
-				// deltax of volume
-    d kappa;			// Artificial viscosity tuning factor,
-    // typically between 0.25 and 0.75
-
-    gases::Gas gas;
-
-    Globalconf(us Nf=0,d freq=100,const string& gasstring="air",d T0=293.15,d p0=101325.0,d kappa=1.0,bool driven=true);
-    const us& Nf() const;
-    const us& Ns() const;
-    static Globalconf airSTP(us Nf,d freq,d kappa=1.0);
-    ~Globalconf(){TRACE(-5,"~Globalconf()");}
-    bool isDriven() const;
-    void setDriven(bool d);
-    d getomg() const;
-    d getfreq() const;
-    d c0() const;
-    d rho0() const;
-    d deltanu0() const;
-    vd omgvec;    
-    void setNf(us);
-    /* %rename set setNffreq */
-    /* void set(us Nf,d freq);	// Set data for new frequency and */
-    // number of samples
-    void setomg(d omg);
-    void setfreq(d freq);
-
-    void setGas(const string& mat);
-    const string& getGas() const;
-
-    void setMass(d mass);
-    d getMass() const;
-
-    void show() const;
-    //    void setgas(string g){ gas(g);}
-
-  }; /* Class Globalconf */
-
-} /* namespace tasystem */
-
+%include "gc.i"
 
 namespace td{
 
+  tasystem::Globalconf gc;
+
   class SolutionInstance{
   public:
     SolutionInstance(int gp,d rho=1.2);
@@ -111,18 +69,10 @@ namespace td{
     void setTime(d t);
     void setrho(d rho);
   };
-
-  tasystem::Globalconf gc;
-
-  class Tube {
-  public:
-    d dx, L;              // Grid spacing, total length
-    int gp;               // Number of gridpoints
-    SolutionInstance sol; // Solutions at time instances
-    d t = 0;              // Current time
-    d pleft(d t);         // Compute pressure bc
-    virtual void Integrate(d dt) = 0;
-
+  // %feature("abstract") Tube;	// Define Tube abstract
+ class Tube {
+ protected:
+    virtual SolutionInstance Integrate(d dt)=0;
   public:
     Tube(double L, int gp) throw(int);
     virtual ~Tube() {}
@@ -132,9 +82,16 @@ namespace td{
     d getTime() { return t; }
   };
   class TubeLF:public Tube{
-    virtual void Integrate(d dt);
+  protected:
+    virtual SolutionInstance Integrate(d dt);
   public:
     TubeLF(d L,int gp);
   };
+  class TubeMCM:public Tube{
+  protected:
+    virtual SolutionInstance Integrate(d dt);
+  public:
+    TubeMCM(d L,int gp);
+  };
 }