Daniel, Jonathan, and FiPyer's,
I'm a little confused about results I'm getting from getDivergence in a
FiPy script.
I've created a triangular mesh with gmsh and extruded it with fipy.
I'm solving a simple diffusion equation, with a diffusion coefficient that
varies spatially, and running the model to a steady-state condition. In my
current test, I'm just varying the diffusion coefficient assigned to a
cell based on a uniform random distribution (from 0.5e-5 to 1.5e-5).
At the end to the simulation, I calculate the the diffusive fluxes across
cell faces based on Fick's Law. To check for mass-balance errors, I
calculate the divergence in the cells. As you know, the divergence
under steady-state conditions should be zero (or a very tiny number).
However, when I do this in FiPy the divergences are relatively large
(about 1e-7). In the script I've separated the divergence calculated at
the ends of the model with Type 1 boundaries (which shouldn't be zero due
to the influence of the boundary conditions), from the 'internal' cells.
The internal cells are only one or two orders of magnitude smaller than
those impacted by boundary conditions.
I realize I don't need to use sweeps to solve this set of equations, but
have left this in the script as I believe I'll need this when I add
anisotropy. I don't think it is impacting the results.
Am I missing something here, or does the large divergence in this
scenario indicate there are large mass-balance errors in each cell? Am
I calculating the diffusive flux correctly? I've tried a variety of
different methods of calculating this Flux term, but what is in the script
now makes the most sense to me.
I've also run (and commented out in the script) a scenario with a regular
mesh and this seems to yield better results (internal cells with a
divergence of 1.e-18). The equation also works fine if Diff doesn't very
spatially.
Hopefully I'll remember to attach the script, but here's a quick rundown
of what I'm doing:
1) create mesh
mesh=GmshImporter2D('./semicircle.msh').extrude(layers=slices,
extrudeFunc=addLayer)
2) define cellVariables: conc, por, Diff
3) vary Diff spatially:
Dist=random.uniform(.5,1.5,(mesh.getNumberOfCells()))
Diff.setValue(Dist*1.e-5)
4) define diffusion Equation:
Eq=TransientTerm(Por)==DiffusionTerm(Diff.getHarmonicFaceValue())
5) solve for conc. using sweeps, with const conc boundaries set at two
ends
6) calculate diffusive flux
Flux=conc.getFaceGrad()*Diff.getHarmonicFaceValue()
7) calculate divergence and make list with divergences at 'internal'
cells.
DivFlux=Flux.getDivergence()
X,Y,Z=mesh.getCellCenters()
InDivFx=[DivFlux()[i] for i in xrange(mesh.getNumberOfCells()) if
(Z[i]>delZ) and (Z[i]<length-delZ)]
Eventually, I'd like to use a diffusion-type equation to calculate
and assign fluxes to the advection-diffusion equation. If there are large
mass-balance errors in the initial set of calculations, I'll get
unreasonable results for the ADE equation...
Andy
----------------------
Andrew Reeve
Dept. of Earth Science
University of Maine
Orono, ME 04469
from fipy import *
###create regular mesh
## length=40. ##cm
## slices=40
## delZ=length/slices
## mesh=Grid3D(nx=8,ny=4,nz=slices,dx=.5,dy=.5,dz=delZ)
###create irregular mesh
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
###calculate concentration
conc=CellVariable(name='concentration ',mesh=mesh,value=100.,hasOld=1.)
Por=CellVariable(name='effective porosity',mesh=mesh,rank=0,value=.2)
Diff=CellVariable(name='concentration ',mesh=mesh)
###generate random numbers
Dist=random.uniform(.5,1.5,(mesh.getNumberOfCells()))
###Set diffusion coeff
Diff.setValue(Dist*1.e-5)
Eq=TransientTerm(Por)==DiffusionTerm(Diff.getHarmonicFaceValue())
###Set constant conc. 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.e12
timeSteps=10
for step in range(timeSteps):
conc.updateOld()
resid=1.
print 'conc',step
while resid>1.e-14:
resid=Eq.sweep(
var=conc,
solver=LinearGMRESSolver(iterations=1000,tolerance=1.e-16),
boundaryConditions=bndry,
dt=delTime)
print resid
###calculate diffusive flux
Flux=conc.getFaceGrad()*Diff.getHarmonicFaceValue()
DivFlux=Flux.getDivergence()
X,Y,Z=mesh.getCellCenters()
###divergence at internal cells
InDivFx=[DivFlux()[i] for i in xrange(mesh.getNumberOfCells()) if (Z[i]>delZ)
and (Z[i]<length-delZ)]
OutDivFx=[DivFlux()[i] for i in xrange(mesh.getNumberOfCells()) if (Z[i]<delZ)
or (Z[i]>length-delZ)]
$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
