Added postprocessor to this repo. Bugfixes and a lot of stuff

doc/usg.lyx

@@ -290,7 +290,7 @@ TaSMET
 Open source and free to use.
  As the author is aware of the possibilities created by using open source
  codes, the choice of publishing this computer code as open source was not
- a hard choice.
+ a hard one.
 \end_layout
 
 \begin_layout Itemize
@@ -370,6 +370,166 @@ TaSMET
 \noun default
 , the periodic steady-state of a TA system is simulated using a numerical
  model.
+ The convergence speed is much accelerated by assuming a periodic steady
+ state of the system.
+ When a system is in periodic steady state, every physical quantity 
+\begin_inset Formula $\xi$
+\end_inset
+
+ can be described by a Fourier series:
+\begin_inset Formula 
+\begin{equation}
+\xi(t)=\sum_{n=0}^{\infty}\Re\left[\hat{\xi}_{n}e^{in\omega t}\right],
+\end{equation}
+
+\end_inset
+
+where 
+\begin_inset Formula $\hat{\xi}_{n}$
+\end_inset
+
+ denote the Fourier coefficients of the quantity 
+\begin_inset Formula $\xi$
+\end_inset
+
+, which can be computed as
+\begin_inset Note Note
+status collapsed
+
+\begin_layout Plain Layout
+\begin_inset Formula $\int\limits _{0}^{T}\xi(t)e^{-im\omega t}\mathrm{d}t=\int\limits _{0}^{T}\sum\limits _{n=0}^{\infty}\Re\left[\hat{\xi}_{n}e^{in\omega t}\right]e^{-im\omega t}\mathrm{d}t$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $\int\limits _{0}^{T}\xi(t)e^{-im\omega t}\mathrm{d}t=\int\limits _{0}^{T}\sum\limits _{n=0}^{\infty}\frac{1}{2}\left(\hat{\xi}_{n}e^{in\omega t}+\hat{\xi}_{n}^{*}e^{-in\omega t}\right)e^{-im\omega t}\mathrm{d}t$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+For 
+\begin_inset Formula $m=0$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $\int\limits _{0}^{T}\xi(t)\mathrm{d}t=T\hat{\xi}_{0}$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+And for 
+\begin_inset Formula $m\neq0$
+\end_inset
+
+:
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $\hat{\xi}_{n}=\frac{2}{T}\int\limits _{0}^{T}\xi(t)e^{-im\omega t}\mathrm{d}t$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Formula 
+\begin{equation}
+\hat{\xi}_{n}=\frac{2}{T}\int\limits _{0}^{T}\xi(t)e^{-im\omega t}\mathrm{d}t
+\end{equation}
+
+\end_inset
+
+ Now, to solve the Fourier coefficients, only a finite number of terms (
+\begin_inset Formula $N$
+\end_inset
+
+) can be taken into account.
+ This also results in a finite time resolution.
+ Writing both the time and frequency domain in a discrete form:
+\begin_inset Note Note
+status collapsed
+
+\begin_layout Plain Layout
+\begin_inset Formula $t=m\Delta t$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $T=M\Delta t$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $M=\Delta t/T=\Delta t\omega$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\begin_inset Formula 
+\begin{equation}
+\xi_{m}\equiv\xi(t_{m})=\sum_{n=0}^{N-1}\Re\left[\hat{\xi}_{n}e^{inm/M}\right].
+\end{equation}
+
+\end_inset
+
+This can be written in matrix-vector form as
+\begin_inset Formula 
+\begin{equation}
+\xi_{m}=\boldsymbol{f}_{m}^{-1}\cdot\hat{\boldsymbol{\xi}},
+\end{equation}
+
+\end_inset
+
+or, for all time-instances 
+\begin_inset Formula 
+\begin{equation}
+\boldsymbol{\xi}=\boldsymbol{\mathcal{F}}^{-1}\cdot\sum_{n=0}^{N-1}\Re\left[in\omega\hat{\xi}_{n}e^{in\omega t}\right]
+\end{equation}
+
+\end_inset
+
+The time-derivative can be easily obtained
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial\boldsymbol{\xi}}{\partial t}=\boldsymbol{\mathcal{F}}^{-1}\cdot\boldsymbol{\text{\omega}}\cdot\boldsymbol{\mathcal{F}}\cdot\boldsymbol{\xi},
+\end{equation}
+
+\end_inset
+
+where
+\begin_inset Formula 
+\begin{equation}
+\boldsymbol{\omega}=\left[\begin{array}{ccccc}
+0 & 0 & 0 & 0 & \dots\\
+0 & i\omega & 0 & 0 & \dots\\
+0 & 0 & 2i\omega & 0 & \dots\\
+0 & 0 & 0 & 3i\omega & \dots\\
+\vdots & \vdots & \vdots & \vdots & \ddots
+\end{array}\right]
+\end{equation}
+
+\end_inset
+
+
 \end_layout
 
 \begin_layout Chapter

src/duct/duct.cpp

@@ -111,14 +111,18 @@ void Duct::residualJac(SegPositionMapper& residual,
         pp   = getvart(constants::p,gp+1);
 
         cont  = ((rhop%up)-(rho%u))/dx;
-        // cont += sys.DDTtd()*rho;
+        cont += sys.DDTtd*rho;
         
         *residual.at(i) = cont;
 
-        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] =
-            -diagmat(rho)/dx;
         jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
             -diagmat(u)/dx;
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] +=
+            sys.DDTtd;
+        
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] =
+            -diagmat(rho)/dx;
+        
         jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] =
             diagmat(up)/dx;
         jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] =
@@ -127,19 +131,33 @@ void Duct::residualJac(SegPositionMapper& residual,
         i++;
 
         mom = (rhop%up%up - rho%u%u + pp - p)/dx;
+        mom += sys.DDTtd *(rho%u);
+
         jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] =
             -diagmat(rho%u)/dx;
+
+        // Time-derivative term
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::u}] +=
+            sys.DDTtd*diagmat(rho);
+
         jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] =
             -diagmat(u%u)/dx;
+
+        // Time-derivative term
+        jacrows.at(i)[{_id,gp*nvars_per_gp+constants::rho}] +=
+            sys.DDTtd*diagmat(u);
+
         jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::rho}] =
-            diagmat(u%u)/dx;
+            diagmat(up%up)/dx;
         jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::u}] =
             diagmat(up%rhop)/dx;
         jacrows.at(i)[{_id,(gp)*nvars_per_gp+constants::p}] =
             -eye(Ns,Ns)/dx;
         jacrows.at(i)[{_id,(gp+1)*nvars_per_gp+constants::p}] =
             eye(Ns,Ns)/dx;
-        *residual.at(i) = mom;
+ 
+
+       *residual.at(i) = mom;
 
         i++;
         
@@ -181,7 +199,7 @@ void Duct::residualJac(SegPositionMapper& residual,
         *residual.at(i) = st; 
 
         i++;
-    }
+    } // end of for loop over gridpoints
 
     // Equation of state for the last node
     st = gas.rho(Tp,pp) - rhop;
@@ -190,9 +208,9 @@ void Duct::residualJac(SegPositionMapper& residual,
     jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] =
         -eye(Ns,Ns);
     jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] =
-        diagmat(gas.drhodp(T,p));
+        diagmat(gas.drhodp(Tp,pp));
     jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::T}] =
-        diagmat(gas.drhodT(T,p));
+        diagmat(gas.drhodT(Tp,pp));
 
     i++;
 
@@ -205,12 +223,12 @@ void Duct::residualJac(SegPositionMapper& residual,
         d rho0 = gas.rho0();
         d gamma0 = gas.gamma(T0,p0);
 
-        en = p/p0 - pow(rho/rho0,gamma0);
+        en = pp/p0 - pow(rhop/rho0,gamma0);
 
         jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::p}] =
             eye(Ns,Ns)/p0;
         jacrows.at(i)[{_id,(ngp()-1)*nvars_per_gp+constants::rho}] =
-            -(gamma0/rho0)*diagmat(pow(rho/rho0,gamma0-1));
+            -(gamma0/rho0)*diagmat(pow(rhop/rho0,gamma0-1));
 
         *residual.at(i) = en; i++;
     }
@@ -231,6 +249,11 @@ PosId Duct::getDof(int varnr,int gp) const {
 
 vd Duct::getvart(int varnr,int gp) const {
     TRACE(14,"Duct::getvart()");
+
+    tasmet_assert(varnr >=0 &&
+                  varnr < constants::varnr_max,
+                  "Illegal variable number");
+
     const SegPositionMapper& sol = sys.getSolution(_id);
 
     us Ns = sys.Ns();
@@ -268,11 +291,11 @@ void Duct::initialSolution(SegPositionMapper& sol) const {
     for(us i=0;i<ngp();i++){
         VARTRACE(15,i);
         // Initial density
-        *sol.at(segdof) += gas.rho0()+0.1;
+        *sol.at(segdof) += gas.rho0();
         segdof++;
         
         // Initial velocity
-        *sol.at(segdof) += 0.1;
+        // *sol.at(segdof) += 0.1;
         segdof++;
 
         // Initial Temperature
@@ -396,6 +419,27 @@ void Duct::exportHDF5(const hid_t group_id) const {
     }
     p.exportHDF5(group_id);
 
+    TXData rho;
+    rho.name = "Fluid density";
+    rho.unit = "kg/m$^3$";
+    rho.symbol = "rho";
+    rho.x = dmat(sys.Ns(),ngp());
+    for(us gp=0;gp<ngp();gp++){
+        rho.x.col(gp) = rhot(gp);
+    }
+    rho.exportHDF5(group_id);
+
+    TXData m;
+    m.name = "Mass flow";
+    m.unit = "kg/s";
+    m.symbol = "m";
+    m.x = dmat(sys.Ns(),ngp());
+    for(us gp=0;gp<ngp();gp++){
+        m.x.col(gp) = this->S(gp) * rho.x.col(gp) % u.x.col(gp);
+    }
+    m.exportHDF5(group_id);
+    
+
     // TXData Ts;
     // Ts.name = "Solid temperature";
     // Ts.unit = "K";

src/duct/geom.cpp

@@ -12,6 +12,7 @@
 #include <memory>
 #include "tasmet_exception.h"
 #include "tasmet_tracer.h"
+#include "tasmet_assert.h"
 
 Geom::Geom(const pb::Duct& duct) {
 
@@ -98,16 +99,16 @@ d Geom::getFluidVolume() const {
   }
   d Geom::getSolidVolume() const {
     d Vs=0;
-    int i=0;
 
-    Vs+=0.5*(x(i+1)-x(i))*Ss(i); // First 'cell'
+    tasmet_assert(false,"Not updated function");
+    // Vs+=0.5*(x(i+1)-x(i))*Ss(i); // First 'cell'
 
-    for(int i=1;i<ngp()-1;i++){    // Middle cells
+    for(us i=1;i<ngp()-1;i++){    // Middle cells
-        Vs+=0.5*(x(i)-x(i-1))*Ss(i);
+        // Vs+=0.5*(x(i)-x(i-1))*Ss(i);
-        Vs+=0.5*(x(i+1)-x(i))*Ss(i);
+        // Vs+=0.5*(x(i+1)-x(i))*Ss(i);
     }
 
-    Vs+=0.5*(x(i)-x(i-1))*Ss(i); // Last 'cell'
+    // Vs+=0.5*(x(i)-x(i-1))*Ss(i); // Last 'cell'
     return Vs;
   }
 

src/ductbc/adiabaticwall.cpp

@@ -16,7 +16,6 @@ AdiabaticWall::AdiabaticWall(const TaSystem& sys,
                              const us id,
                              const pb::DuctBc& dbc):
     DuctBc(sys,id,dbc)
-    
 {
     TRACE(15,"AdiabaticWall(id,sys,dbc)");
     tasmet_assert(dbc.type() == pb::AdiabaticWall,"Wrong type given to constructor");
@@ -47,7 +46,7 @@ void AdiabaticWall::residualJac(SegPositionMapper& residual,
 
     us i=0;
 
-    if(_side == pb::left) {
+    if(_dbc.side() == pb::left) {
         *residual.at(i) = duct.ut(0);
         jac.at(i)[duct.getDof(constants::u,0)].eye(Ns,Ns);
 
@@ -57,7 +56,8 @@ void AdiabaticWall::residualJac(SegPositionMapper& residual,
             tasmet_assert(false,"");
         }
     }
-    else {
+    else if(_dbc.side() == pb::right) {
+        TRACE(25,"AdWall RIGHT");
         *residual.at(i) = duct.ut(-1);
         jac.at(i)[duct.getDof(constants::u,-1)].eye(Ns,Ns);
         if(dpb.htmodel() != pb::Isentropic) {
@@ -66,6 +66,9 @@ void AdiabaticWall::residualJac(SegPositionMapper& residual,
             tasmet_assert(false,"");
         }
     }
+    else {
+        throw TaSMETError("Illegal side parameter given in AdiabaticWall.");
+    }
 
 }
 

src/ductbc/ductbc.h

@@ -15,6 +15,7 @@ class TaSystem;
 class Duct;
 
 class DuctBc :public Segment {
+protected:
     pb::DuctBc _dbc;
 public:
     DuctBc(const TaSystem& sys,

src/ductbc/pressurebc.cpp

@@ -41,6 +41,7 @@ PressureBc::PressureBc(const TaSystem& sys,
     }
 
     else {
+        TRACE(25,"PressureBc not Isentropic");
         EvaluateFun Tfun(dbc.temperature(),
                          "Error in evaluating prescribed temperature",
                          "t");
@@ -58,10 +59,10 @@ PressureBc::PressureBc(const TaSystem& sys,
                        const PressureBc& o):
     DuctBc(sys,o),
     _p(o._p),
-    _T(o._T) {
+    _T(o._T),
-
+    _Ts(o._Ts)
+{
     TRACE(15,"PressureBc(o)");
-
 }
 PressureBc::~PressureBc() {
     
@@ -78,21 +79,30 @@ void PressureBc::residualJac(SegPositionMapper& residual,
     const pb::Duct& dpb = getDuct().getDuctPb();
     us Ns = sys.Ns();
 
+    // VARTRACE(25,_T-sys.gas().T0);
+    
     const Duct& duct = getDuct();
-    if(_side == pb::left) {
+    if(_dbc.side() == pb::left) {
+        TRACE(25,"PressureBc LEFT");
         *residual.at(0) = duct.pt(0) - _p;
-        jacrows.at(0)[duct.getDof(constants::p,0)] = eye(Ns,Ns);
+        jacrows.at(0)[duct.getDof(constants::p,0)].eye(Ns,Ns);
         if(dpb.htmodel() != pb::Isentropic) {
+            tasmet_assert(false,"");
             // residual.subvec(Ns,2*Ns-1) = duct.Tt(0) - _T;
         }
     }
-    else {
+    else if(_dbc.side() == pb::right) {
+        TRACE(25,"PB RIGHT");
         *residual.at(0) = duct.pt(-1) - _p;
-        jacrows.at(0)[duct.getDof(constants::p,-1)] = eye(Ns,Ns);
+        jacrows.at(0)[duct.getDof(constants::p,-1)].eye(Ns,Ns);
         if(dpb.htmodel() != pb::Isentropic) {
+            tasmet_assert(false,"");
             // residual.subvec(Ns,2*Ns-1) = duct.Tt(-1) - _T;
         }
     }
+    else {
+        throw TaSMETError("Illegal side parameter given in PressureBc.");
+    }
 }
 
 us PressureBc::getNEqs() const {

src/ductbc/pressurebc.h

@@ -25,7 +25,6 @@ class Variable;
 class PressureBc: public DuctBc {
     vd _p,_T,_Ts;               /**< Prescribed values for pressure,
                                    temperature and solid temperature */
-    pb::DuctSide _side;         /**< Duct side at which this b.c. works */
 protected:
     PressureBc(const TaSystem&,const PressureBc&);
 public:

src/gui/add_ductbc_dialog.cpp

@@ -49,6 +49,7 @@ AddDuctBcDialog::AddDuctBcDialog(const std::string& name,QWidget* parent):
 }
 void AddDuctBcDialog::changed() {
     if(_init) return;
+    TRACE(16,"AddDuctBcDialog::changed() continue");
     pb::DuctBcType type = (pb::DuctBcType) _dialog->type->currentIndex();
     bool isentropic =  _dialog->isentropic->isChecked();
     VARTRACE(15,type);

src/gui/add_ductbc_dialog.h

@@ -30,7 +30,11 @@ private:
 
 private slots:
     void on_type_currentIndexChanged(int) {changed();}
+    void on_duct_id_valueChanged(int) { changed();}
+    void on_side_currentIndexChanged(int) {changed();}    
     void on_isentropic_stateChanged(int) {changed();}
+    void on_pressure_textChanged() {changed();}
+    void on_temperature_textChanged() {changed();}
     void accept();
     void reject();
 

src/gui/mainwindow.cpp

@@ -204,6 +204,12 @@ void TaSMETMainWindow::saveAsModel() {
                                                     filetype);
 
     if(fileName.size()!=0) {
+
+        // Add extension if necessary
+        if(!fileName.contains(constants::model_fileext)) {
+            fileName += constants::model_fileext;
+        }
+
         string fn = fileName.toStdString();
 
         saveModel(&fn);
@@ -327,13 +333,14 @@ int add_edit_segment(QWidget* parent,
     dialogtype dialog(name,parent);
     
     if(segmap.find(id) != segmap.end()) {
+        TRACE(16,"Found the segment in the list");
         dialog.set(segmap[id]);
     }
 
     int exitcode = dialog.exec();
     
     if(exitcode == QDialog::Accepted) {
-
+        TRACE(16,"Updating existing segment");
         segmap[id] = dialog.get();
         
     }
@@ -424,7 +431,7 @@ void TaSMETMainWindow::on_actionSolve_triggered() {
         
         if(_model.solution_size()>0) {
             vd solution(_model.solution_size());
-            for(us i=0;i<_model.solution_size();i++) {
+            for(int i=0;i<_model.solution_size();i++) {
                 solution(i) = _model.solution(i);
             }
 

src/gui/solver_dialog.cpp

@@ -25,8 +25,8 @@ SolverDialog::SolverDialog(QWidget* parent,
                            const GradientNonlinearSystem& sys,
                            pb::SolverParams& sparams):
     QDialog(parent),
-    _sys(sys),
     _sparams(sparams),
+    _sys(sys),
     _dialog(new Ui::solver_dialog())
 {
     

src/gui/solver_dialog.h

@@ -29,18 +29,22 @@ class SolverProgress;
 class SolverDialog: public QDialog {
     Q_OBJECT
 
-    pb::SolverParams& _sparams;
+    pb::SolverParams& _sparams; /**< Reference to solver params as
+                                   given in constructor. */
+
+    const GradientNonlinearSystem& _sys;
 
     Ui::solver_dialog* _dialog;
 
-    bool _init = true;
+    bool _init = true;          /**< Disables Qt callbacks */
 
-    QCustomPlot* _plot;
+    QCustomPlot* _plot = nullptr;
-    QCPGraph *_funer,*_reler,*_funtol,*_reltol;
+    QCPGraph *_funer=nullptr,*_reler=nullptr,
+        *_funtol=nullptr,*_reltol=nullptr;
 
     SolverWorker* _solver_worker = nullptr;
 
-    const GradientNonlinearSystem& _sys;
+
 
     /// Place where the final solution will be stored
     vd _solution;

src/material/perfectgas.h

@@ -67,7 +67,7 @@ public:
         return -p/Rs()/pow(T,2);
     }
     d drhodp(d T,d p) const {
-        return Rs()/T;
+        return 1/(Rs()*T);
     }
 
     d cv(d T,d p) const {

src/solver/solver.h

@@ -97,7 +97,8 @@ public:
 
 protected:
     /// This member fcn should be implemented by the Solver instance.
-    virtual void start_implementation(system_T& sys,progress_callback*)=0;
+    virtual void start_implementation(system_T& sys,
+                                      progress_callback*)=0;
 };
 
 

src/sys/globalconf.cpp

@@ -14,9 +14,9 @@ GlobalConf::GlobalConf(us Nf,d freq):
     _Nf(Nf)
 {
 
-    TRACE(10,"GlobalConf constructor done");
+    TRACE(10,"GlobalConf()");
 
-    if(Nf>=constants::max_Nf)
+    if(Nf>constants::max_Nf)
         throw TaSMETError("Too large number of frequencies given");
 
 
@@ -38,7 +38,7 @@ GlobalConf::GlobalConf(us Nf,d freq):
     }
 
     _iDFT.col(0).fill(1.0);	// Steady part
-    for(us k=0;k<Ns;k++){
+    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);
@@ -71,10 +71,10 @@ void GlobalConf::setomg(d omg){
     this->_omg=omg;
 
     for(us i=1;i<=_Nf;i++){
-        _DDTfd(2*i-1,2*i  )=-double(i)*_omg;
+        _DDTfd(2*i-1,2*i  )=- ((d) i)*_omg;
-        _DDTfd(2*i  ,2*i-1)=double(i)*_omg;
+        _DDTfd(2*i  ,2*i-1)= ((d) (i))*_omg;
     }
-    _DDTtd = _iDFT * _DDTfd * _fDFT;
+    _DDTtd = _iDFT * (_DDTfd * _fDFT);
 }
 
 //////////////////////////////////////////////////////////////////////

src/sys/tasystem.cpp

@@ -5,7 +5,7 @@
 // Description:
 //
 //////////////////////////////////////////////////////////////////////
-
+#define TRACERPLUS (-10)
 #include "segment.h"
 #include "tasystem.h"
 #include "triplets.h"

src/tasmet_constants.h

@@ -101,7 +101,8 @@ namespace constants {
     const int T=2;
     const int p=3;
     const int Ts=4;
-  
+    const int varnr_max = 5;
+    
     const char* const model_fileext = ".tasmet";
 
 } // namespace constants

src/tasmet_solvemodel.cpp

@@ -16,17 +16,35 @@
 
 void usage(const char* fn) {
     cout << "Usage: " << endl;
-    cout << fn << " [Input file] " << endl;
+    cout << fn << " [Input file] ";
+    #ifdef TASMET_DEBUG
+    cout << "[Tracer_Level]";
+    #endif
+    cout << endl;
 }
 
 int main(int argc, char *argv[]) {
 
-    INITTRACE(10);
 
-    if(argc != 2) {
+
+    if(argc < 2 || argc > 3) {
         usage(argv[0]);
         return -1;
     }
+    #ifdef TASMET_DEBUG
+    int trace_level = 20;
+
+    if(argc >= 3) {
+        try {
+            trace_level = std::stoi(argv[2]);
+        }
+        catch(...) {
+            usage(argv[0]);
+            return -1;
+        }
+    }
+    INITTRACE(trace_level);
+    #endif
 
     pb::Model model;
     

tasmet_postprocessor.py

@@ -0,0 +1,576 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Feb  2 08:33:53 2017
+
+@author: anne
+"""
+
+# Here we provide the necessary imports.
+# The basic GUI widgets are located in QtGui module. 
+import sys, h5py
+import numpy as np
+from PyQt5.QtCore import *
+from PyQt5.QtWidgets import *
+from PyQt5.QtGui import QIcon,QDoubleValidator
+from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
+from matplotlib.backends.backend_qt5 import NavigationToolbar2QT as NavigationToolbar
+from matplotlib.figure import Figure
+from gui.about_dialog import Ui_about_dialog
+
+scale = [('nano',-9,'n'),
+         ('micro',-6,'$\mu$'),
+         ('milli:',-3,'m'),
+         ('centi',-2,'c'),
+         ('deci',-1,'d'),
+         ('-',0,''),
+         ('deca',1,'da'),
+         ('hecto',2,'h'),
+         ('kilo',3,'k'),
+         ('mega',6,'M'),
+         ('giga',9,'G'),
+         ('tera',12,'T')]
+
+appname = 'TaSMET Postprocessor'
+def getGuiSettings():
+    return QSettings(appname,appname)
+
+ 
+class MaxCounter(object):
+    def __init__(self,max=-1):
+        self.a = 0
+        self.max = max            
+    def __call__(self): 
+        rval = self.a
+        self.a += 1
+        if self.a == self.max+1:
+            self.a = 0
+        return rval
+
+class StepCounter(object):
+    def __init__(self,step=1):
+        self.a = 0
+        self.step = step
+        self.inner = 0
+        
+    def __call__(self): 
+        rval = self.a
+        
+        self.inner += 1
+        
+        if self.inner == self.step:
+            self.inner = 0
+            self.a += 1
+     
+        return rval    
+
+class TaSMETFigure(Figure):
+    
+    def __init__(self, width=5, height=4, dpi=100):
+        
+        Figure.__init__(self,figsize=(width, height), dpi=dpi)
+        
+        self._axes = self.add_subplot(111)
+
+        # We want the axes cleared every time plot() is called
+        
+        self._plots = []
+        self._grid = False
+    def addPlot(self,plot,one):
+        if one:
+            self._plots = [plot]
+        else:
+            self._plots.append(plot)
+        
+        self.rePlot()
+        
+    def rePlot(self):
+        self._axes.clear()
+        for plot in self._plots:
+            self._axes.plot(plot['x'],plot['y'],label=plot['ylabel'])
+        
+            self._axes.set_ylabel(plot['ylabel'])
+            self._axes.set_xlabel(plot['xlabel'])
+            
+        self._axes.grid(self._grid,'both')
+        
+        self.canvas.draw_idle()
+    def setGrid(self,grid):
+        self._grid = grid
+        self.rePlot()
+        
+    def clearLastPlot(self):
+        if len(self._plots) > 0:
+            del self._plots[-1]
+        self.rePlot()
+    def clearAllPlot(self):
+        self._plots = []
+        self.rePlot()
+
+        
+class TaSMETCanvas(FigureCanvas):
+
+    def __init__(self, parent, figure):
+        
+        FigureCanvas.__init__(self, figure)
+        self.setParent(parent)
+
+        FigureCanvas.setSizePolicy(self,
+                                   QSizePolicy.Expanding,
+                                   QSizePolicy.Expanding)
+        
+        FigureCanvas.updateGeometry(self)
+
+
+class MainWindow(QMainWindow):
+    
+    def __init__(self, file_):
+        QMainWindow.__init__(self)
+
+        # Open the postprocessing file
+        try:
+            self._file = h5py.File(file_,'r')
+        except:
+            QMessageBox.warning(self,
+            'Error opening file',
+            'File %s is not a valid TaSMET Postprocessing file. Closing.' %file_)
+            exit(1)
+
+        # Suppress triggering of the GUI changed() callbacks
+        self._suppressed = True
+
+        self.setWindowTitle(appname)
+        
+        # Window icon
+        self.setWindowIcon(QIcon('images/tasmet_small.png'))        
+        
+        self._figure = TaSMETFigure()
+        self._canvas= TaSMETCanvas(self,figure=self._figure)        
+
+        self._seg_box = QComboBox()
+        self._seg_box.currentIndexChanged.connect(self.segItemChanged)
+        
+        self._qty_box = QComboBox()
+        self._qty_box.currentIndexChanged.connect(self.qtyItemChanged)
+        
+        
+        # Horizontal offset
+        self._hor_offset = QLineEdit()
+        self._hor_offset.setAlignment(Qt.AlignRight)
+        self._hor_offset.setText('0.0')
+        self._hor_offset.setMinimumWidth(40)
+
+        # Vertical offset
+        self._ver_offset = QLineEdit()
+        self._ver_offset.setAlignment(Qt.AlignRight)
+        self._ver_offset.setText('0.0')
+        self._ver_offset.setMinimumWidth(40)
+                             
+        offset_validator = QDoubleValidator()
+        self._hor_offset.setValidator(offset_validator)
+        self._ver_offset.setValidator(offset_validator)
+        
+        
+        self._scale = QComboBox()
+        for val in scale:
+            self._scale.addItem(val[0] + ' (%0.0e)' %(10**val[1]))
+        self._scale.setCurrentIndex(5)
+        
+        self._hold_button = QPushButton('Hold on')
+        self._hold_button.setCheckable(True)
+        
+        self._vstime = QRadioButton('Plot vs time')
+        self._vstime.clicked.connect(self.qtyItemChanged)
+        
+        self._vspos = QRadioButton('Plot vs position')
+        self._vspos.clicked.connect(self.qtyItemChanged)
+        
+        self._vspos.setChecked(True)
+        self._vs_instance = QSpinBox()
+        
+        self._plot_button = QPushButton('Plot')
+        self._plot_button.clicked.connect(self.plot)
+
+        self._output_button = QPushButton('Output')
+        self._output_button.clicked.connect(self.output)
+        
+        self._animate_button = QPushButton('Animate')
+        self._animate_button.clicked.connect(self.animate)
+                
+        
+        top = QGridLayout()
+
+
+                
+        col = StepCounter(step=2)
+        row = MaxCounter(max=1)
+        
+        top.addWidget(QLabel('Segment:'),row(),col())
+        top.addWidget(QLabel('Quantity:'),row(),col())
+        top.addWidget(self._seg_box,row(),col())
+        top.addWidget(self._qty_box,row(),col())
+        
+        top.addWidget(QLabel('Horizontal offset'),row(),col())
+        top.addWidget(QLabel('Vertical offset'),row(),col())
+
+        top.addWidget(self._hor_offset,row(),col())
+        top.addWidget(self._ver_offset,row(),col())
+              
+        top.addWidget(QLabel('Scale:'),row(),col())
+        row(); col()
+        top.addWidget(self._scale,row(),col())
+        row(); col()
+
+        radiobox = QHBoxLayout()
+        radiobox.addWidget(self._vstime)
+        radiobox.addWidget(self._vspos)
+
+        top.addLayout(radiobox,row(),col())
+        
+        instance_box = QHBoxLayout()
+        instance_box.addWidget(QLabel('Time/position instance:'))
+        instance_box.addWidget(self._vs_instance)
+        
+        top.addLayout(instance_box,row(),col())
+        
+        top.addWidget(self._plot_button,row(),col())
+        top.addWidget(self._hold_button,row(),col())
+        
+        top.addWidget(self._output_button,row(),col())
+        top.addWidget(self._animate_button,row(),col())
+        
+        htop = QHBoxLayout()
+        htop.addLayout(top)
+        htop.addStretch()
+
+        vlayout = QVBoxLayout()
+        vlayout.addLayout(htop)
+        navbar = NavigationToolbar(self._canvas, self)
+        
+        vlayout.addWidget(navbar)
+        vlayout.addWidget(self._canvas)
+        
+        hlayout = QHBoxLayout()
+        
+        output_layout = QVBoxLayout()        
+        output_layout.addWidget(QLabel('Output'))
+        
+        self._output = QPlainTextEdit()
+        self._output.setFixedWidth(300)
+        
+        self._output.setStyleSheet('font-family: monospace; background-color: white; color: black;')
+        
+        output_layout.addWidget(self._output)
+        
+        hlayout.addLayout(vlayout)        
+        hlayout.addLayout(output_layout)
+
+        # Put stuff in a widget and put it in the main window
+        widget = QWidget()        
+        widget.setLayout(hlayout)
+        self.setCentralWidget(widget)
+    
+        # File menu
+        file_menu = QMenu('&File', self)
+        file_menu.addAction('&Open',self.loadFile)
+        file_menu.addAction('&Exit',self.close)        
+
+        # Plot menu
+        plot_menu = QMenu('&Plot', self)
+        self._grid_option = QAction('Enable grid')
+        self._grid_option.setCheckable(True)
+        self._grid_option.triggered.connect(self._figure.setGrid)
+        plot_menu.addAction(self._grid_option)
+        plot_menu.addAction('&Delete last line',self._figure.clearLastPlot)
+        plot_menu.addAction('&Delete all lines',self._figure.clearAllPlot)
+
+        # Options menu
+        options_menu = QMenu('&Options', self)
+        options_menu.addAction('Clear output',self.clearOutput)
+        options_menu.addAction('Default output',self.defaultOutput)
+        
+        # Help menu
+        help_menu = QMenu('&Help', self)
+        help_menu.addAction('&About', self.about)
+        
+        # Add menus to the menubar        
+        self.menuBar().addMenu(file_menu)
+        self.menuBar().addMenu(plot_menu)
+        self.menuBar().addMenu(options_menu)
+        self.menuBar().addMenu(help_menu)
+
+        # Set window sizes
+        settings = getGuiSettings()
+        if settings.value('Geometry') is not None:
+            self.restoreGeometry(settings.value("Geometry"))
+        
+        self.initGui()
+        self._suppressed = False
+    
+    def clearOutput(self):
+        self._output.setPlainText('')
+    
+    def defaultOutput(self):
+        txt = ''
+        nsegments = 0
+        for x in self._file:
+            if isinstance(self._file[x], h5py.Group):
+                nsegments+=1
+        basestr = '{:22}: {:>12} [{:2}]\n'
+        txt += basestr.format('Number of segments',nsegments,'-')
+        txt += basestr.format('Fundamental frequency',self._file.attrs['freq'][0],'Hz')
+        txt += basestr.format('Numer of time samples',self._file.attrs['Ns'][0],'-')
+        txt += basestr.format('Reference temperature',self._file.attrs['T0'][0],'K')
+        txt += basestr.format('Reference pressure',self._file.attrs['p0'][0],'Pa')
+        txt += basestr.format('Gas type',u''+self._file.attrs['gas'].decode('utf-8'),'-')
+        tinst = self._file.attrs['t']
+        txt += '\nTime instances:\n'
+        for i in range(tinst.shape[0]):
+            txt += '{:2}: {:2e}\n'.format(i,tinst[i])
+            
+        
+        
+        self._output.setPlainText(txt)
+    
+    def output(self):
+        print('output')
+    
+    def animate(self):
+        print('animate')
+    
+    def getCurrentData(self):
+        """
+        Returns the current datatype, None if invalid or not selected        
+        """
+        segitemno = self._seg_box.currentText()
+        itemtxt = self._qty_box.currentText()
+        if len(segitemno) == 0:
+            return None
+        if len(itemtxt)==0:
+            return None
+        else:
+            data = self._file[segitemno][itemtxt]
+            datatype = data.attrs['datatype'].decode('utf-8')
+            return data,datatype
+        
+  
+        
+    
+    def segItemChanged(self,item):
+        """
+        Update the combobox with quantities
+        """
+        self._suppressed = True
+        self._qty_box.clear()
+        one = False
+        for dataset in self._file[str(item)]:
+            self._qty_box.addItem(dataset)
+            one = True
+        
+        # Not at least one plottable item in the list of the selected segment
+        self._plot_button.setEnabled(one)
+        self._output_button.setEnabled(one)
+            
+        self._suppressed = False
+        self._qty_box.setCurrentIndex(0)
+        self.qtyItemChanged(0)
+    
+    def qtyItemChanged(self,item):
+        """
+        The quantity item has been changed
+        """
+        
+        if self._suppressed: return
+        try:
+            data,datatype = self.getCurrentData()
+        except:
+            data = None
+        
+        segitemno = self._seg_box.currentText()
+        
+        if data is None:
+            self._vstime.setChecked(False)
+            self._vspos.setEnabled(False)
+            self._vstime.setEnabled(False)
+            self._vs_instance.setEnabled(False)
+            self._plot_button.setEnabled(False)
+            self._output_button.setEnabled(False)
+            self._animate_button.setEnabled(False)
+            return
+        
+        self._plot_button.setEnabled(True)
+        self._output_button.setEnabled(True)
+        
+        
+        if datatype == 'time':
+            self._vstime.setChecked(True)
+            self._vspos.setEnabled(False)
+            self._vstime.setEnabled(False)
+            self._vs_instance.setEnabled(False)
+            self._animate_button.setEnabled(False)
+            
+        elif datatype == 'pos':
+            self._vspos.setChecked(True)
+            self._vspos.setEnabled(False)
+            self._vstime.setEnabled(False)
+            self._vs_instance.setEnabled(False)
+            self._animate_button.setEnabled(False)
+            
+        elif datatype == 'postime':
+            vspos = self._vspos.isChecked()
+            if vspos:
+                t = self._file.attrs['t']
+                self._vs_instance.setRange(0,len(t)-1)
+            else:
+                x = self._file[segitemno]['x']                
+                self._vs_instance.setRange(0,len(x)-1)
+            self._vs_instance.setValue(0)
+            
+            self._vspos.setEnabled(True)
+            self._vstime.setEnabled(True)
+            self._vs_instance.setEnabled(True)
+        
+    def initGui(self):
+        self.defaultOutput()
+        self._suppressed = True
+        for x in self._file:
+            if isinstance(self._file[x], h5py.Group):
+                self._seg_box.addItem(x)
+            else:
+                print(self._file[x])
+        self._suppressed = False
+        self.segItemChanged(0)
+
+    def loadFile(self, file=None):
+        pass
+
+    def getx(self):
+        """
+        Returns the x-axis for the current data selected. Returns None if invalid.
+        Returns a tuple containing: x-values, x-legend
+        """
+        try:
+            data,datatype = self.getCurrentData()
+        except:
+            return None
+        
+        segitemno = self._seg_box.currentText()
+        vspos = self._vspos.isChecked()
+        
+        x = None
+        if datatype == 'time' or (datatype == 'postime' and not vspos):
+            x = self._file.attrs['t']
+            xsymbol = 't'
+            xunit = 's'
+        elif datatype == 'pos' or (datatype == 'postime' and vspos):
+            xset = self._file[segitemno]['x']
+            x = xset[:]
+            xsymbol = xset.attrs['symbol'].decode('utf-8')
+            xunit = xset.attrs['unit'].decode('utf-8')
+        
+        xlabel = xsymbol + ' [' + xunit + ']'
+        
+        return x,xlabel
+        
+    def gety(self):
+        """
+        Returns the y-axis for the current data selected. Returns None if invalid.
+        Returns a tuple containing x-values, x-legend
+        """
+        try:
+            data,datatype = self.getCurrentData()
+        except:
+            return None
+        
+        symbol = data.attrs['symbol'].decode('utf-8')
+        unit = data.attrs['unit'].decode('utf-8')
+        
+        vspos = self._vspos.isChecked()
+                                                       
+        ylabel = symbol + ' [' + unit + ']'
+                            
+        if datatype == 'time' or datatype == 'pos':
+            y = data[:]
+        elif datatype == 'postime':
+            data_inst  = self._vs_instance.value()
+            if self._vspos.isChecked():
+                y = data[:,data_inst]
+            else:
+                y = data[data_inst,:]
+                
+        return y,ylabel
+        
+    
+    def plot(self):
+
+        plot = {}
+
+        
+        hor_offset = float(self._hor_offset.text())
+        ver_offset = float(self._ver_offset.text())
+        try:
+            x,xlabel = self.getx()
+        except TypeError:
+            return
+        
+        try:
+            y,ylabel = self.gety()
+        except TypeError:
+            return
+        
+        plot['x'] = x+hor_offset
+        plot['xlabel'] = xlabel
+            
+        thescale = scale[self._scale.currentIndex()]
+        scaleval = thescale[1]
+        scaleylabel = '' if scaleval == 0 else ' x $10^{%i}$' %scaleval
+        
+        plot['y'] = (y+ver_offset)/(10**scaleval)
+        plot['ylabel'] = ylabel + ' ' + scaleylabel
+        
+        self._figure.addPlot(plot, not self._hold_button.isChecked())
+        
+    def about(self):
+        dialog = QDialog(self)
+        about_dialog = Ui_about_dialog()
+        about_dialog.setupUi(dialog)
+        dialog.exec_()
+        
+    def closeEvent(self,event):
+        settings = getGuiSettings()
+        settings.setValue("Geometry", self.saveGeometry())
+        
+        event.accept()
+
+def usage():
+    print("Usage:")
+    if hasattr(sys,'frozen'):
+        pass
+    else:
+        print('python ' + sys.argv[0] + ' <file>.tasmet.h5')
+
+if __name__ == '__main__':
+    if(len(sys.argv) < 2):
+        usage()
+        file_ = '../second_duct_sinomgt.tasmet.h5'
+#        return -1
+    else:
+        file_ = sys.argv[1]
+        
+    # Every PyQt4 application must create an application object.
+    # The application object is located in the QtGui module.
+#    try:
+#        app = QApplication.instance()
+#    except:
+    app = QApplication(sys.argv)
+    
+    # The QWidget widget is the base class of all user interface objects in PyQt4.
+    # We provide the default constructor for QWidget. The default constructor has no parent.
+    # A widget with no parent is called a window. 
+#    w = MainWindow(sys.argv[1])
+    w = MainWindow(file_)    
+    w.resize(640, 480)  # The resize() method resizes the widget.
+    
+    w.show()  # The show() method displays the widget on the screen.
+    
+    sys.exit(app.exec_())  # Finally, we enter the mainloop of the application.