Daniel,
I've simplified script a little more and have attached it. I think the
answers to all you're questions are listed below. Let me know if I've
missed something.
I changing the sign on the diffusion term and it did not change the
results. I've tried changing the line:
gwFlowEq=TransientTerm(spStore)==DiffusionTerm(kond.getHarmonicFaceValue())
to
gwFlowEq=TransientTerm(spStore)==-DiffusionTerm(kond.getHarmonicFaceValue())
or
gwFlowEq=TransientTerm(spStore)==DiffusionTerm(-kond.getHarmonicFaceValue())
and the calculated values (head) remains the same. I agree that the way you
wrote it is correct.
vel.shape does return scalers but vel().shape doesn't:
In [2]: vel.shape
Out[2]: (161,)
In [3]: vel().shape
Out[3]: (3, 161)
To get the matrix results I sent before,
I do the following in ipython after running the attached script:
In [37]: A=trMtrx.getNumpyArray()
In [38]: A[:3,:3]
Out[38]:
array([[-0.0064516, 0.0016129, 0. ],
[ 0. , -0.0080645, 0.0016129],
[ 0. , 0. , -0.0080645]])
I need to write out the equations to see if the magnitudes and signs
make sense. They won't be equal in magnitude, as I indicated before
and I didn't think about the impact of the transient term on the
matrix, which is part of the main diagonal. This is the matrix fipy is
using to solve the PDE, right?
If I substitute +DiffusionTerm(coeff=Disp) for
-DiffusionTerm(coeff=Disp) in transportEq2, the solver doesn't converge.
The most sensible result occurs when I set the transport equation up
as:
transportEq2=(TransientTerm()==
-UpwindConvectionTerm(coeff=vel())
-DiffusionTerm(coeff=Disp)
+fluxIn
+ImplicitSourceTerm(coeff=fluxOut)
)
comparing this to the analytic solution results in curves that at
least have similar shapes. Flipping the sign on the diffusion equation
should be the right thing to do, but results in negative
concentrations. I need to check the boundary conditions again...
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
$MeshFormat
2 0 8
$EndMeshFormat
$Nodes
4
1 0 0 0
2 2.54 0 0
3 5.08 0 0
4 2.53999999999888 2.54 0
$EndNodes
$Elements
7
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
4 1 3 0 1 0 3 4
5 1 3 0 1 0 4 1
6 1 3 0 2 0 1 3
7 2 3 0 4 0 3 4 1
$EndElements
from fipy import *
import pylab
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 column of peat
head=CellVariable(name='hydraulic head',
mesh=mesh,value=100.,hasOld=1) ## in cm
Kxx=1.e-6*100*ones([mesh.getNumberOfCells()]) ##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=.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:
###don't need to use sweep now, but will need this in more complex
models
resid=gwFlowEq.sweep(var=head,
solver=LinearPCGSolver(iterations=500,
tolerance=1.e-13),
boundaryConditions=bndry,
dt=delTime)
print 'hyd.head',resid
print head
###calculate spec. discharge and velocity
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',value=D,mesh=mesh,rank=0)
## bndry2=(FixedFlux(mesh.getExteriorFaces(),00.))
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.))
X,Y,Z=mesh.getCellCenters()
###Do I need to explicity include influx term at boundary???
fluxIn=CellVariable(name='inflow',mesh=mesh,value=0.0)
ZValue=length/(slices*2.)
fluxIn.setValue(101.*discharge.getDivergence(),where=(Z==ZValue))
### Need term at outflow boundary to account for advective efflux
ZValue=length-length/(slices*2.)
fluxOut=CellVariable(name='outflow',mesh=mesh,value=0)
fluxOut.setValue(discharge.getDivergence(), where=(Z==ZValue))
transportEq2=(TransientTerm()==
-UpwindConvectionTerm(coeff=vel())
-DiffusionTerm(coeff=Disp)
+fluxIn
+ImplicitSourceTerm(coeff=fluxOut)
)
t=0
delTime=1.e3
timeSteps=10
for step in range(timeSteps):
t+=delTime ##need this for analytic solution at end of script
conc.updateOld()
resid=1.
transportEq2.solve(var=conc,
solver=LinearCGSSolver(iterations=500,
tolerance=1.e-15),
boundaryConditions=bndry2,
dt=delTime)
transportEq2.cacheMatrix() ##need this to look at matix built by fipy
transportEq2.cacheRHSvector()
trMtrx=transportEq2.getMatrix() ##returns matrix built by fipy
trRHS=transportEq2.getRHSvector() ##returns RHS term built by fipy
###Analytic solution for problem
###Ogata-Banks solution for 1-D ADE with constant influx term
###modified slightly (superposed) to account for baselinc conc =1.
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]
