Description

EX_COMPRESSIBLEEULER1 SOD Shock tube example.

[ FEA, OUT ] = EX_COMPRESSIBLEEULER1( VARARGIN ) Sets up and solves an instationary SOD shock tube compressible Euler equation example on the unit line. Initial data are rhoL, uL=0, pL=1 and rhoR=0.125, uR=0, pR=0.1 with an inital discontinuity at t=0 and x=0.5. The computed solution with shocks, expansion, and compression waves is compared with the analytical Riemann solution.

Reference

[1] G.A. Sod, A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws, JCP, 1978.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.002}         Max grid cell size
tmax        scalar {0.2}           Maximum time
tstep       scalar {0.0025}        Time step size
ischeme     scalar {2}             Time discretization scheme
sfun        string {sflag1}        Shape function
iplot       scalar 0/{1}           Plot solution and error (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { 'hmax',     0.002;
             'tmax',     0.2;
             'tstep',    0.0025;
             'sfun',     'sflag1';
             'ischeme',  2;
             'iplot',    1;
             'tol',      0.05;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


 fea.sdim = { 'x' };
 fea.geom.objects = { gobj_line(0,1) };
 fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );

 fea = addphys(fea,@compressibleeuler);

 rhoL  = 1;
 uL    = 0;
 pL    = 1;

 rhoR  = 0.125;
 uR    = 0;
 pR    = 0.1;

 x     = fea.grid.p(:);
 x0    = (max(x) - min(x))/2;

 init0 = { sprintf('%g+(x>%g)*%g',rhoL,x0,rhoR-rhoL), ...
           sprintf('%g+(x>%g)*%g',uL,x0,uR-uL), ...
           sprintf('%g+(x>%g)*%g',pL,x0,pR-pL) };

 fea = parsephys(fea);
 fea = parseprob(fea);

 [fea.sol.u,fea.sol.t] = solvetime( fea, 'init', init0, 'fid', fid, ...
                                    'ischeme', 2, 'tstep', opt.tstep, 'tmax', opt.tmax );


% Error checking.
 [r_ref,u_ref, p_ref] = analytic_sod( opt.tmax, x );
 r = evalexprp( fea.dvar{1}, fea );
 u = evalexprp( fea.dvar{2}, fea );
 p = evalexprp( fea.dvar{3}, fea );
 out.err = [ sum(abs(r_ref-r))/length(x), ...
             sum(abs(u_ref-u))/length(x), ...
             sum(abs(p_ref-p))/length(x) ];
 out.pass = all(out.err<opt.tol);


% Postprocessing.
 if( opt.iplot>0 )
   figure
   subplot(2,2,1)
   plot( x, r, 'k-' )
   hold on
   plot( x, r_ref, 'b--' )
   legend('Computed','Analytic solution')
   title('Density (rho)')
   subplot(2,2,2)
   plot( x, u, 'k-' )
   hold on
   plot( x, u_ref, 'b--' )
   title('Velocity (u)')
   subplot(2,2,3)
   plot( x, p, 'k-' )
   hold on
   plot( x, p_ref, 'b--' )
   title('Pressure (p)')
   subplot(2,2,4)
   plot( x, p./r, 'k-' )
   hold on
   plot( x, p_ref./r_ref, 'b--' )
   title('Non-dimensional temperature (p/rho)')
 end


 if( nargout==0 )
   clear fea out
 end


%