#include <ShockTube.h>

Collaboration diagram for ShockTube< T >:

Public Types
typedef Vector< T, 3 >	Vec3

Public Member Functions
	ShockTube (const int nCells)

T	getPressure (const Vec3 &q)

Vec3	getConservedVars (const T rho, const T p, const T u)

Vec3	getFlux (const Vec3 &q)

Vec3	getRoeFlux_Frink (const Vec3 &q0, const Vec3 &q1)

Vec3	getRoeFlux (const Vec3 &q0, const Vec3 &q1)

void	init ()

void	solve (const double dt, const int nsteps)

ArrayBase &	getSolution ()

ArrayBase &	getPressure ()

ArrayBase &	getDensity ()

ArrayBase &	getVelocity ()

ArrayBase &	getMachNumber ()

Public Attributes
T	pL

T	pR

T	rhoL

T	rhoR

T	uL

T	uR

const int	_nCells

T	_gmm

const int	_nStages

Array< double >	_stageCoeffs

Array< Vec3 >	q

Array< T >	pressure

Array< T >	rho

Array< T >	u

Array< T >	M

Detailed Description

template<class T>
class ShockTube< T >

Definition at line 12 of file ShockTube.h.

Member Typedef Documentation

template<class T >

typedef Vector<T,3> ShockTube< T >::Vec3

Definition at line 16 of file ShockTube.h.

Constructor & Destructor Documentation

template<class T >

ShockTube< T >::ShockTube ( const int nCells )

inline

Definition at line 18 of file ShockTube.h.

References ShockTube< T >::_stageCoeffs.

                               :
     pL(10.0),
     pR(1.0),
     rhoL(8.0),
     rhoR(1.0),
     uL(0),
     uR(0),
     _nCells(nCells),
     _gmm(0.4),
     _nStages(5),
     _stageCoeffs(5),
     q(_nCells),
     pressure(_nCells),
     rho(nCells),
     u(nCells),
     M(nCells)
   {
     _stageCoeffs[0] = 0.25;
     _stageCoeffs[1] = 1.0/6.0;
     _stageCoeffs[2] = 0.375;
     _stageCoeffs[3] = 0.5;
     _stageCoeffs[4] = 1;
   }

Member Function Documentation

template<class T >

Vec3 ShockTube< T >::getConservedVars	(	const T	rho,
		const T	p,
		const T	u
	)

inline

Definition at line 48 of file ShockTube.h.

References ShockTube< T >::_gmm, ShockTube< T >::q, ShockTube< T >::rho, and ShockTube< T >::u.

Referenced by ShockTube< T >::init().

   {
     Vec3 q;
     q[0] = rho;
     q[1] = rho*u;
     q[2] = p/_gmm + 0.5*rho*u*u;
     return q;
   }

template<class T >

ArrayBase& ShockTube< T >::getDensity ( )

inline

Definition at line 234 of file ShockTube.h.

References ShockTube< T >::rho.

234 {return rho;}

ShockTube::rho

Array< T > rho

Definition: ShockTube.h:252

template<class T >

Vec3 ShockTube< T >::getFlux ( const Vec3 & q )

inline

Definition at line 57 of file ShockTube.h.

References ShockTube< T >::getPressure(), and ShockTube< T >::u.

Referenced by ShockTube< T >::getRoeFlux(), and ShockTube< T >::getRoeFlux_Frink().

   {
     const T u = q[1]/q[0];
     const T p = getPressure(q);
     Vec3 flux;
     flux[0] = q[1];
     flux[1] = q[0]*u*u + p;
     flux[2] = (q[2]+p)*u;
     return flux;
   }

template<class T >

ArrayBase& ShockTube< T >::getMachNumber ( )

inline

Definition at line 236 of file ShockTube.h.

References ShockTube< T >::M.

236 {return M;}

ShockTube::M

Array< T > M

Definition: ShockTube.h:254

template<class T >

T ShockTube< T >::getPressure ( const Vec3 & q )

inline

Definition at line 42 of file ShockTube.h.

References ShockTube< T >::_gmm, and ShockTube< T >::u.

   {
     const T u = q[1]/q[0];
     return (q[2] - 0.5*(q[0]*u*u))*_gmm;
   }

template<class T >

ArrayBase& ShockTube< T >::getPressure ( )

inline

Definition at line 233 of file ShockTube.h.

References ShockTube< T >::pressure.

Referenced by ShockTube< T >::getFlux(), ShockTube< T >::getRoeFlux(), ShockTube< T >::getRoeFlux_Frink(), and ShockTube< T >::solve().

233 {return pressure;}

ShockTube::pressure

Array< T > pressure

Definition: ShockTube.h:251

template<class T >

Vec3 ShockTube< T >::getRoeFlux	(	const Vec3 &	q0,
		const Vec3 &	q1
	)

inline

Definition at line 120 of file ShockTube.h.

References ShockTube< T >::_gmm, fabs(), ShockTube< T >::getFlux(), ShockTube< T >::getPressure(), and sqrt().

Referenced by ShockTube< T >::solve().

   {
     const T half(0.5);
     
     const T sr0 = sqrt(q0[0]);
     const T sr1 = sqrt(q1[0]);
 
     const T sr0psr1 = sr0 + sr1;
     const T u0 = q0[1]/q0[0];
     const T u1 = q1[1]/q1[0];
 
     const T p0 = getPressure(q0);
     const T p1 = getPressure(q1);
 
     const T h0 = (q0[2]+p0) / q0[0];
     const T h1 = (q1[2]+p1) / q1[0];
     
     const T rho_f = sr0 * sr1;
     const T u_f = (sr0*u0 + sr1*u1) / sr0psr1;
     const T h_f = (sr0*h0 + sr1*h1) / sr0psr1;
 
     const T c_f = sqrt( _gmm * (h_f - half*u_f*u_f));
 
     const T abs_uf = fabs(u_f);
     const T abs_upc = fabs(u_f + c_f);
     const T abs_umc = fabs(u_f - c_f);
 
     
     const Vec3 flux0 = getFlux(q0);
     const Vec3 flux1 = getFlux(q1);
 
     const T dr = q1[0]-q0[0];
     const T rdu = rho_f*(u1-u0);
 
     const T dp = p1-p0;
 
     const T dpbyc2 = dp/(c_f*c_f);
 
     const T a0 = half*(dpbyc2 + rdu/c_f)*abs_upc;
     const T a4 = half*(dpbyc2 - rdu/c_f)*abs_umc;
 
     const T drmdpbyc2 = dr-dpbyc2;
     Vec3 dFlux;
 
     dFlux[0] = abs_uf*drmdpbyc2 + a0 + a4;
     dFlux[1] = abs_uf*drmdpbyc2*u_f + a0*(u_f+c_f) + a4*(u_f-c_f);
     dFlux[2] = abs_uf*(half*drmdpbyc2*u_f*u_f) + a0*(h_f+u_f*c_f)
       + a4*(h_f-u_f*c_f);
     
     Vec3 flux(flux0+flux1-dFlux);
     return half*flux;
   }

template<class T >

Vec3 ShockTube< T >::getRoeFlux_Frink	(	const Vec3 &	q0,
		const Vec3 &	q1
	)

inline

Definition at line 68 of file ShockTube.h.

References ShockTube< T >::_gmm, fabs(), ShockTube< T >::getFlux(), ShockTube< T >::getPressure(), and sqrt().

   {
     const T half(0.5);
     
     const T sr0 = sqrt(q0[0]);
     const T sr1 = sqrt(q1[0]);
 
     const T sr0psr1 = sr0 + sr1;
     const T u0 = q0[1]/q0[0];
     const T u1 = q1[1]/q1[0];
 
     const T p0 = getPressure(q0);
     const T p1 = getPressure(q1);
 
     const T h0 = (q0[2]+p0) / q0[0];
     const T h1 = (q1[2]+p1) / q1[0];
     
     const T rho_f = sr0 * sr1;
     const T u_f = (sr0*u0 + sr1*u1) / sr0psr1;
     const T h_f = (sr0*h0 + sr1*h1) / sr0psr1;
 
     const T c_f = sqrt( _gmm * (h_f - half*u_f*u_f));
 
     const T abs_uf = fabs(u_f);
     const T abs_upc = fabs(u_f + c_f);
     const T abs_umc = fabs(u_f - c_f);
 
     
     const Vec3 flux0 = getFlux(q0);
     const Vec3 flux1 = getFlux(q1);
 
     const T dr = q1[0]-q0[0];
     const T rdu = rho_f*(u1-u0);
 
     const T dp = p1-p0;
 
     const T dpc2 = dp*c_f*c_f;
 
     const T a0 = half*(dpc2 + rdu*c_f)*abs_upc;
     const T a4 = half*(dpc2 - rdu*c_f)*abs_umc;
 
     const T drmdpc2 = dr-dpc2;
     Vec3 dFlux;
 
     dFlux[0] = abs_uf*drmdpc2 + a0 + a4;
     dFlux[1] = abs_uf*drmdpc2*u_f + a0*(u_f+c_f) + a4*(u_f-c_f);
     dFlux[2] = abs_uf*(half*drmdpc2*u_f*u_f) + a0*(h_f+u_f*c_f) + a4*(h_f-u_f*c_f);
     
     Vec3 flux(flux0+flux1-dFlux);
     return half*flux;
   }

template<class T >

ArrayBase& ShockTube< T >::getSolution ( )

inline

Definition at line 232 of file ShockTube.h.

References ShockTube< T >::q.

232 {return q;}

ShockTube::q

Array< Vec3 > q

Definition: ShockTube.h:250

template<class T >

ArrayBase& ShockTube< T >::getVelocity ( )

inline

Definition at line 235 of file ShockTube.h.

References ShockTube< T >::u.

235 {return u;}

ShockTube::u

Array< T > u

Definition: ShockTube.h:253

template<class T >

void ShockTube< T >::init ( )

inline

Definition at line 174 of file ShockTube.h.

References ShockTube< T >::_nCells, ShockTube< T >::getConservedVars(), ShockTube< T >::pL, ShockTube< T >::pR, ShockTube< T >::q, ShockTube< T >::rhoL, ShockTube< T >::rhoR, ShockTube< T >::uL, and ShockTube< T >::uR.

   {
     Vec3 qL(getConservedVars(rhoL, pL, uL));
     Vec3 qR(getConservedVars(rhoR, pR, uR));
 
     for(int i=0; i<_nCells/2; i++)
       q[i] = qL;
     
     for(int i=_nCells/2; i<_nCells; i++)
       q[i] = qR;
   }

template<class T >

void ShockTube< T >::solve	(	const double	dt,
		const int	nsteps
	)

inline

Definition at line 186 of file ShockTube.h.

References ShockTube< T >::_gmm, ShockTube< T >::_nCells, ShockTube< T >::_nStages, ShockTube< T >::_stageCoeffs, ShockTube< T >::getPressure(), ShockTube< T >::getRoeFlux(), ShockTube< T >::M, ShockTube< T >::pressure, ShockTube< T >::q, ShockTube< T >::rho, sqrt(), ShockTube< T >::u, and Array< T >::zero().

   {
     const T dx = T(1.0)/_nCells;
 
     Array<Vec3> R(_nCells);
 
     Array<Vec3> qN(_nCells);
 
     for(int n=0; n<nsteps; n++)
     {
         cout << "step " << n << endl;
         qN = q;
         for(int s=0; s<_nStages; s++)
         {
             R.zero();
 
             for(int f=1; f<_nCells; f++)
             {
                 int c0 = f-1;
                 int c1 = f;
                 const Vec3 flux = getRoeFlux(q[c0], q[c1]);
                 R[c0] += flux;
                 R[c1] -= flux;
             }
 
             R[0] -= getRoeFlux(q[0],q[0]);
             R[_nCells-1] += getRoeFlux(q[_nCells-1],q[_nCells-1]);
             
             const T sCoeff = _stageCoeffs[s]*dt/dx;
             for(int c=0; c<_nCells; c++)
             {
                 q[c] = qN[c] - sCoeff*R[c];
             }
         }
     }
 
     for(int c=0; c<_nCells; c++)
     {
         rho[c] = q[c][0];
         pressure[c] = getPressure(q[c]);
         u[c] = q[c][1]/q[c][0];
         T soundSpeed = sqrt((_gmm+1)*pressure[c]/rho[c]);
         M[c] = u[c]/soundSpeed;
     }
   }