FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_compressibleeuler1.m File Reference

## 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.

Ref. 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 );

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

%