Jonathan,
Thanks for the suggests in the last e-mail. I've reworked and
simplified my problem, hopefully casting it into something that's
easier to understand. I've included my script below.
It seems my problem was I was not muliplying the flux terms by
mesh.faceNormals to get only fluxes normal to the cell faces. Once I
started doing this along with your other suggestions, the mass
balances all worked out.
Thanks again for your help.
P.S. Is there an introductory paper on phase field methods that you
could recommend? Something simple to get me started?
Andrew Reeve
###Mass balance example script###
from fipy import *
from fipy.variables.fixedBCFaceGradVariable import _FixedBCFaceGradVariable
'''
Calculate mass balance in entire model domain based on a simple fickian
diffusion model using a 3-D grid to solve a 2-D problem. Assign constant
concentration boundary conditions along left and right edge of model and
include a constant mass flux across top face of model.
constant flux
|| ||
\/ \/
Const. |--------------------|
Conc. | 10x5 grid | Const.
Bndry->| | Conc.
| |<-Bndry
|--------------------|
'''
### Set up very simple grid
nx=10;ny=1;nz=5
dx=[10. for i in range(nx)]
dy=[1. for i in range(ny)]
dz=[1. for i in range(nz)]
mesh=Grid3D(dx=dx,dy=dy,dz=dz,nx=nx,ny=ny,nz=nz)
conc=CellVariable(name='concentration [g/cm^3]',mesh=mesh,value=0.,hasOld=1)
diffuse=CellVariable(name='diffusion rate [cm^2/sec]',mesh=mesh,rank=0,value=1.e-5)
### define type 1 boundaries
Xf,Yf,Zf=mesh.getFaceCenters()
Xc,Yc,Zc=mesh.getCellCenters()
Bndry=[FixedValue((mesh.getExteriorFaces()&(Xf==Xf.min())),100.),
FixedValue((mesh.getExteriorFaces()&(Xf==Xf.max())),0.)]
### Assign a fixed flux of 1.e-5 grams/second/cm^2 across top of model
### Note that the flux is based on the surface area of the top cell face.
Source=FaceVariable(name='flux across top of model',mesh=mesh, rank=1,value=0.)
Source.setValue([Zf*0.,Zf*0,(Zf==Zf.max())*1.e-5])
Eqn=(TransientTerm(coeff=1.)
==DiffusionTerm(coeff=diffuse.getHarmonicFaceValue())
+Source.getDivergence())
delT=1.e3
for i in range(5):
conc.updateOld()
resid=1.
while resid>1.e-12:
resid=Eqn.sweep(var=conc,
solver= LinearGMRESSolver(tolerance=1.0e-16,
iterations=500),
dt=delT,
underRelaxation=None,
boundaryConditions=Bndry
)
print i,resid
type1Flux=(diffuse.getHarmonicFaceValue()*
_FixedBCFaceGradVariable(conc,boundaryConditions=Bndry))*mesh.faceNormals
###calculate mass change in model based on change in mass stored in cells
storedMass=(conc-conc.getOld())*mesh.getCellVolumes()
print 'change in stored mass,',sum(storedMass())
print 'mass flux due to
source,',sum(Source.getDivergence()*mesh.getCellVolumes())*delT
print 'mass fluc due to type1
bndry,',sum(sum(type1Flux*mesh.faceAreas*mesh.getExteriorFaces()))*delT
