1D coupled fluid equations in FiPy

[Dan Brunner]

All,

First, thanks for developing and supporting FiPy. I am just starting to
learn to use it with the hopes that it will become my standard tool for
solving PDEs.

I've hit a bottle-neck in setting up my problem and cannot piece it all
together from the examples and FAQs.

I'm working on solving (building up to more complex, two-fluid plasma
equations):
[image: \frac{\partial}{\partial t}n + \frac{\partial}{\partial
x}\left(nv\right)=S_n]

[image: \frac{\partial}{\partial
t}\left(mnv\right)+\frac{\partial}{\partial x}\left(mnv^2+P\right)=S_m]

I've tried splitting up the derivatives, as suggested in the FAQs, and
including them as source terms. But cannot figure out how to properly do
this. E.g., including a variable that I am solving for as the coefficient
to a convection term. I can provide clarification as needed.

Thanks
~Dan Brunner

[Daniel Wheeler]

Post by Dan Brunner
All,
First, thanks for developing and supporting FiPy. I am just starting to
learn to use it with the hopes that it will become my standard tool for
solving PDEs.

I hope it works out.

Post by Dan Brunner
I've hit a bottle-neck in setting up my problem and cannot piece it all
together from the examples and FAQs.

Have you tried solving for "mnv" as the independent variable and
making "m" simply "mnv / (n * v)" assuming you need "m" for one of the
P or S variables? That might work. You could also then try doing that
as a fully coupled equation.

These equiation have no diffusion terms so they will be hard to solve.
Depending on what you're looking for in the solution, you really need
Riemann coupling and high order advection terms, which FiPy doesn't
currently have. Actually, there is no explicit coupling of the
equations within the derivatives unless "P" has "n" dependence so
Riemann does nothing.

I probably need to see how your implementing things to be more helpful.

[Dan Brunner]

Thanks for you help so far. I will be more explicit in what I am doing. The
equations are:

[image: \frac{\partial}{\partial t}\left(n\right)+\frac{\partial}{\partial
x}\left(nv\right)=S_n]
[image: \frac{\partial}{\partial t}\left(m_i n
v\right)+\frac{\partial}{\partial x}\left(m_i
nv^2\right)=S_m-\frac{\partial}{\partial x}\left(nT_e+nT_i\right)]
[image: \frac{\partial}{\partial
t}\left(\frac{3}{2}nT_e\right)+\frac{\partial}{\partial
x}\left(\frac{5}{2}nvT_e\right)=S_{q,e}]
[image: \frac{\partial}{\partial t}\left(\frac{3}{2}nT_i+\frac{1}{2}m_i
nv^2\right)+\frac{\partial}{\partial
x}\left(\frac{5}{2}nvT_i+\frac{1}{2}m_i nv^3\right)=S_{q,i}]

With n, v, Te, and Ti the variables, x the spatial dimension, t the
temporal dimension, and m a constant. The S-terms are source terms, which
will get more complicated as I progress (including a thermal conductivity
with a T-dependent coefficient). To solve, for now ignoring the terms with
m in the last equation for simplicity, I have tried recasting them in the
form:

[image: \frac{\partial}{\partial
t}\left(f_0\right)+\frac{\partial}{\partial
x}\left(\frac{1}{m}f_1\right)=S_0]

[image: \frac{\partial}{\partial
t}\left(f_1\right)+\frac{\partial}{\partial
x}\left(\frac{1}{m}\frac{f_1^2}{f_0}\right)=S_1-\frac{\partial}{\partial
x}\left(\frac{2}{3}f_2+\frac{2}{3}f_3\right)]

[image: \frac{\partial}{\partial
t}\left(f_2\right)+\frac{\partial}{\partial x}\left(\frac{3}{5m}\frac{f_1
f_2}{f_0}\right)=S_{2}]

[image: \frac{\partial}{\partial
t}\left(f_3\right)+\frac{\partial}{\partial x}\left(\frac{3}{5m}\frac{f_1
f_3}{f_0}\right)=S_{3}]

I roughly follow the process laid out in the coupled diffusion equations
example:

from fipy import *

mesh = Grid1D(nx = 100, Lx = 1.0)

m = 1.0
f0 = CellVariable(mesh = mesh, hasOld = True, value = 0.5, name = 'f0')

f1 = CellVariable(mesh = mesh, hasOld = True, value = 0.5, name = 'f1')
f2 = CellVariable(mesh = mesh, hasOld = True, value = 0.5, name = 'f2')
f3 = CellVariable(mesh = mesh, hasOld = True, value = 0.5, name = 'f3')

f0.constrain(0, mesh.facesLeft)
f0.constrain(1, mesh.facesRight)
f1.constrain(0, mesh.facesLeft)
f1.constrain(1, mesh.facesRight)
f2.constrain(0, mesh.facesLeft)
f2.constrain(1, mesh.facesRight)
f3.constrain(0, mesh.facesLeft)
f3.constrain(1, mesh.facesRight)
s0 = 0.0
s1 = 0.0
s2 = 0.0
s3 = 0.0

eq0 = TransientTerm(var = f0) + ExponentialConvectionTerm(coeff = (1.0/m,),
var = f1) == s0
eq1 = TransientTerm(var = f1) + ExponentialConvectionTerm(coeff = (1.0/m,),
var = f1**2/f0) == s1 - ExponentialConvectionTerm(coeff = (2/3.0,), var =
f2) - ExponentialConvectionTerm(coeff = (2/3.0,), var = f3)
eq2 = TransientTerm(var = f2) + ExponentialConvectionTerm(coeff =
(3.0/(5*m),), var = f1*f2/f0) == s2
eq3 = TransientTerm(var = f3) + ExponentialConvectionTerm(coeff =
(3.0/(5*m),), var = f1*f3/f0) == s3
eqs = eq0 & eq1 & eq2 &eq3
for t in range(1):

f0.updateOld()

f1.updateOld()

f2.updateOld()

f3.updateOld()
eqs.solve(dt = 1.0e-3)

viewer = Viewer((f0, f1, f2, f3))

Which results in the error:

fipy.terms.SolutionVariableNumberError: Different number of solution
variables and equations.

I believe that I am ignorant to something fundamental about FiPy and would
appreciate any guidance as to how to cast these equations in a form that
FiPy will take in.

~Dan

Post by Daniel Wheeler

Post by Dan Brunner
All,
First, thanks for developing and supporting FiPy. I am just starting to
learn to use it with the hopes that it will become my standard tool for
solving PDEs.

I hope it works out.

Post by Dan Brunner
I've hit a bottle-neck in setting up my problem and cannot piece it all
together from the examples and FAQs.

Have you tried solving for "mnv" as the independent variable and
making "m" simply "mnv / (n * v)" assuming you need "m" for one of the
P or S variables? That might work. You could also then try doing that
as a fully coupled equation.
These equiation have no diffusion terms so they will be hard to solve.
Depending on what you're looking for in the solution, you really need
Riemann coupling and high order advection terms, which FiPy doesn't
currently have. Actually, there is no explicit coupling of the
equations within the derivatives unless "P" has "n" dependence so
Riemann does nothing.
I probably need to see how your implementing things to be more helpful.
--
Daniel Wheeler
_______________________________________________
fipy mailing list
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

[Daniel Wheeler]

I think FiPy is really badly suited for coupled hyperbolic equations.
It is essential to have a "Riemann-type solver" for this type of
problem to correctly get the cross-coupling in the waves across the
cell faces. I would really recommend using CLAWPACK for this. FiPy is
much better suited for convection-diffusion problems and not so good
for pure hyperbolic equations. It would be great to have that
functionality in FiPy, but unfortunately it isn't available.
--
Daniel Wheeler