Hello Daniel,
Thank you for the paper and the script - I am afraid it will take me some
time though, its impressive and long work !!
Regarding the equations, I think I was not clear enough in my previous
explanations, I apologize.
In substance, I have 4 variables : ro, roU, roV and roE. They are density,
velocity along X, velocity along Y and energy.
I also have 4 differential equations :
dro/dt + nabla.(ro*[U,V]) = 0
droU/dt + nabla.(roU*[U,V]) = -dp / dx
droV/dt + nabla.(roV*[U,V]) = -dp / dy
droE/dt + nabla.(roE*[U,V]) = - d(p*(roU/ro))/dx - d(p*(roV/ro))/dy
As you can see, they are written with a fifth variable, p, for pressure,
that is related to the others by :
p = (gamma-1.0)*roE - 0.5*(gamma-1.0)*(roU**2 + roV**2)/ro
this is what I call the equation of state, it is not differential it is
just algebraic.
I wager the differential equations above are not the most complex you have
seen or implemented in FiPy, and I think I succeeded, although of course,
since the whole thing does not work, I cannot be sure. But then I am stuck
with that non-differential equation that I would have to solve with the
four differential (?) to close the system …
Can you see what I am stuck with ?
Also, could you please elaborate on the last paragraph of your previous
e-mail ? « If you do include … » I am not sure I get why you speak of
linearization and relaxation ?
I attach the latest version of my non-working script ..
Thanks for the help
Best,
Thibault
2017-05-17 17:06 GMT+02:00 Daniel Wheeler <daniel.wheel...@gmail.com>:
> Hi Thibault,
>
> The paper and your attempts at solving it are very impressive, nice
> work. I'm not sure I understand the issue with the equation of state /
> 2D issue though. Is it that you want to couple in the equation of
> state with the other equations instead of solving explicitly?
>
> Here is a link for some work I did with FiPy quite a long time ago,
>
> https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec2d
>
> It resulted in a publication, see https://arxiv.org/pdf/1006.4881.pdf.
> In that work, everything is being solved in one matrix, see
>
> https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec
> 2d#file-input-py-L245
>
> Look at,
>
> https://gist.github.com/wd15/c28ab796cb3d9781482b01fb67a7ec
> 2d#file-input-py-L286
>
> That particular equation is like an equation of state (it has no
> transient term). The point is that with the use of
> ImplicitSourceTerm's, you can include an equation of state implicitly.
> I'm not sure whether this is your issue or not, but it might be
> related. The above set of equations are dissipative, which yours are
> not, but also, I wasn't concerned with resolving the shocks
> numerically so a Roe solver wasn't used. Also, this work was done
> before FiPy had coupled solutions.
>
> If you do include the equation of state in the matrix, then think
> carefully about how you linearize the "roU" and "roV" terms. Also, as
> the equation of state is like an immediate update, it might be
> necessary to relax the updates a bit during the non-linear iteration
> cycle.
>
> Cheers,
>
> Daniel
>
>
>
> On Mon, May 15, 2017 at 11:46 AM, Thibault Bridel-Bertomeu
> <thibault.bridellel...@gmail.com> wrote:
> > 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 <daniel.wheel...@gmail.com>:
> >>
> >> 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
> >> <thibault.bridellel...@gmail.com> 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
> >> > fipy@nist.gov
> >> > http://www.ctcms.nist.gov/fipy
> >> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
> ]
> >> >
> >>
> >>
> >>
> >> --
> >> Daniel Wheeler
> >>
> >> _______________________________________________
> >> fipy mailing list
> >> fipy@nist.gov
> >> http://www.ctcms.nist.gov/fipy
> >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
> >
> >
> > _______________________________________________
> > fipy mailing list
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
>
>
>
> --
> Daniel Wheeler
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> 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
from fipy import parallel
import numpy as np
#==============================
nx = 128
ny = 128
# dx = 0.1
# dy = 0.1
# Lx = nx*dx
# Ly = ny*dy
Lx = 0.1
Ly = 0.1
dx = Lx / nx
dy = Ly / ny
mesh = fipy.PeriodicGrid2D(nx=nx, ny=ny, dx=dx, dy=dy)
X, Y = mesh.cellCenters
x, y = mesh.faceCenters
if parallel.procID == 0:
print Y
print y
#==============================
gamma = 1.4
gm1 = gamma - 1.0
Cv = 717.875
Cp = gamma * Cv
Rgas = Cp - Cv
#==============================
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")#, hasOld=True)
X0 = 0.05
Y0 = 0.05
beta = 5.0
Rsq = (X-X0)**2 + (Y-Y0)**2
Rc = 0.005
pinf = 1.0
roinf = 1.0
uinf = 1.0
vinf = 0.0
# ro_init = (pinf/roinf + gm1*beta**2/(8*gamma*numerix.pi**2)*numerix.exp((1.0-Rsq)/Rc**2))**(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
U_init = uinf - beta/(2*np.pi)*numerix.exp(-0.5*Rsq/Rc**2)*(Y-Y0)/Rc
V_init = vinf + beta/(2*np.pi)*numerix.exp(-0.5*Rsq/Rc**2)*(X-X0)/Rc
temp_init = 1.0 - gm1*beta**2/(8*gamma*np.pi**2)*numerix.exp(-Rsq/Rc**2)
ro_init = temp_init**(1.0/gm1)
press_init = ro_init**gamma
roE_init = press_init / gm1 + 0.5 * ro_init * (U_init**2 + V_init**2)
ro.setValue( ro_init )
roU.setValue( ro_init * U_init )
roV.setValue( ro_init * V_init )
roE.setValue( roE_init )
p.setValue( press_init )
# X0 = 0.05
# Y0 = 0.05
# Gamma = 34.7279426399
# Rc = 0.005
# cut = 8.0
# press_inf = 100000
# Ma_inf = 0.0 #0.5
# temp_inf = 300.0
# alpha_inf = 0.0
# beta_inf = 0.0
# ro_inf = press_inf / (Rgas * temp_inf)
# c_inf = numerix.sqrt(gamma * press_inf / ro_inf)
# unorm_inf = c_inf * Ma_inf
# U_inf = unorm_inf * numerix.cos(alpha_inf*numerix.pi/180.0) * numerix.cos(beta_inf*numerix.pi/180.0)
# V_inf = unorm_inf * numerix.sin(beta_inf*numerix.pi/180.0)
# roE_inf = press_inf / gm1 + 0.5 * ro_inf * unorm_inf**2
# Rsq = (X - X0)**2 + (Y-Y0)**2
# expon = numerix.exp(-0.5*Rsq/Rc**2)
# U_init = U_inf - (Y-Y0) / Rc * Gamma * expon
# V_init = V_inf + (X-X0) / Rc * Gamma * expon
# temp_init = temp_inf - 0.5 * (Gamma * Gamma) * numerix.exp(-Rsq / Rc**2) / Cp
# ro_init = ro_inf * (temp_init / temp_inf)**(1.0/gm1)
# pres_init = ro_init * Rgas * temp_init
# roE_init = pres_init / gm1 + 0.5*(U_init**2 + V_init**2)*ro_init
# ro.setValue( ro_init*(np.sqrt(Rsq)<=Rc*cut)+ro_inf*(np.sqrt(Rsq)>Rc*cut) )
# roU.setValue( ro_init*U_init*(np.sqrt(Rsq)<=Rc*cut)+ro_inf*U_inf*(np.sqrt(Rsq)>Rc*cut) )
# roV.setValue( ro_init*V_init*(np.sqrt(Rsq)<=Rc*cut)+ro_inf*V_inf*(np.sqrt(Rsq)>Rc*cut) )
# roE.setValue( roE_init*(np.sqrt(Rsq)<=Rc*cut)+roE_inf*(np.sqrt(Rsq)>Rc*cut) )
# p.setValue( pres_init*(np.sqrt(Rsq)<=Rc*cut)+press_inf*(np.sqrt(Rsq)>Rc*cut) )
ro_init = fipy.CellVariable(mesh=mesh, value=ro_init)
if parallel.procID == 0:
vi2D = fipy.Viewer(vars=p) #(ro,roU,roV,roE,p))
vi2D.plot()
# raw_input("Initialization. Press <return> to continue ...")
#==============================
umax = max(max(U_init), max(V_init))
CFL = 0.75
dx = numerix.sqrt(min(mesh.cellVolumes))
dt = dx * CFL / umax
if parallel.procID == 0:
print "Timestep for CFL=%.2f, t=%.8fs"%(CFL, dt)
#==============================
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.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=ro) == 0.0
#----
eqnMx = fipy.TransientTerm(var=roU) + \
fipy.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=roU) == \
- p.grad.dot([1., 0.])
# gm1*fipy.ExponentialConvectionTerm(var=(roE-0.5*(roU**2+rov**2)/ro)) == 0.0
# (1.0-0.5*gm1)*fipy.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=roU) + \
# gm1*roE.grad.dot([1., 0.]) + \
# 0.5*gm1*fipy.ExponentialConvectionTerm(coeff=coeffConvMinus.faceValue, var=roV) == 0.0
#----
eqnMy = fipy.TransientTerm(var=roV) + \
fipy.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=roV) == \
- p.grad.dot([0., 1.])
# (1.0-0.5*gm1)*fipy.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=roV) + \
# gm1*roE.grad.dot([0., 1.]) - \
# 0.5*gm1*fipy.ExponentialConvectionTerm(coeff=coeffConvMinus.faceValue, var=roU) == 0.0
#----
eqnE = fipy.TransientTerm(var=roE) + \
fipy.ExponentialConvectionTerm(coeff=coeffConv.faceValue, var=roE) == \
-(roU*p/ro).grad.dot([1.,0.]) - (roV*p/ro).grad.dot([0.,1.])
# -fipy.CentralDifferenceConvectionTerm(coeff=coeffConv.faceValue, var=p)
# (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
#----
# eqnP = fipy.TransientTerm(var=p) == (1.0/dt) * ( gm1*roE - 0.5*gm1*(roU**2 + roV**2)/ro - gm1*roE.old + 0.5*gm1*(roU.old**2 + roV.old**2)/ro.old )
eqnP = fipy.ImplicitSourceTerm(coeff=1.0, var=p) == fipy.ImplicitSourceTerm(coeff=gm1, var=roE) - 0.5*gm1*(roU**2 + roV**2)/ro
# gm1*(roE - roE.old) / dt + 0.5*gm1*( (roU-roU.old)**2/dt**2 + (roV-roV.old)**2/dt**2 )/(ro)
#----
eqn = eqnC & eqnMx & eqnMy & eqnE & eqnP
#==============================
mesh1D = fipy.Grid1D(dx=dx, nx=nx)
xx = mesh1D.cellCenters[0]
imin1 = np.argmin(np.fabs(Y - Ly/2.0))
ro_cut = fipy.CellVariable(name="ro_cut Y=Ly/2", mesh=mesh1D)
ro_cut.setValue( ro((mesh1D.cellCenters[0], Y[imin1]*numerix.ones(ro_cut.mesh.numberOfCells)), order=0) )
ro_cut_init = fipy.CellVariable(name="ro_cut_init Y=Ly/2", mesh=mesh1D)
ro_cut_init.setValue( ro_init((mesh1D.cellCenters[0], Y[imin1]*numerix.ones(ro_cut_init.mesh.numberOfCells)), order=0) )
#ro_inf*((temp_inf - 0.5 * (Gamma * Gamma) * numerix.exp(-( (xx-X0)**2 - (Y[imin1]-Y0)**2 ) / Rc**2) / Cp) / temp_inf)**(1.0/gm1)
imin2 = np.argmin(np.fabs(Y - 2.0*Ly/3.0))
ro_cut_higher = fipy.CellVariable(name="ro_cut Y=2Ly/3", mesh=mesh1D)
ro_cut_higher.setValue( ro((mesh1D.cellCenters[0], Y[imin2]*numerix.ones(ro_cut_higher.mesh.numberOfCells)), order=0) )
ro_cut_higher_init = fipy.CellVariable(name="ro_cut_init Y=2Ly/3", mesh=mesh1D)
ro_cut_higher_init.setValue( ro_init((mesh1D.cellCenters[0], Y[imin2]*numerix.ones(ro_cut_higher_init.mesh.numberOfCells)), order=0) )
if parallel.procID == 0:
vi1D = fipy.Viewer(vars=(ro_cut_init,ro_cut))
# vars=(np.fabs(ro_cut-ro_cut_init)/np.fabs(ro_cut_init), np.fabs(ro_cut_higher-ro_cut_higher_init)/np.fabs(ro_cut_higher_init)), ylog=True
vi1D.plot()
# raw_input("1D cut, before computation. Press <return> to continue ...")
# vi2D = fipy.Viewer(vars=ro)
tcur = 0.0
steps = 100
for step in range(steps):
ro.updateOld()
roU.updateOld()
roV.updateOld()
roE.updateOld()
# p.updateOld()
#
# eqn.solve(dt=dt)
eqnC.sweep(dt=dt)#, underRelaxation=0.2)
eqnMx.sweep(dt=dt)#, underRelaxation=0.2)
eqnMy.sweep(dt=dt)#, underRelaxation=0.2)
eqnE.sweep(dt=dt)#, underRelaxation=0.2)
#
p.setValue( gm1*roE - 0.5*gm1*((roU/ro)**2 + (roV/ro)**2) )
#
if parallel.procID == 0:
vi2D.plot()
ro_cut.setValue( ro((mesh1D.cellCenters[0], Y[imin1]*numerix.ones(ro_cut.mesh.numberOfCells)), order=0) )
if parallel.procID == 0:
vi1D.plot()
#
tcur += dt
#
if parallel.procID == 0:
print "Time t=%.2e, distance from exact 1D=%.6e, distance from exact 2D=%.6e"%(tcur, np.linalg.norm(ro_cut-ro_cut_init, ord=2)/np.linalg.norm(ro_cut_init, ord=2), np.linalg.norm(ro-ro_init, ord=2)/np.linalg.norm(ro_init, ord=2))
raw_input("Final state ...")
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
