Hello Daniel,
Thanks for the answer
An implementation of the Roe scheme would certainly go a long way in
simplifying the resolution of hyperbolic equations but I think FiPy already
has was it takes to do that - terms wise.
I know of CLAWPACK yes, but I really would like to see if FiPy can handle
fluid dynamics equations. I work at CERFACS, a lab in France centered
around the use of Computational Fluid Dynamics for academic and industrial
purpose, and although we do have a production code, I would like to let my
colleagues know about FiPy because I think it has the potential to be a
great first-approach solver for simple to mildly complex problems.
I actually already succeeded in solving (and validating with the known Sod
shock tube test case) the 1D compressible Euler equations (script attached
with an exact solver for reference).
But then in 1D, things are pretty easy - equation wise.
If you have time to check the code, can you tell me if you see any kind of
mis-use of FiPy ? Are there any lines you would have written differently ?
In 2D everything becomes harder and I can’t get it right - although I
suspect it is because I don’t know FiPy very well yet.
If you don’t mind, I have a couple questions to try and make it work.
1/ I do not get how to include the resolution of an equation of state in
the system of differential equations. I wrote a stub of code (attached,
called 7-covo2D.py) in which I solve for (rho, rhoU, rhoV, rhoE) as usual.
An easy way to write the Euler equations however is to include a notion of
pressure, which is related to the unknown by an equation of state (it is
not a differential equation though) - see for instance
http://bender.astro.sunysb.edu/hydro_by_example/compressible/Euler.pdf.
Because of this, I have to inject the EoS directly into the equations and
try and make do with that, but the Finite Volume method is not designed to
have the fluxes computed like I do in the script, so naturally everything
crashes or gives nonsensical results.
Do you see a possibility in FiPy to solve 4 differential equations and 1
algebraic equation (the EoS) in a coupled manner ?
2/ If it is not possible to include the equation of state in the system,
how would you take the -dp/dx and the -dp/dy from the momentum equations
(with p expanded as a function of rhoE, rhoU and rhoV) into account ?
Thank you very much for your help,
Best regards,
Thibault
2017-05-15 15:38 GMT+02:00 Daniel Wheeler <[email protected]>:
> Hi Thibault,
>
> I started down the road of implementing a Riemann solver in FiPy, see
>
> https://github.com/usnistgov/fipy/blob/riemann/fipy/terms/
> baseRoeConvectionTerm.py
>
> Unfortunately, I never completed and merged this work. You might want
> to look into CLAWPACK for a code to better handle compressible flow.
>
> I've given a few answers about this in the past as well, see
>
> - http://git.net/ml/python-fipy/2012-02/msg00024.html
>
> - http://fipy.nist.narkive.com/A0gJrSl2/semilinear-wave-
> equation-in-fipy
>
> - http://fipy.nist.narkive.com/YVZTRM0G/1d-coupled-fluid-
> equations-in-fipy
>
> Cheers,
>
> Daniel
>
> On Sun, May 14, 2017 at 9:33 AM, Thibault Bridel-Bertomeu
> <[email protected]> wrote:
> > Good afternoon,
> >
> > I was wondering if anyone had ever tried implementing the compressible
> euler
> > equations in FiPy ?
> > I have been trying for a while now but even in 1D I can’t get a proper
> shock
> > tube solution.
> >
> > If anyone has ever tried, I would be infinitely grateful if they were to
> > share their knowledge !!
> >
> > Thank you very much,
> >
> > T. BRIDEL-BERTOMEU
> >
> >
> >
> > _______________________________________________
> > fipy mailing list
> > [email protected]
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
>
>
>
> --
> Daniel Wheeler
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
#!/usr/bin/env python
import fipy
from fipy import numerix
#==============================
nx = 1000
dx = 1./nx
Lx = nx * dx
gamma = 1.4
gm1 = gamma-1.0
#------------------------------
mesh = fipy.Grid1D(nx=nx, dx=dx)
X = mesh.cellCenters[0]
x = mesh.faceCenters[0]
#------------------------------
ro = fipy.CellVariable(mesh=mesh, name="ro", hasOld=True)
roU = fipy.CellVariable(mesh=mesh, name="roU", hasOld=True)
roE = fipy.CellVariable(mesh=mesh, name="roE", hasOld=True)
ro.setValue( 1.0*(X<=0.5)+0.125*(X>0.5) )
roU.setValue( 0.0 )
roE.setValue( 1.0/gm1*(X<=0.5)+0.1/gm1*(X>0.5) )
ro.constrain(value=1.0, where=mesh.facesLeft)
ro.constrain(value=0.125, where=mesh.facesRight)
roU.constrain(value=0.0, where=mesh.facesLeft|mesh.facesRight)
roE.constrain(value=1.0/gm1, where=mesh.facesLeft)
roE.constrain(value=0.1/gm1, where=mesh.facesRight)
vi = fipy.Viewer(vars=(ro, roU, roE), datamin=-0.1, datamax=2.6)
vi.plot()
raw_input("Initialization ...")
U = fipy.CellVariable(mesh=mesh, name="U")
U.setValue( roU / ro )
p = fipy.CellVariable(mesh=mesh, name="p")
p.setValue( gm1*(roE - 0.5*roU*U) )
vi = fipy.Viewer(vars=(ro, U, p), datamin=-0.1, datamax=1.1, )
vi.plot()
raw_input("Primitive variables ...")
#------------------------------
eqnC = fipy.TransientTerm(var=ro) + fipy.ExponentialConvectionTerm(coeff=(roU[None,:]/ro[None,:]).faceValue, var=ro) == 0.0
eqnM = fipy.TransientTerm(var=roU) + \
(1.0-(gamma-1.0)/2.0)*fipy.ExponentialConvectionTerm(coeff=(roU[None,:]/ro[None,:]).faceValue, var=roU) + \
(gamma-1.0)*fipy.CentralDifferenceConvectionTerm(coeff=((1.,),), var=roE) == 0.0
# (gamma-1.0)*roE.grad[0] == 0.0
eqnE = fipy.TransientTerm(var=roE) + \
fipy.ExponentialConvectionTerm(coeff=(roU[None,:]/ro[None,:]).faceValue, var=roE) + \
(gamma-1.0)*fipy.ExponentialConvectionTerm(coeff=(roU[None,:]/ro[None,:]).faceValue, var=roE) - \
0.5*(gamma-1.0)*fipy.ExponentialConvectionTerm(coeff=(roU[None,:]/ro[None,:]).faceValue, var=roU) == 0.0
eqn = eqnC & eqnM & eqnE
#------------------------------
umax = 1.0
CFL = 0.75
dt = dx * CFL / umax
tmax = 0.2
steps = int(tmax / dt)
for step in range(steps):
ro.updateOld()
roU.updateOld()
roE.updateOld()
#
eqn.solve(dt=dt)
#
U.setValue( roU / ro )
p.setValue( gm1*(roE - 0.5*roU*U) )
vi.plot()
raw_input("Final state ...")
#------------------------------
import sod
positions, regions, values = sod.solve(left_state=(1, 1, 0), right_state=(0.1, 0.125, 0.),
geometry=(0., 1., 0.5), t=tmax, gamma=gamma, npts=nx+1)
press_exact = fipy.CellVariable(mesh=mesh, value=(values['p'][:-1] + values['p'][1:])/2.0, name="p_exact")
rho_exact = fipy.CellVariable(mesh=mesh, value=(values['rho'][:-1] + values['rho'][1:])/2.0, name="rho_exact")
u_exact = fipy.CellVariable(mesh=mesh, value=(values['u'][:-1] + values['u'][1:])/2.0, name="u_exact")
rhovi = fipy.Viewer(vars=(ro, rho_exact), datamin=-0.1, datamax=1.1, limits={'xmin':0.2, 'xmax':0.9})
pvi = fipy.Viewer(vars=(p, press_exact), datamin=-0.1, datamax=1.1, limits={'xmin':0.2, 'xmax':0.9})
Uvi = fipy.Viewer(vars=(U, u_exact), datamin=-0.1, datamax=1.1, limits={'xmin':0.2, 'xmax':0.9})
rhovi.plot()
pvi.plot()
Uvi.plot()
raw_input("Comparison with exact solution ...")
import numpy as np
import scipy
import scipy.optimize
def shock_tube_function(p4, p1, p5, rho1, rho5, gamma):
"""
Shock tube equation
"""
z = (p4 / p5 - 1.)
c1 = np.sqrt(gamma * p1 / rho1)
c5 = np.sqrt(gamma * p5 / rho5)
gm1 = gamma - 1.
gp1 = gamma + 1.
g2 = 2. * gamma
fact = gm1 / g2 * (c5 / c1) * z / np.sqrt(1. + gp1 / g2 * z)
fact = (1. - fact) ** (g2 / gm1)
return p1 * fact - p4
def calculate_regions(pl, ul, rhol, pr, ur, rhor, gamma=1.4):
"""
Compute regions
:rtype : tuple
:return: returns p, rho and u for regions 1,3,4,5 as well as the shock speed
"""
# if pl > pr...
rho1 = rhol
p1 = pl
u1 = ul
rho5 = rhor
p5 = pr
u5 = ur
# unless...
if pl < pr:
rho1 = rhor
p1 = pr
u1 = ur
rho5 = rhol
p5 = pl
u5 = ul
# solve for post-shock pressure
p4 = scipy.optimize.fsolve(shock_tube_function, p1, (p1, p5, rho1, rho5, gamma))[0]
# compute post-shock density and velocity
z = (p4 / p5 - 1.)
c5 = np.sqrt(gamma * p5 / rho5)
gm1 = gamma - 1.
gp1 = gamma + 1.
gmfac1 = 0.5 * gm1 / gamma
gmfac2 = 0.5 * gp1 / gamma
fact = np.sqrt(1. + gmfac2 * z)
u4 = c5 * z / (gamma * fact)
rho4 = rho5 * (1. + gmfac2 * z) / (1. + gmfac1 * z)
# shock speed
w = c5 * fact
# compute values at foot of rarefaction
p3 = p4
u3 = u4
rho3 = rho1 * (p3 / p1)**(1. / gamma)
return (p1, rho1, u1), (p3, rho3, u3), (p4, rho4, u4), (p5, rho5, u5), w
def calc_positions(pl, pr, region1, region3, w, xi, t, gamma):
"""
:return: tuple of positions in the following order ->
Head of Rarefaction: xhd, Foot of Rarefaction: xft,
Contact Discontinuity: xcd, Shock: xsh
"""
p1, rho1 = region1[:2] # don't need velocity
p3, rho3, u3 = region3
c1 = np.sqrt(gamma * p1 / rho1)
c3 = np.sqrt(gamma * p3 / rho3)
if pl > pr:
xsh = xi + w * t
xcd = xi + u3 * t
xft = xi + (u3 - c3) * t
xhd = xi - c1 * t
else:
# pr > pl
xsh = xi - w * t
xcd = xi - u3 * t
xft = xi - (u3 - c3) * t
xhd = xi + c1 * t
return xhd, xft, xcd, xsh
def region_states(pl, pr, region1, region3, region4, region5):
"""
:return: dictionary (region no.: p, rho, u), except for rarefaction region
where the value is a string, obviously
"""
if pl > pr:
return {'Region 1': region1,
'Region 2': 'RAREFACTION',
'Region 3': region3,
'Region 4': region4,
'Region 5': region5}
else:
return {'Region 1': region5,
'Region 2': region4,
'Region 3': region3,
'Region 4': 'RAREFACTION',
'Region 5': region1}
def create_arrays(pl, pr, xl, xr, positions, state1, state3, state4, state5, npts, gamma, t, xi):
"""
:return: tuple of x, p, rho and u values across the domain of interest
"""
xhd, xft, xcd, xsh = positions
p1, rho1, u1 = state1
p3, rho3, u3 = state3
p4, rho4, u4 = state4
p5, rho5, u5 = state5
gm1 = gamma - 1.
gp1 = gamma + 1.
dx = (xr-xl)/(npts-1)
x_arr = np.arange(xl, xr, dx)
rho = np.zeros(npts, dtype=float)
p = np.zeros(npts, dtype=float)
u = np.zeros(npts, dtype=float)
c1 = np.sqrt(gamma * p1 / rho1)
if pl > pr:
for i, x in enumerate(x_arr):
if x < xhd:
rho[i] = rho1
p[i] = p1
u[i] = u1
elif x < xft:
u[i] = 2. / gp1 * (c1 + (x - xi) / t)
fact = 1. - 0.5 * gm1 * u[i] / c1
rho[i] = rho1 * fact ** (2. / gm1)
p[i] = p1 * fact ** (2. * gamma / gm1)
elif x < xcd:
rho[i] = rho3
p[i] = p3
u[i] = u3
elif x < xsh:
rho[i] = rho4
p[i] = p4
u[i] = u4
else:
rho[i] = rho5
p[i] = p5
u[i] = u5
else:
for i, x in enumerate(x_arr):
if x < xsh:
rho[i] = rho5
p[i] = p5
u[i] = -u1
elif x < xcd:
rho[i] = rho4
p[i] = p4
u[i] = -u4
elif x < xft:
rho[i] = rho3
p[i] = p3
u[i] = -u3
elif x < xhd:
u[i] = -2. / gp1 * (c1 + (xi - x) / t)
fact = 1. + 0.5 * gm1 * u[i] / c1
rho[i] = rho1 * fact ** (2. / gm1)
p[i] = p1 * fact ** (2. * gamma / gm1)
else:
rho[i] = rho1
p[i] = p1
u[i] = -u1
return x_arr, p, rho, u
def solve(left_state, right_state, geometry, t, gamma=1.4, npts=500):
"""
Solves the Sod shock tube problem (i.e. riemann problem) of discontinuity across an interface.
:rtype : tuple
:param left_state: tuple (pl, rhol, ul)
:param right_state: tuple (pr, rhor, ur)
:param geometry: tuple (xl, xr, xi): xl - left boundary, xr - right boundary, xi - initial discontinuity
:param t: time for which the states have to be calculated
:param gamma: ideal gas constant, default is air: 1.4
:param npts: number of points for array of pressure, density and velocity
:return: tuple of: dicts of positions,
constant pressure, density and velocity states in distinct regions,
arrays of pressure, density and velocity in domain bounded by xl, xr
"""
pl, rhol, ul = left_state
pr, rhor, ur = right_state
xl, xr, xi = geometry
# basic checking
if xl >= xr:
print('xl has to be less than xr!')
exit()
if xi >= xr or xi <= xl:
print('xi has in between xl and xr!')
exit()
# calculate regions
region1, region3, region4, region5, w = \
calculate_regions(pl, ul, rhol, pr, ur, rhor, gamma)
regions = region_states(pl, pr, region1, region3, region4, region5)
# calculate positions
x_positions = calc_positions(pl, pr, region1, region3, w, xi, t, gamma)
pos_description = ('Head of Rarefaction', 'Foot of Rarefaction',
'Contact Discontinuity', 'Shock')
positions = dict(zip(pos_description, x_positions))
# create arrays
x, p, rho, u = create_arrays(pl, pr, xl, xr, x_positions,
region1, region3, region4, region5,
npts, gamma, t, xi)
val_names = ('x', 'p', 'rho', 'u')
val_dict = dict(zip(val_names, (x, p, rho, u)))
return positions, regions, val_dict#!/usr/bin/env python
import fipy
from fipy import numerix
#==============================
nx = 64
ny = 64
dx = 0.25
dy = 0.25
Lx = nx*dx
Ly = ny*dy
mesh = fipy.PeriodicGrid2D(nx=nx, ny=ny, dx=dx, dy=dy)
X, Y = mesh.cellCenters
#==============================
gamma = 1.4
gm1 = gamma - 1.0
#==============================
ro = fipy.CellVariable(mesh=mesh, name="ro", hasOld=True)
roU = fipy.CellVariable(mesh=mesh, name="roU", hasOld=True)
roV = fipy.CellVariable(mesh=mesh, name="roV", hasOld=True)
roE = fipy.CellVariable(mesh=mesh, name="roE", hasOld=True)
p = fipy.CellVariable(mesh=mesh, name="p")
X0 = 0.5*Lx
Y0 = 0.5*Ly
Rsq = (X-X0)**2 + (Y-Y0)**2
uinf = 1.0
vinf = 0.0
pinf = 1.0
roinf = 1.0
beta = 5.0
ro_init = (pinf/roinf + gm1*beta**2/(8*gamma*numerix.pi)*numerix.exp(1.0-Rsq))**(1/gm1)
roU_init = ro_init * (uinf - beta*numerix.exp(0.5*(1-Rsq))*0.5/numerix.pi)
roV_init = ro_init * (vinf + beta*numerix.exp(0.5*(1-Rsq))*0.5/numerix.pi)
press_init = ro_init**gamma
roE_init = press_init / gm1 + 0.5 * ro_init * ( roU_init**2 + roV_init**2 ) / ro_init**2
ro.setValue( ro_init )
roU.setValue( roU_init )
roV.setValue( roV_init )
roE.setValue( roE_init )
vi = fipy.Viewer(vars=roU)
vi.plot()
raw_input("Initialization. Press <return> to continue ...")
#==============================
coeffConv = fipy.CellVariable(mesh=mesh, rank=1)
coeffConv[0] = roU/ro
coeffConv[1] = roV/ro
coeffConvMinus = fipy.CellVariable(mesh=mesh, rank=1)
coeffConvMinus[0] = -roV/ro
coeffConvMinus[1] = roU/ro
coeffConvSquare = fipy.CellVariable(mesh=mesh, rank=1)
coeffConvSquare[0] = roU**2/ro**2
coeffConvSquare[1] = roU*roV/ro**2
coeffConvSquareBis = fipy.CellVariable(mesh=mesh, rank=1)
coeffConvSquareBis[0] = roU*roV/ro**2
coeffConvSquareBis[1] = roV**2/ro**2
#----
eqnC = fipy.TransientTerm(var=ro) + \
fipy.UpwindConvectionTerm(coeff=coeffConv.faceValue, var=ro) == 0.0
#----
eqnMx = fipy.TransientTerm(var=roU) + \
(1.0-0.5*gm1)*fipy.UpwindConvectionTerm(coeff=coeffConv.faceValue, var=roU) + \
gm1*roE.grad.dot([1., 0.]) + \
0.5*gm1*fipy.UpwindConvectionTerm(coeff=coeffConvMinus.faceValue, var=roV) == 0.0
#----
eqnMy = fipy.TransientTerm(var=roV) + \
(1.0-0.5*gm1)*fipy.UpwindConvectionTerm(coeff=coeffConv.faceValue, var=roV) + \
gm1*roE.grad.dot([0., 1.]) - \
0.5*gm1*fipy.UpwindConvectionTerm(coeff=coeffConvMinus.faceValue, var=roU) == 0.0
#----
eqnE = fipy.TransientTerm(var=roE) + \
(1.0+gm1)*fipy.UpwindConvectionTerm(coeff=coeffConv.faceValue, var=roE) - \
0.5*gm1*fipy.UpwindConvectionTerm(coeff=coeffConvSquare.faceValue, var=roU) - \
0.5*gm1*fipy.UpwindConvectionTerm(coeff=coeffConvSquareBis.faceValue, var=roV) == 0.0
#----
eqn = eqnC & eqnMx & eqnMy & eqnE
#==============================
umax = max(max(roU_init/ro_init), max(roV_init/ro_init))
CFL = 0.2
dx = numerix.sqrt(min(mesh.cellVolumes))
dt = dx * CFL / umax
print 'TIMESTEP =', dt
steps = 10000
for step in range(steps):
ro.updateOld()
roU.updateOld()
roV.updateOld()
roE.updateOld()
#
p.setValue( gm1*(roE - 0.5*(roU**2 + roV**2)/ro ) )
#
eqn.solve(dt=dt)
#
vi.plot()
raw_input("Final state ...")
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
