Hi,
I'm trying to use fipy to solve some differential equations, and I'm
not sure how to work with the particular equation I've got at hand. I'm
a total newbie at this, so I'm probably doing something silly, but I
can't figure out how to proceed so I thought I'd ask for some help.
The attached Python file demonstrates my problem and what I've done so
far. The key issue is that I have an equation with diffusion terms that
have spatially-dependent prefactors; in the simplest case, something
like
L^2 d/dL (L df/dL)
where L is a spatial variable and I'm trying to determine f(L). As far
as I can tell, this kind of term is not directly allowed by fipy, and I
get a TermMultiplyError if I try to construct it.
In the simplest case, demonstrated in the attachment, I can divide out
the prefactor and fipy works great. However, in my full problem, there
will be multiple diffusion terms (in multiple dimensions) with
differing prefactors, so I won't be able to make them all disappear.
I'm pretty sure that I can generically re-express any prefactored
diffusion term as the sum of a prefactor-free diffusion term and an
advection term. However, I've tried to implement this with fipy, and it
doesn't work. The result from the solver is as if I didn't include the
advection term at all. Again this is demonstrated in the attached file.
So, I have two questions:
(1) I'm probably doing something wrong with the advection term in the
attached sample. Can anyone spot the problem?
(2) It'd be nice not to have to do the error-prone elimination of my
diffusion prefactors. Are there any ways to avoid this?
Thanks for any tips,
Peter#! /usr/bin/env python
# -*- mode: python; coding: utf-8 -*-
# Copyright 2015 Peter Williams. Licensed under the MIT License.
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
import fipy
""" ============================================================
We seek to find f(L) satisfying the differential equation
0 = L^2 d/dL[ L df/dL ] - x0^2 f / 4
(This comes from modeling of Van Allen belt particles; see Goertz et al
(1979), DOI 10.1029/JA084iA01p00087.)
For f(1) = 0 and f(4) = 1, the analytic solution is
f(L) = (K0(x0) I0(x) - I0(x0) K0(x)) / (K0(x0) I0(x2) - I0(0) K0(x2))
where I0, K0 are the modified Bessel functions of first/second kind, x = x0 /
sqrt(L), and x2 = x0 / 2. We implement the analytic solution here:
"""
l_range = (1., 4.)
l_count = 32
x0 = 15.
def analytic ():
from scipy.special import iv # modified Bessel function of first kind
from scipy.special import kv # modified Bessel function of second kind
offset = (l_range[1] - l_range[0]) / (2 * l_count)
l = np.linspace (l_range[0] + offset, l_range[1] - offset, l_count)
x = x0 * l**-0.5
kv0 = kv (0, x0)
iv0 = iv (0, x0)
kv2 = kv (0, x0 / 2)
iv2 = iv (0, x0 / 2)
f = (kv0 * iv (0, x) - iv0 * kv (0, x)) / (kv0 * iv2 - iv0 * kv2)
return f
print ('analytic:', analytic ())
""" ============================================================
fipy doesn't let us put variable coefficients in front of diffusion terms.
In this simple example, we can just divide out the L^2 and solve:
0 = d/dL[ L df/dL ] - x0^2 f / (4 L^2)
This works.
"""
def only_possible_for_simple_cases ():
mesh = fipy.Grid1D (
dx=(l_range[1] - l_range[0]) / l_count,
nx=l_count,
)
l = mesh.x + l_range[0]
f = fipy.CellVariable (name='f', mesh=mesh, value=0.)
f.constrain (0., mesh.facesLeft) # f=0 at L=1
f.constrain (1., mesh.facesRight) # f=1 at L=4
diffusion = fipy.DiffusionTerm (coeff=l)
losses = fipy.ImplicitSourceTerm (coeff=-x0**2 / (4 * l**2))
eq = (diffusion + losses)
eq.solve (var=f)
return np.array (f.value)
print ('approach 1:', only_possible_for_simple_cases ())
"""============================================================
For more complicated problems, however, we have multiple diffusion terms with
differing prefactors, so they can't all be divided out. I think it should
generally be possible to convert them to linear combinations of prefactor-free
diffusion terms and advection terms, though. In this example,
L^2 d/dL[ L df/dL ] = d/dL[ L^3 df/dL ] - 2 L^2 df/dL
Unfortunately, my best-effort attempt to implement this does not succeed. The
advection term seemingly has no effect, as shown below. I get identical
numbers from the solver regardless of whether I include it in the equation or
not.
"""
def im_doing_something_wrong (use_advection=True):
mesh = fipy.Grid1D (
dx=(l_range[1] - l_range[0]) / l_count,
nx=l_count,
)
l = mesh.x + l_range[0]
f = fipy.CellVariable (name='f', mesh=mesh, value=0.)
f.constrain (0., mesh.facesLeft) # f=0 at L=1
f.constrain (1., mesh.facesRight) # f=1 at L=4
diffusion = fipy.DiffusionTerm (coeff=l**3)
advection = fipy.AdvectionTerm (coeff=-2 * l**2)
losses = fipy.ImplicitSourceTerm (coeff=-x0**2 / 4)
if use_advection:
eq = (diffusion + advection + losses)
else:
eq = (diffusion + losses)
eq.solve (var=f)
return np.array (f.value)
print ('approach 2a:', im_doing_something_wrong (use_advection=True))
print ('approach 2b:', im_doing_something_wrong (use_advection=False))
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
