Hi Andy,
I have a file in 2D that I used for a kinetics class this year. I
attach the code here for your
use. I will upload it on the Wiki in the near future.
Note that I've tried to thoroughly comment, however, I make no
guarantee that I've made this
file easy to use.
There is a particular feature of they way FiPy stores cell variables
that leads to the use of a
function called "take". I trolled the user manual to find an example
of this from the phase
field code, but, I don't know all the details for why this particular
function works correctly.
An update to FiPy that would make this easier to implement would be
welcome.
If you use this code and or make improvements (or find errors), I'd
be happy to hear about it
and/or use it next year.
Dan Lewis
Assistant Professor
Materials Research Center, Room 110
Rensselaer Polytechnic Institute
110 8th Street
Troy, NY 12180
[EMAIL PROTECTED]
518-276-2297
http://www.rpi.edu/~lewisd2
#!/usr/bin/env python
##
# ###################################################################
# MTLE-6963: Advanced Kinetics
# Rensselaer Polytechnic Institute
# Professor Dan Lewis, Spring 2007
#
# Laboratory #2: Anisotropic Diffusion Solutions
#
# This lab
#
#
#
#
#
# Objectives:
#
# Read through the comments in the code to
#
# 1)
#
# 2)
#
# 3)
#
# ###################################################################
##
# NOTE: The following lines import parts of fipy and numerical
routines. You
# should not need to change these items.
#
import Numeric
import fipy.viewers
from fipy.meshes.grid2D import Grid2D
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.boundaryConditions.fixedValue import FixedValue
from fipy.variables.cellVariable import CellVariable
from fipy.variables.vectorCellVariable import VectorCellVariable
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.powerLawConvectionTerm import PowerLawConvectionTerm
from fipy.terms.transientTerm import TransientTerm
# NOTE: The parameters of our system.
#
# 1) L is the LENGTH of our box. You can set this to any value you
choose
# however, 1 is the easiest.
#
# 2) nx is the number of grid points we wish to have in our
solution. If we
# have fewer we sacrifice accuracy, if we have more, the
computational time
# increases. You should always check that your solution does not
depend on
# the number of grid points and the grid spacing!
#
# 3) dx is the spacing between grid points. Similar comments as
above in #2
#
# 4) timeStepDuration is the amount of 'time' at each step of the
solution.
# Accuracy and stability of the solution depend on choices for the
timeStepDuration
# and the grid point spacing and the diffusion coefficient. We will
not discuss
# stability any further, just know that it is something that needs to
be considered.
#
# 5) The number of steps to run before finishing the program.
#
# 6) D{ij} are the anisotropic diffusivity entries in the diffusion
tensor.
L = 1.
nx = 60
ny = 60
dx = L / nx
dy = L / ny
timeStepDuration = 0.0001
steps = 30
D11 = 1.0
D12 = 0.0
D22 = 10.0
# NOTE: The 'mesh' command creates the mesh (gridpoints) on which we
will
# solve the equation. This is specific to FiPy.
mesh = Grid2D(dx = dx, dy = dy, nx = nx, ny = ny)
# NOTE: We assign to 'c' objects that are 'CellVariables'. This is
a special
# type of variable used by FiPy to hold the values for the concentration
# in our solution.
c = CellVariable(name = "c", mesh = mesh)
D = VectorCellVariable(name = "D", mesh=mesh, value=(D11,D22))
dcdydx = VectorCellVariable(name = "dcdydx", mesh=mesh, value=(0,0))
# NOTE: This command sets 'x' to contain a list of numbers that
define the
# x position of the grid-points. (The centers of the volume elements.)
x = mesh.getCellCenters()[...,0]
y = mesh.getCellCenters()[...,1]
c.setValue(0.0)
#
# NOTE: Use the following syntax if you wish to have diffusant in
your initial
# conditions. Don't forget to comment out the line above if you use
this form.
c.setValue(1.0, where=((x-0.5)**2.0 + (y-0.5)**2.0 < 0.03))
# NOTE: Set scaling for our output graphs.
viewer = fipy.viewers.make(vars = (c),
limits = {'datamin': 0., 'datamax': 1.})
# NOTE: Boundary conditions come in two types. Fixed flux and fixed
value.
# The syntax is:
#
# FixedValue(mesh.getFacesLeft(), VALUE)
# FixedFlux(mesh.getFacesLeft(), FLUX)
#
# Choose the edge where the boundary conditions are applied by using
mesh.getFacesLeft
# or mesh.getFacesRight
#
boundaryConditions=(FixedFlux(mesh.getFacesLeft(),0.0),
FixedFlux(mesh.getFacesRight(),0.0),
FixedFlux(mesh.getFacesTop(),0.0),
FixedFlux(mesh.getFacesBottom(),0.0)
)
viewer.plot()
raw_input("Initial condition. Press <return> to proceed...")
# To build of off diagonal terms we need a vector pointing in the "j"
direction with magnitude
# of the gradient of "c" in the "i" direction. This python code does
this exactly.
dcdydx = c.getFaceGrad()._take((1, 0), axis = 1)
# NOTE: If you want a transient solution you build your equation
this way:
eqn1 = TransientTerm() == ImplicitDiffusionTerm(D) + (D12 + D12)
*dcdydx.getDivergence()
#
# NOTE: The off diagonal terms require more complicated FiPy/Python/
Numeric interactions. I think
# these are right!
#
# NOTE: Note the use of '=' vs. '==' in the above line. You must
assign
# to 'eqn1' a BOOLEAN expression of the form A == B + C + ...
#
if __name__ == '__main__':
for step in range(steps):
c.updateOld()
eqn1.solve(c, boundaryConditions = boundaryConditions, dt =
timeStepDuration)
# Uncomment the two lines below if you want the solution to pause
after every 10th step.
#if step%10 == 0:
# raw_input('pause')
viewer.plot()
raw_input('finished')
# These lines print out a text file with all the data for your
plotting enjoyment.
# from fipy.viewers.tsvViewer import TSVViewer
# TSVViewer(vars=(c)).plot(filename="output.txt")
On Apr 9, 2007, at 6:40 PM, ASReeve wrote:
FiPy users,
Can tensors be assigned to Cell Variables in FiPy? I'd like to
assigning diffusion coefficients to cells with different Values for
D_xx and D_yy. I believe this was answered about a year ago on the
mailing list and wondered if this functionality might have been
added since then.
Andy
