Daniel,
vel is a 3 by number of faces array... x,y,z components of velocity at
each face. In the array I'm using, vel[0] and vel[1] are all nearly zero,
and vel[2] is + or - 2.5e-4. The magnitude of these numbers is
correct.
passing 'vel' to the conductance term doesn't work. I get an error:
<type 'exceptions.TypeError'>: The coefficient must be a vector value.
WARNING: Failure executing file: <PeatCylSimple.py>
Both of the other options don't result in an error, just a wrong answer.
Should the coeff in the convection term be a vector normal to the face?
As an aside, the equation I'm solving is the standard Advection-diffusion
equation:
dC/dx=d/dX(D dC/dX) - d(v*C)/dX
v=velocity
C=concentration
D=diffusion coef
I get a sensible answer only when I switch the sign on the diffusion term
in my transport equation.
Here's the script I'm using (running anisotropy Branch of FiPy). The last
couple of lines calculate the analytic solution.:
from scipy.special import erfc
def addLayer(verts):
addArray = verts.copy()
addArray[:] = 0
addArray[2] = length/slices
return verts + addArray
###create mesh, hemi-cylinder of peat, 40 cm long, 2.54 cm radius
length=40. ##cm
slices=40
mesh=GmshImporter2D('./HalfCirc.msh').extrude(layers=slices,
extrudeFunc=addLayer)
###calculate hydraulic head in peat
head=CellVariable(name='hydraulic head',
mesh=mesh,value=100.,hasOld=1) ## in cm
Dist=random.uniform(1.,1.,(mesh.getNumberOfCells()))
Kxx=1.e-6*100*Dist ##hyd. cond. in cm/sec
Ss=1.e-2/100. ##Sp.Storage in 1/cm
kond=CellVariable(name='hydraulic conductivity',mesh=mesh,rank=0,value=Kxx)
spStore=CellVariable(name='spec. storage',mesh=mesh,rank=0,value=Ss)
effPor=CellVariable(name='effective porosity',mesh=mesh,rank=0,value=0.2)
gwFlowEq=TransientTerm(spStore)==DiffusionTerm(kond.getHarmonicFaceValue())
##hydraulic gradients across peat are 0.5 or 0.8, units in cm, 40 cm long tube
X,Y,Z=mesh.getFaceCenters()
ZMax=max(Z)-0.05
bndry=(FixedValue(mesh.getFaces().where((Z<0.05)),120.),
FixedValue(mesh.getFaces().where((Z>ZMax)),100.))
delTime=1.e12
timeSteps=3
for step in range(timeSteps):
head.updateOld()
resid=1.
while resid>1.e-12:
resid=gwFlowEq.sweep(var=head,
solver=LinearPCGSolver(iterations=500,
tolerance=1.e-13),
boundaryConditions=bndry,
dt=delTime)
print 'hyd.head',resid
discharge=dot(head.getFaceGrad(),kond.getHarmonicFaceValue())
vel=discharge/effPor.getArithmeticFaceValue()
conc=CellVariable(name='solute conc',mesh=mesh,
value=1.,hasOld=1)
D=.0002
Disp=FaceVariable(name='dispersion tensor',
mesh=mesh,rank=0)
Disp.setValue(D)
X,Y,Z=mesh.getFaceCenters()
ZMax=max(Z)-0.05
bndry2=(FixedValue(mesh.getFaces().where((Z<0.05)),101.),
FixedValue(mesh.getFaces().where((Z>ZMax)),1.))
transportEq2=(TransientTerm()==
-PowerLawConvectionTerm(coeff=vel())
-DiffusionTerm(coeff=Disp))
t=0
delTime=1.e3
timeSteps=40
for step in range(timeSteps):
t+=delTime
conc.updateOld()
resid=1.
while resid>1.e-12:
#for i in range(10):
resid=transportEq2.sweep(var=conc,
solver=LinearGMRESSolver(iterations=500,
tolerance=1.e-15),
boundaryConditions=bndry2,
dt=delTime)
print 'conc',resid
X,Y,Z=mesh.getCellCenters()
v=2.5e-4
Dis=D
OB=[1+100*.5*(erfc((x-v*t)/(2*sqrt(Dis*t)))+exp(v*x/Dis)*erfc((x+v*t)/(2*sqrt(Dis*t))))
for x in Z]
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
On Jan 22, 2008 1:41 PM, asreeve <[EMAIL PROTECTED]> wrote:
Hello,
Hi, hope all is well.
I've having some difficulty figuring out how to use the convection terms
in a 3-D problem.
I haven't solved a lot of 3D convection problems so I am not sure how well
tested everything is. Should work though.
I've calculated a velocity term for a 3-D simulations:
vel=dot(head.getFaceGrad(),kond.getHarmonicFaceValue())
I'm confused about which axis is getting dotted. What are the shapes of "kond"
and "vel"?
and would like to pass these velocities to a convection term.
I assume that "vel" is rank 1. Is it?
I've simplified my problem to the point that there are non-zero velocities
only in the z direction (but I need to solve a problem where vel varies in
space and direction).
Shouldn't be a problem.
I've tried using :
PowerLawConvectionTerm(coeff=vel())
and
PowerLawConvectionTerm(coeff=(vel()[0],vel()[1],vel()[2]))
in the equation term, but both of these yield an incorrect result.
However, when I assign the following (values taken from vel variable),
PowerLawConvectionTerm(coeff=(0,0,2.5e-4))
the output is correct.
What's the proper syntax to assign the calculated 'vel' to a convection
term?
Have you tried?
PowerLawConvectionTerm(coeff=vel)
Cheers
--
Daniel Wheeler
