Hi Tony,
This is an interesting issue. Currently in fipy, when a fourth order
equation is specified there is a limitation in the combination of
boundary conditions that can be applied to each face. Certainly, the
problem you are trying to solve breaks this paradigm.
For a fourth order equation, the limitation is that any given face
must have two boundary conditions specified. Of the two boundary
condition on each face, one must be either zero or first order and the
other must be second or third order. In the example you sent there is
a zero and first order boundary condition on the left and a second and
third order boundary condition on the right. This doesn't fit the
paradigm.
I believe this problem is related to the underlying discretization and
whether the problem is categorized as a boundary value problem or an
initial value problem. Certainly, with a second order diffusion
equation it is not possible to specify both zero and first order on
the same face. This transforms the problem to an IVP. Fourth order
equations are less clear. The way you have specified the boundary
condition seems like a BVP at first glance, but at least with the
discretization that we use this specification does not work. This is a
thorny issue so I haven't really documented it. FiPy should at least
raise an informative error. I will submit a ticket for that at least.
Below is an example of a way to specify the boundary conditions that
does "work". There is a zero and second order boundary condition on
the left and a first and third order boundary condition on the right.
It matches the analytical result.
I hope you find this useful.
Cheers
from fipy import *
Lx = 1.
nx = 10000
dx = Lx / nx
mesh = Grid1D(dx=dx, nx=nx)
x = mesh.getCellCenters()[0]
k = 5.
f = 5.
y_analytical = CellVariable(mesh=mesh, name='ana',
value = f / k**3 * (-sin(k*x) / cos(k) + k*x))
# fixed at x = 0; forced at x = Lx
BCs = (FixedValue(mesh.getFacesLeft(), value=0.),
FixedFlux(mesh.getFacesRight(), value=0.),
NthOrderBoundaryCondition(mesh.getFacesLeft(), value=0.,
order=2),
NthOrderBoundaryCondition(mesh.getFacesRight(), value=f,
order=3))
# y_xxxx + k^2 * y_xx = 0
ElasticityTerm = ImplicitDiffusionTerm(coeff=(1.0, 1.0))
InertiaTerm = ImplicitDiffusionTerm(coeff=k**2)
eq = ElasticityTerm + InertiaTerm == 0.
solver = LinearLUSolver()
y = CellVariable(mesh=mesh, name='fipy', value=0.)
viewer = viewers.make(vars=(y, y_analytical))
eq.solve(var=y, boundaryConditions=BCs, solver=solver)
viewer.plot()
raw_input()
On Mon, May 26, 2008 at 5:03 PM, Tony S Yu <[EMAIL PROTECTED]> wrote:
Hi,
I wrote a simple test script for a 1D ODE. I'm having trouble
getting the
fipy solution to match the analytical solution. I'm pretty sure the
analytical solution is correct (but I guess I could be wrong). One
thing
that really bothers me is that changing the number of points (nx in
the
script below) drastically changes the solution fipy returns.
If you care about the physics of the problem: This a linearized
beam bending
problem with the position and slope fixed at one end (x = 0) and a
dimensionless force 'f' applied at x = L. I was planning to add a
time
dependent term in later, but I wanted to make sure I could match a
simpler
problem with its analytical solution.
Thanks in advance for your help.
-Tony
Specs:
OS X 10.5.2
FiPy trunk
Numpy 1.0.4
Python 2.5.1
#===========
from __future__ import division
from fipy import *
Lx = 1.
nx = 100
dx = Lx / nx
mesh = Grid1D(dx=dx, nx=nx)
x = mesh.getCellCenters()[0]
k = 5
f = 5
y_analytical = CellVariable(mesh=mesh, name='analytical',
value = f / k**3 * (sin(k)*(cos(k*x) - 1) + cos(k)*(-sin(k*x) +
k*x)))
# fixed at x = 0; forced at x = Lx
BCs = (FixedValue(mesh.getFacesLeft(), value=0.),
FixedFlux(mesh.getFacesLeft(), value=0.),
NthOrderBoundaryCondition(mesh.getFacesRight(), value=0.,
order=2),
NthOrderBoundaryCondition(mesh.getFacesRight(), value=f,
order=3))
# y_xxxx + k^2 * y_xx = 0
ElasticityTerm = ImplicitDiffusionTerm(coeff=(1.0, 1.0))
InertiaTerm = ImplicitDiffusionTerm(coeff=k**2)
eq = ElasticityTerm + InertiaTerm == 0.
solver = LinearLUSolver()
y = CellVariable(mesh=mesh, name='fipy', value=0.)
viewer = viewers.make(vars=(y, y_analytical))
eq.solve(var=y, boundaryConditions=BCs, solver=solver)
viewer.plot('./test')
--
Daniel Wheeler
