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.
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.
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?
Should the ways described be equivalent? If though, why do they give
different solutions?
May my problems in the second case depend on solver type? I noticed some
dependency ...
Thanx to all who has any ideas!
I attach the script I use:
from fipy.tools.numerix import exp
from fipy.tools.numerix import sqrt
from fipy.tools.numerix import log10
L = 1000.
nx = 1000
dx = L / nx
from fipy.meshes.grid1D import Grid1D
mesh = Grid1D(dx = dx, nx = nx)
axis = 0
x = mesh.getCellCenters()[:,axis]
Xmin = 0.01
Xmax = L
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 Lder(E):
Loss_derE = (1e-12)*(0.65e-3)*B*B*(1/(Ee0*Ee0))*2*E+ \
(1e-12)*(2.66e-14)*Ucmb*(1/(Ee0*Ee0))*2*E+ \
(1.37e-16)*n*(log10(E/Ee0+0.36)+E/(E+0.36*Ee0))
return Loss_derE*yr*1/1.
def Jin(x):
E=x
return Ke*A*E**alpha*exp(-E/Ecut)
# return 50*exp(-(x-50)**2/25)
from fipy.variables.faceVariable import FaceVariable
from fipy.variables.cellVariable import CellVariable
from fipy.variables.vectorFaceVariable import VectorFaceVariable
from fipy.variables.vectorCellVariable import VectorCellVariable
X = mesh.getFaceCenters()[...,0]
DiffCoeff = FaceVariable(mesh=mesh)
#DiffCoeff.setValue(0)
DiffCoeff = 0
ConvCoeff = VectorCellVariable(mesh=mesh)
ConvCoeff.setValue(L(X))
#ConvCoeff = 1
AdvCoeff = CellVariable(mesh=mesh)
AdvCoeff.setValue(L(x))
#initial conditions
phi = CellVariable(mesh=mesh, name = "variable")
phi.setValue(Jin(x))
phiInitial = CellVariable(name="variable initial", mesh=mesh)
phiInitial.setValue(Jin(x))
Source = CellVariable(mesh=mesh, name = "source")
Source.setValue(phi*Lder(x))
#Source = 0
#boundary conditions
#valueLeft = Jin(Xmin)
#valueRight = Jin(Xmax)
valueLeft = 0
valueRight = 0
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.boundaryConditions.fixedValue import FixedValue
boundaryConditions = (FixedFlux(faces=mesh.getFacesLeft(),
value=valueLeft), FixedFlux(faces=mesh.getFacesRight(), value=valueRight))
#boundaryConditions = (FixedValue(faces=mesh.getFacesLeft(),
value=valueLeft), FixedValue(faces=mesh.getFacesRight(), value=valueRight))
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.powerLawConvectionTerm import PowerLawConvectionTerm
from fipy.models.levelSet.advection.advectionEquation import
buildAdvectionEquation
advEqn = buildAdvectionEquation(advectionCoeff=AdvCoeff)
diffTermI = ImplicitDiffusionTerm(coeff=DiffCoeff)
eq = advEqn + Source
eqI = TransientTerm() == diffTermI +
PowerLawConvectionTerm(coeff=ConvCoeff, diffusionTerm=diffTermI)
timeStepDuration = 0.01
steps = 1000
from fipy.viewers import make
viewer = make(vars=(phi,phiInitial),limits={'datamin': -1e-15,
'datamax': 1e3})
viewer.xlog=1
viewer.ylog=1
from fipy.viewers.tsvViewer import TSVViewer
#TSVViewer(vars=(phi)).plot(filename="res.tsv")
from fipy.solvers.linearLUSolver import LinearLUSolver
for step in range(steps):
# eqI.solve(var=phi, boundaryConditions=boundaryConditions,
dt=timeStepDuration, solver = LinearLUSolver(tolerance = 1.e-15))
eq.solve(var=phi, boundaryConditions=boundaryConditions,
dt=timeStepDuration, solver = LinearLUSolver(tolerance = 1.e-15))
viewer.plot()
raw_input("Press enter to close")
