From ca73a502f6f99f1d8dc4dc6f13cb533af3f07893 Mon Sep 17 00:00:00 2001
From: "J.A. de Jong @ vulgaris" <anne@amdj.nl>
Date: Sun, 13 Nov 2016 20:21:55 +0100
Subject: [PATCH] Put T0 and p in Gas class. Updated Globalconf

---
 src/CMakeLists.txt        |  2 +-
 src/material/air.h        |  2 +-
 src/material/gas.cpp      | 16 +++++++---
 src/material/gas.h        | 32 +++++++++++++-------
 src/material/helium.h     |  2 +-
 src/material/nitrogen.h   |  2 +-
 src/material/perfectgas.h |  2 +-
 src/sys/globalconf.cpp    | 63 +++++++++++++++++++++++++++++++++++++++
 src/sys/globalconf.h      | 54 +++++++++++++++++++++++++++++++++
 9 files changed, 156 insertions(+), 19 deletions(-)
 create mode 100644 src/sys/globalconf.cpp
 create mode 100644 src/sys/globalconf.h

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f8d7204..2966575 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -29,7 +29,7 @@ add_subdirectory(material)
 # add_subdirectory(mech)
 # add_subdirectory(seg)
 # add_subdirectory(sol)
-# add_subdirectory(sys)
+add_subdirectory(sys)
 # add_subdirectory(var)
 
 add_library(tasmet_src tasmet_tracer.cpp tasmet_exception.cpp tasmet_assert.cpp)
diff --git a/src/material/air.h b/src/material/air.h
index 1fa221d..8177768 100644
--- a/src/material/air.h
+++ b/src/material/air.h
@@ -15,7 +15,7 @@ class Air : public PerfectGas {
 protected:
     d Rs() const {return 287;}
 public:
-    Air():PerfectGas(air){}
+    Air(d T0,d p0):PerfectGas(air,T0,p0){}
     d cp(d T,d p) const;
     vd cp(const vd& T,const vd& p) const;
     d h(d T,d p) const;
diff --git a/src/material/gas.cpp b/src/material/gas.cpp
index f398d8a..943d860 100644
--- a/src/material/gas.cpp
+++ b/src/material/gas.cpp
@@ -12,16 +12,24 @@
 #include "air.h"
 #include "nitrogen.h"
 
-Gas* Gas::newGas(const GasType gastype) {
+#include "tasmet_constants.h"
+
+Gas* Gas::newGas(const GasType gastype,d T0,d p0) {
+
+    if(T0>constants::maxT || T0< constants::minT)
+        throw TaSMETError("Illegal reference temperature given");
+    if(p0>constants::maxp || p0< constants::minp)
+        throw TaSMETError("Illegal reference pressure given");
+    // End sanity checks
 
     switch (gastype) {
     case helium:
-        return new Helium();
+        return new Helium(T0,p0);
         break;
     case air:
-        return new Air();
+        return new Air(T0,p0);
     case nitrogen:
-        return new Nitrogen();
+        return new Nitrogen(T0,p0);
     default:
         FATAL("Gas type not known");
         break;
diff --git a/src/material/gas.h b/src/material/gas.h
index 3a401bc..3aba24b 100644
--- a/src/material/gas.h
+++ b/src/material/gas.h
@@ -12,7 +12,7 @@
 #include <armadillo>
 #include "tasmet_enum.h"
 #include "tasmet_types.h"
-
+#include "tasmet_constants.h"
 
 
 class Gas{
@@ -21,19 +21,31 @@ public:
 private:    
     GasType _gastype;
 protected:
-    Gas(GasType gastype):_gastype(gastype){}
+    d _T0,_p0;			/* Reference temperature and pressure */
+    
+    Gas(GasType gastype,d T0,d p0):_gastype(gastype),_T0(T0),_p0(p0){}
 public:
-
-
-
     Gas(const Gas& other) =delete;
     Gas& operator=(const Gas&)=delete;
     virtual ~Gas(){}
 
     // Static method to generate a Gas
-    static Gas* newGas(const GasType gastype);
+
+
+    static Gas* newGas(const GasType gastype,d T0=constants::T0,
+                     d p0=constants::p0);
+
     operator GasType() { return _gastype;}
 
+    d T0() const {return _T0;}
+    d p0() const {return _p0;}
+
+    // Speed of sound at the reference temperature and pressure
+    d c0() const {return cm(_T0,_p0);}
+    d rho0() const {return rho(_T0,_p0);}
+    d deltanu0(d freq) const{ return sqrt(2*mu(_T0,_p0)/(rho0()*2*number_pi*freq));}
+    d deltanu(d freq,d T,d p) const { return sqrt(2*mu(T,p)/(rho(T,p)*2*number_pi*freq)); }
+
     // Dimensionless numbers:
 
     // Specific heat ratio
@@ -49,6 +61,10 @@ public:
     virtual d rho(d T,d p) const=0;
     virtual vd rho(const vd& T,const vd& p) const=0;
 
+    // Adiabatic speed of sound
+    virtual d cm(d T,d p) const=0;
+    virtual vd cm(const vd& T,const vd& p) const=0;    
+
     // Internal energy [J/kg]
     virtual d e(d T,d p) const=0;
     virtual vd e(const vd& T,const vd& p) const=0;
@@ -68,10 +84,6 @@ public:
     virtual d beta(d T,d p) const=0;
     virtual vd beta(const vd& T,const vd& p) const=0;
 
-    // Adiabatic speed of sound [m/s]
-    virtual d cm(d T,d p) const=0;
-    virtual vd cm(const vd& T,const vd& p) const=0;
-
     // Dynamic viscosity [Pa*s]
     virtual d mu(d T,d p) const=0;
     virtual vd mu(const vd& T,const vd& p) const=0;
diff --git a/src/material/helium.h b/src/material/helium.h
index a422dc5..714e254 100644
--- a/src/material/helium.h
+++ b/src/material/helium.h
@@ -14,7 +14,7 @@ class Helium :public PerfectGas {
 protected:
     d Rs() const {return 2070;}
 public:
-    Helium():PerfectGas(helium){}
+    Helium(d T0,d p0):PerfectGas(helium,T0,p0){}
 
     d cp(d T,d p) const { return 5195;}
     vd cp(const vd& T,const vd& p) const;
diff --git a/src/material/nitrogen.h b/src/material/nitrogen.h
index c0fefb3..9e5068f 100644
--- a/src/material/nitrogen.h
+++ b/src/material/nitrogen.h
@@ -16,7 +16,7 @@ class Nitrogen : public PerfectGas {
 protected:
     d Rs() const {return 297;}
 public:
-    Nitrogen():PerfectGas(nitrogen){}
+    Nitrogen(d T0,d p0):PerfectGas(nitrogen,T0,p0){}
 
     d cp(d T,d p) const;
     vd cp(const vd& T,const vd& p) const;
diff --git a/src/material/perfectgas.h b/src/material/perfectgas.h
index 5ec3258..226be17 100644
--- a/src/material/perfectgas.h
+++ b/src/material/perfectgas.h
@@ -23,7 +23,7 @@
   
 class PerfectGas: public Gas {
 protected:
-    PerfectGas(GasType gastype): Gas(gastype){}
+    PerfectGas(GasType gastype,d T0,d p0): Gas(gastype,T0,p0){}
 
     // Not implemented for a perfect gas:
     // mu, kappa, h, cp, Rs
diff --git a/src/sys/globalconf.cpp b/src/sys/globalconf.cpp
new file mode 100644
index 0000000..9a5557d
--- /dev/null
+++ b/src/sys/globalconf.cpp
@@ -0,0 +1,63 @@
+#include "globalconf.h"
+#include "tasmet_constants.h"
+#include "tasmet_exception.h"
+#include "tasmet_io.h"
+
+Globalconf::Globalconf(us Nf,d freq)
+
+{
+
+    set(Nf,freq);
+    TRACE(10,"Globalconf constructor done");
+}
+void Globalconf::show() const {
+    cout << "------- Global configuration ------ \n";        
+    cout << "------- Number of harmonics to solve for: "<< _Nf <<"\n";
+    cout << "------- Fundamental frequency           : " << _omg/2/number_pi << " Hz\n";           
+}
+void Globalconf::set(us Nf,d freq){
+    TRACE(15,"Globalconf::set(_Nf,freq)");
+    d omg = 2*number_pi*freq;
+    //ctor
+    // Sanity checks
+    if(omg<constants::minomg && omg>constants::maxomg)
+        throw TaSMETError("Illegal frequency given");
+    if(Nf>=constants::maxNf)
+        throw TaSMETError("Too large number of frequencies given");
+
+    this->_Nf=Nf;
+    this->_omg=omg;
+
+    us Ns=this->Ns();
+    // Reinitialize all operators
+    iDFT_=zeros<dmat>(Ns,Ns);
+    fDFT_=zeros<dmat>(Ns,Ns);
+
+    DDTfd_=zeros<dmat>(Ns,Ns);
+    fDFT_.row(0).fill(1.0/double(Ns));
+
+    for(us i=1;i<=_Nf;i++){
+        for(us j=0; j<Ns;j++){
+            //Row i+1 (cosine components)
+            fDFT_(2*i-1,j)=2.0*cos(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
+            //Row i (sine components)
+            fDFT_(2*i,j)=-2.0*sin(2.0*number_pi*double(i)*double(j)/double(Ns))/double(Ns);
+        }
+    }
+
+    iDFT_.col(0).fill(1.0);	// Steady part
+    for(us k=0;k<Ns;k++){
+        for (us n=1;n<=_Nf;n++){
+            iDFT_(k,2*n-1)=cos(2.0*number_pi*double(n)*double(k)/Ns);
+            iDFT_(k,2*n)=-sin(2.0*number_pi*double(n)*double(k)/Ns);
+        }
+    }
+
+    for(us i=1;i<=_Nf;i++){
+        DDTfd_(2*i-1,2*i  )=-double(i)*_omg;
+        DDTfd_(2*i  ,2*i-1)=double(i)*_omg;
+    }
+
+}
+
+
diff --git a/src/sys/globalconf.h b/src/sys/globalconf.h
new file mode 100644
index 0000000..58a5782
--- /dev/null
+++ b/src/sys/globalconf.h
@@ -0,0 +1,54 @@
+// globalconf.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Global configuration options
+//////////////////////////////////////////////////////////////////////
+#pragma once
+
+#ifndef _GLOBALCONF_H_
+#define _GLOBALCONF_H_
+
+#include "tasmet_types.h"
+#include "tasmet_tracer.h"
+
+class Globalconf{
+    d _omg;		// The "base" frequency in rad/s
+    us _Nf;			// Number of frequencies to solve for
+
+    dmat fDFT_,iDFT_,DDTfd_;
+    // d Wfo_=0;			// First order 'upwind' factor. If
+    // Wfo=-1, interpolation is done from
+    // the left side. If Wfo=0,
+    // interpolation is second order. If
+    // Wfo=1, interpolation is done from
+    // the right side
+public:
+    Globalconf(us Nf,d freq);
+    
+    us Nf() const {return _Nf;}
+    us Ns() const {return 2*_Nf+1;}    
+
+    ~Globalconf(){TRACE(-5,"~Globalconf()");}
+    d getomg() const {return _omg;}
+    d getfreq() const {return _omg/2/number_pi;}
+    // d meshPeclet(const Gas& gas,d dx,d u) const {return u*dx*gas.rho0()*gas().cp(T0())/gas().kappa(T0());}
+
+    void set(us Nf,d freq);	// Set data for new frequency and
+    // number of samples
+
+    void setNf(us Nf){set(Nf,getfreq());}  
+    void setfreq(d freq){set(Nf(),freq);}
+    void setomg(d omg)  {set(Nf(),omg/(2*number_pi));}
+
+    const dmat& iDFT() const {return iDFT_;} //inverse discrete Fourier transform matrix
+    const dmat& fDFT() const {return fDFT_;} //forward discrete Fourier transform matrix
+    const dmat& DDTfd() const {return DDTfd_;}//Derivative in frequency domain
+
+    void show() const;
+
+}; /* Class Globalconf */
+
+#endif /* _GLOBALCONF_H_ */
+//////////////////////////////////////////////////////////////////////