Dear developers,
thank you for the answers to my previous posts. Now, as I have new
fipy version installed, I want to ask you some questions related to
it.
I have a problem with multiple convection terms and diffusion term and
all of this in 2D (non-zero components of del operator are denoted by
x and y). Could you tell me if it is correct to break the problem into
subproblems in the following way:
"basic" problem:
eq = TransTerm == DiffTerm_x - ConvTerm_x + ConvTermA_y - ConvTermB_y + Source
..
while time <= timeStop:
eq.solve(var=phi,...)
"broken" problem:
eqX = TransTerm == DiffTerm_x - ConvTerm_x + Source
eqAY = TransTerm == ConvTermA_y
eqBY = TransTerm == -ConvTermB_y
..
while time <= timeStop:
eqX.solve(var=phi,...)
eqAY.solve(var=phi,...)
eqBY.solve(var=phi,...)
I tested it and the results agreed up to 4th digit after coma.
I thought it would be some way out of the situation because in fact my
ConvTermA_y is needed to be multiplied by some function f(param)
(independent on y) and in previous fipy version I couldn't do it
directly. The only way to multiply a term by some function so that del
operator would not work on it was the introduction of the coefficient
in transient term, so knowing the above examples agree very well, I
introduced multiplication of the needed term by "breaking" the problem
in the above way, but eqAY changed:
eqAY = TransientTerm(coeff=1./f(param)) == ConvTermA_y
New fipy lets me do multiplication directly in "basic" problem:
eq = TransTerm == DiffTerm_x - ConvTerm_x + ConvTermA_y*f(param) -
ConvTermB_y + Source
so I tried it and checked the results with "broken" problem results.
Surprisingly, they didn't agree! So, I questioned the way the problem
is broken!
HOWEVER, It turned out that changing f(param) does not change ANYTHING
in the solution of the "basic" problem, while in "broken" problem I
clearly see the difference when changing f(param) and this change
seems reasonable: drifting in y direction. MOREOVER, it seems that
simply introduction of f(param) (even if f(param)==1.0) to "basic"
problem changes the result, whereas in "broken" problem the result as
expected doesn't chenge.
Could you, please, help me understand this issue?
Below I give the code.
NB: I tested that in new fipy
eqAY = TransientTerm(coeff=1./f(param)) == ConvTermA_y
and
eqAY = TransientTerm() == ConvTermA_y*f(param)
are equivalent so I write it in the new way in the code I give you.
the bottom function mkslice make a slice at x=1 and writes down to
file to compare results.
Thanks in advance,
Igor.
Code:
-----------------------------------------------
from fipy import *
def V(r):
return r
def f(param):
return param
Xmax=1.0
Ymax=1.0
nsteps=100
dx = Xmax / float(nsteps)
dy = Ymax / float(nsteps)
mesh_xy = Grid2D(dx = dx, dy = dy, nx = nsteps, ny = nsteps) +
((dx/2.,),(-Ymax/4.,))
xc,yc = mesh_xy.getCellCenters()
xf,yf = mesh_xy.getFaceCenters()
valueLeft = 0
valueRight = 0
BCs = (FixedValue(faces=mesh_xy.getFacesLeft(), value=valueLeft),
FixedValue(faces=mesh_xy.getFacesRight(), value=valueRight))
phi = CellVariable(mesh=mesh_xy, name = "N", value=0.0)
Source = CellVariable(mesh=mesh_xy, value = 1.0)
diffCoeff_x = FaceVariable(mesh=mesh_xy, value=0.0, rank=2)
diffCoeff_x[0,0] = 1.0
diffTerm_x = ImplicitDiffusionTerm(diffCoeff_x)
convCoeff_x = FaceVariable(mesh=mesh_xy, value=0.0, rank=1)
convCoeff_x[0] = V(xf)
convTerm_x = PowerLawConvectionTerm(convCoeff_x)
convCoeffA_y = FaceVariable(mesh=mesh_xy, value=0.0, rank=1)
convCoeffA_y[1] = V(yf)
convTermA_y = PowerLawConvectionTerm(convCoeffA_y)
convCoeffB_y = FaceVariable(mesh=mesh_xy, value=0.0, rank=1)
convCoeffB_y[1] = V(yf)
convTermB_y = PowerLawConvectionTerm(convCoeffB_y)
param=1.0
#eq = TransientTerm() == diffTerm_x - convTerm_x + convTermA_y -
convTermB_y + Source
eq = TransientTerm() == diffTerm_x - convTerm_x + convTermA_y*f(param)
- convTermB_y + Source
eqX = TransientTerm() == diffTerm_x - convTerm_x + Source
#eqAY = TransientTerm() == convTermA_y
eqAY = TransientTerm() == convTermA_y*f(param)
eqBY = TransientTerm() == -convTermB_y
timeStepDuration = (10.0*dx**2)/(2*1.0)
time = 0
timeStop = 100*timeStepDuration
viewer = Viewer(vars=(phi), limits={'datamin': 1e-2, 'datamax': 0.5})
while time <= timeStop:
timeStep = timeStepDuration
# eq.solve(var=phi, boundaryConditions=BCs, dt=timeStep)
eqX.solve(var=phi, boundaryConditions=BCs, dt=timeStep)
eqAY.solve(var=phi, boundaryConditions=BCs, dt=timeStep)
eqBY.solve(var=phi, boundaryConditions=BCs, dt=timeStep)
time=time+timeStep
viewer.plot()
def mkslice(N=phi,xc=xc,nrsteps=nsteps,dx=dx,yc=yc,npsteps=nsteps,dy=dy):
LOCATIONinX=1.0
dimX=nrsteps
dimY=npsteps
j=0
f=open("edge.txt","w")
for step in range(dimX*dimY):
if (xc.item(j) >= LOCATIONinX-dx/3. and xc.item(j) <=
LOCATIONinX+dx/3.):
X=yc.item(j)
f.write(str(X)+" "+str(N[j])+"\n")
j=j+1
f.close()
mkslice()
