On Thu, Jul 31, 2008 at 10:57 AM, Igor Telezhinsky <[EMAIL PROTECTED]> wrote:
>
> Dear all!
>
> Recently I found fipy is a powerful tool for solving pdes. I have an
> equation to solve of the following type:
>
> \frac{\partial \phi(E,t)}{\partial t} = \frac{\partial}{\partial E}
> \left(L(E) * \phi(E,t) \right) (***)
>
> Apparently, I guess there should be two ways of solving this equation in
> fipy
>
> 1. Using the ConvectionTerm with the convection coefficient (ConvCoeff)
> equal to function L(E), where the DiffusionTerm used to generate the
> ConvectionTerm equals zero. This I do by equating diffusion coefficient
> (DiffCoeff) to zero.
You don't need to specify the diffusion term You just need to write,
eqn = TransientTerm() == PowerLawConvectionTerm(L)
> eqI = TransientTerm() == diffTermI + PowerLawConvectionTerm(coeff=ConvCoeff,
> diffusionTerm=diffTermI)
>
> This way I'm not sure as ConvCoeff must be vector and function L(E) is a
> scalar.
L needs to be a vector to work correctly with the convection term. If
you define E to be mesh.getCellCenters() then L will automatically be
a vector.
> 2. Using the AdvectionEquation, though some analytical operations must be
> done:
> a) rewrite the right-hand side of the equation (***) as
>
> L(E) * \frac{\partial}{\partial E}\phi(E,t) + \phi(E,t) *
> \frac{\partial}{\partial E}L(E)
>
> and take the derivative of L(E) analytically naming it as Lder(E) function;
> b) introduce the Source that would be equal to \phi(E,t)*Lder(E):
>
> eq = advEqn + Source
>
> The difficulties I met are below.
>
> In the first way:
>
> a) the solution is insensitive to the dependence of ConvCoeff on a variable.
> b) the solution doesn't seem physical to me concerning the problem I'm
> solving (electrons cooling and their spectrum deformation).
>
> In the second way:
>
> a) the solution first seems physical, although with time very large negative
> numbers appear.
> b) the solution is insensitive to the boundary conditions
>
> Questions:
> What is the proper way to introduce equation (***) to the fipy?
The first way.
> Should the ways described be equivalent? If though, why do they give
> different solutions?
The AdvectionEquation should only really be used with a level set
variable. Don't use it for this problem.
> May my problems in the second case depend on solver type? I noticed some
> dependency ...
Advection problems require non-symmetric solvers, so yes there is some
dependence.
I updated your script to work with revision 2713. Things are a lot
cleaner. I'm assuming that you are using fipy version 1.2. I will test
your script against 1.2, but before I do that, I would like to know
what the intended boundary conditions are.
from fipy import *
Lx = 1000.
nx = 1000
dx = Lx / nx
mesh = Grid1D(dx = dx, nx = nx)
x, = mesh.getCellCenters()
Xmin = 0.01
Xmax = Lx
Ee0=5.11e-7
n=122.5
B=0.0001
Ucmb=0.25
A=5.075
Ke=1
alpha=-2
Ecut = 10
yr=3.156E+7
def L(E):
Loss = (1e-12)*(0.65e-3)*B*B*(1/(Ee0*Ee0))*E*E+ \
(1e-12)*(2.66e-14)*Ucmb*(1/(Ee0*Ee0))*E*E+ \
(1.370e-16)*n*E*log10(E/Ee0+0.36)+ \
(8.978e-13)*n*Ee0 #*(E/sqrt(E*E-Ee0*Ee0))
return Loss*yr*1/1.
def Jin(x):
E=x
return Ke*A*E**alpha*exp(-E/Ecut)
##X = mesh.getFaceCenters()[...,0]
convCoeff = FaceVariable(mesh=mesh, value=L(mesh.getFaceCenters()), rank=1)
##convCoeff[...,0] = convCoeff[...,-1] = 0
print convCoeff[...,0]
print convCoeff[...,-1]
raw_input()
#initial conditions
phi = CellVariable(mesh=mesh, name = "variable", value=Jin(x))
phiInitial = CellVariable(name="variable initial", mesh=mesh)
phiInitial.setValue(Jin(x))
valueLeft = 0
valueRight = 0
boundaryConditions = (FixedFlux(faces=mesh.getFacesLeft(), value=valueLeft),
FixedFlux(faces=mesh.getFacesRight(), value=valueRight))
eq = TransientTerm() == PowerLawConvectionTerm(convCoeff)
timeStepDuration = 1.0
steps = 1000
viewer = make(vars=(phi,phiInitial),limits={'datamin': -1e-15, 'datamax':
1e10})
viewer.xlog=1
viewer.ylog=1
for step in range(steps):
eq.solve(var=phi, boundaryConditions=boundaryConditions,
dt=timeStepDuration)
viewer.plot()
print min(phi)
raw_input("Press enter to close")
--
Daniel Wheeler
