Daniel and Jonathan,
I've rerun the problem we worked on last Wednesday, and I'm having the
same problems that I described at that time. I'm not sure if the script
I'm running is different than the one we worked on at NIST, or if we were
a little premature in thinking the problem was fixed.
I'm still having problems with the convection term in my 3-D model. The
flow conditions should be at steady state, but some cells divergence of
the fluxes are far enough from zero to cause problems with the solution.
If I've set the problem up correctly, the max and min concentration values
should be 1700 and 180, respectively.
To confirm that the divergence is the problem, I added another
ImplicitSource to the equation (line 82) based on the cell divergence (to
add a flux assigned the cell conc, recirculate the water in the cell back
into the cell and balancing the flux). This produces concentrations that
are within a sensible range, but I've not checked this against an
analytic solution.
I'd note I get these problems what I add heterogeneity to the mesh by
varying the permeability of cells. Things work fine when all the flow is in
the Z direction (along the axis of the extruded tube), and not diverted in
the x or y directions (where the grid is not orthogonal).
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
from fipy import *
def addLayer(verts):
addArray = verts.copy()
addArray[:] = 0
addArray[2] = length/slices
return verts + addArray
length=40. ##cm
slices=40
mesh=GmshImporter2D('./semicircle.msh').extrude(layers=slices,
extrudeFunc=addLayer)
delZ=length/slices ##needed later in script
conc=CellVariable(name='solute conc',mesh=mesh,value=180.,hasOld=1)
###calculate head
head=CellVariable(name='head',mesh=mesh,value=110.,hasOld=1.)
Por=CellVariable(name='effective porosity',mesh=mesh,rank=0,value=.2)
SpStore=CellVariable(name='specific storage',mesh=mesh,rank=0,value=0.)
kond=CellVariable(name='conductivity',mesh=mesh,value=1.e-4)
###generate random numbers
#Dist=random.uniform(-2.,2.,mesh.getNumberOfCells())
###Set diffusion coeff
X,Y,Z=mesh.getCellCenters()
kond.setValue(1.e-6, where=((Z>20.)&(Z<30)&(X>2.54)))
#kond.setValue(1.e-4*(2**Dist))
Eq=TransientTerm(SpStore)==DiffusionTerm(kond.getHarmonicFaceValue())
###Set constant head. boundary conditions
X,Y,Z=mesh.getFaceCenters()
ZMax=max(Z)-0.005
bndry=(FixedValue(mesh.getFaces().where((Z<0.005)),120.),
FixedValue(mesh.getFaces().where((Z>ZMax)),100.))
###Calculate steady state distribution
delTime=1.e10
timeSteps=5
for step in range(timeSteps):
head.updateOld()
resid=1.
print 'head',step,max(head),min(head)
while resid>1.e-14:
resid=Eq.sweep(
var=head,
solver=LinearCGSSolver(iterations=2000,tolerance=1.e-17),
boundaryConditions=bndry,
dt=delTime,
underRelaxation=None)
print 'head',resid
##calculate diffusive flux
Flux=FaceVariable(name='specific discharge', mesh=mesh,
value=-head.getFaceGrad()*kond.getHarmonicFaceValue())
Disp=FaceVariable(name='dispersion tensor',mesh=mesh,rank=1)
Disp.setValue([1.e-5,1.e-5,1.e-5])
X,Y,Z=mesh.getFaceCenters()
bndryOn=(FixedValue(mesh.getFaces().where((Z<0.001)),1700.))
X,Y,Z=mesh.getCellCenters()
fluxInOn=CellVariable(name='inflow on',mesh=mesh)
ZValueIn=length/(slices*2.)+1.e-5
fluxInOn.setValue(1700.*Flux.getDivergence(),where=(Z<ZValueIn))
ZValueOut=length-length/(slices*2.)-1.e-5
fluxOut=CellVariable(name='outflow',mesh=mesh,value=0)
fluxOut.setValue(Flux.getDivergence(), where=(Z>ZValueOut))
fluxError=CellVariable(name='water balance error',mesh=mesh,value=0)
fluxError.setValue(Flux.getDivergence(), where=((Z>ZValueIn)&(Z<ZValueOut)))
conc=CellVariable(name='solute conc',mesh=mesh,value=180.,hasOld=1)
transportEqOn=(TransientTerm(coeff=Por)==
DiffusionTerm(coeff=Disp)
-UpwindConvectionTerm(coeff=Flux)
+fluxInOn
+ImplicitSourceTerm(coeff=fluxOut)
+ImplicitSourceTerm(coeff=fluxError))###should not need this
last term
delTime=600.
t=0
timeSteps=10
for step in range(timeSteps):
t+=delTime
print 'conc ',t,max(conc),min(conc)
conc.updateOld()
resid=1.
print step
while resid>1.e-14:
resid=transportEqOn.sweep(var=conc,
solver=LinearCGSSolver(iterations=500,
tolerance=1.e-16),
boundaryConditions=bndryOn,
dt=delTime)
print 'conc', resid
#viewer=make(conc,limits={'datamin':180.,'datamax':1700.})
#viewer.plot()
$MeshFormat
2 0 8
$EndMeshFormat
$Nodes
22
1 0 0 0
2 2.54 0 0
3 5.08 0 0
4 4.828460924472671 1.102064697357506 0
5 4.123664096722479 1.985851965467763 0
6 3.10520317224885 2.476316896941875 0
7 1.974796827747875 2.476316896941128 0
8 0.9563359032731413 1.985851965464269 0
9 0.2515390755258018 1.102064697354334 0
10 1.015999999999299 0 0
11 2.031999999999298 0 0
12 3.048000000000643 0 0
13 4.064000000001987 0 0
14 4.155014625466476 0.8569626475459882 0
15 3.389915582005075 1.421237035085163 0
16 2.584548425562228 0.9504419447165823 0
17 1.662631170398718 1.24745043415448 0
18 2.642308901053341 1.826504769976305 0
19 1.978148650661269 1.796132000536478 0
20 3.543715653916555 0.5986322973810793 0
21 1.098957930737896 0.9683454136911297 0
22 1.703842637684158 0.5902104244185791 0
$EndNodes
$Elements
43
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 5
6 1 3 0 1 0 5 6
7 1 3 0 1 0 6 7
8 1 3 0 1 0 7 8
9 1 3 0 1 0 8 9
10 1 3 0 1 0 9 1
11 1 3 0 2 0 1 10
12 1 3 0 2 0 10 11
13 1 3 0 2 0 11 12
14 1 3 0 2 0 12 13
15 1 3 0 2 0 13 3
16 2 3 0 4 0 6 7 18
17 2 3 0 4 0 7 8 19
18 2 3 0 4 0 19 8 17
19 2 3 0 4 0 13 14 20
20 2 3 0 4 0 20 14 15
21 2 3 0 4 0 18 16 15
22 2 3 0 4 0 19 18 7
23 2 3 0 4 0 16 18 19
24 2 3 0 4 0 19 17 16
25 2 3 0 4 0 12 13 20
26 2 3 0 4 0 15 16 20
27 2 3 0 4 0 8 9 21
28 2 3 0 4 0 21 17 8
29 2 3 0 4 0 12 20 16
30 2 3 0 4 0 22 10 11
31 2 3 0 4 0 21 10 22
32 2 3 0 4 0 17 21 22
33 2 3 0 4 0 16 17 22
34 2 3 0 4 0 12 16 11
35 2 3 0 4 0 11 16 22
36 2 3 0 4 0 5 15 14
37 2 3 0 4 0 14 4 5
38 2 3 0 4 0 10 21 9
39 2 3 0 4 0 9 1 10
40 2 3 0 4 0 6 18 15
41 2 3 0 4 0 15 5 6
42 2 3 0 4 0 14 13 3
43 2 3 0 4 0 3 4 14
$EndElements
