Hi, guys:
Thanks a lot for all the suggestions. It really kind of you!
As a start, I tried to add a source and a sink term in the equation to
mimic the boundary condition with the following script:
c_source = 0
c_sink = 1
source_value = CellVariable(mesh=mesh, value=1e+10 * ((x-0.5)**2+
(y-1)**2<0.01))
sink_value = CellVariable(mesh=mesh, value=1e+10 * ((x-2.5)**2+
(y-1)**2<0.01))
# pde
D = 1./beta
eq = (TransientTerm() == ImplicitDiffusionTerm(coeff=D)
-
PowerLawConvectionTerm(coeff=pot.getFaceGrad())
+
ImplicitSourceTerm(coeff=(pot.getFaceGrad()).getDivergence())
- ImplicitSourceTerm(coeff=sink_value) +
c_sink * sink_value
- ImplicitSourceTerm(coeff=source_value) +
c_source * source_value
)
However, the result I got does not really make sense. It does not have
the symmetry along the y axis as I expected. I would appreciate if you
could take your professional eye and see if you could catch anything
obvious wrong.
Thanks a lot,
Bin
The complete script is:
########################################
from fipy import *
# system parameter
beta = 1.
# define mesh
n = 1
nx = 100*n
ny = 100*n
dx = 0.03/n
dy = 0.02/n
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
# define variable
phi = CellVariable(name='committor', mesh=mesh, value=0.0)
pot = CellVariable(name='potential', mesh=mesh, value=0.0)
x, y = mesh.getCellCenters()
pot.setValue(2.5*((x-1.5)**2-1)**2 + 5*(y-1)**2)
# boundary condition
#all_faces = mesh.getExteriorFaces()
BCs = (FixedFlux(faces=mesh.getFacesTop(), value=0),
FixedFlux(faces=mesh.getFacesBottom(), value=0),
FixedFlux(faces=mesh.getFacesLeft(), value=0),
FixedFlux(faces=mesh.getFacesRight(), value=0),
)
c_source = 0
c_sink = 1
source_value = CellVariable(mesh=mesh, value=1e+10 * ((x-0.5)**2+
(y-1)**2<0.01))
sink_value = CellVariable(mesh=mesh, value=1e+10 * ((x-2.5)**2+
(y-1)**2<0.01))
# pde
D = 1./beta
eq = (TransientTerm() == ImplicitDiffusionTerm(coeff=D)
-
PowerLawConvectionTerm(coeff=pot.getFaceGrad())
+
ImplicitSourceTerm(coeff=(pot.getFaceGrad()).getDivergence())
- ImplicitSourceTerm(coeff=sink_value) +
c_sink * sink_value
- ImplicitSourceTerm(coeff=source_value) +
c_source * source_value
)
eq.solve(var=phi, boundaryConditions=BCs)
viewer = Viewer(vars=(pot,), datamin=0., datamax=5)
viewer.plot()
raw_input("implicit transient diffusion. press <return> to proceed...")
viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)
viewer.plot()
raw_input("Implicit transient diffusion. Press <return> to proceed...")
On Oct 13, 2010, at 12:56 AM, Benny Malengier wrote:
2010/10/12 BIN ZHANG <[email protected]>
Dear Daniel:
Thanks a lot for your suggestions. Now I'm actually able to solve
the problem with normal boundary conditions. But I still have one
extra question, is it possible for me to use complicated boundaries
like the one shown in the attachment. For the white circle on the
left, I would have phi=0, and for the white circle on the right, I
would have phi=1.
I tried to define these boundaries using mesh.getInteriorFaces(),
but I got an error about: "Face list has interior faces"
Is there a way to enforce these complicated boundaries in fipy?
Yes, eg we use gmsh to create a domain with holes, and then have
code like
face_in = ((self.mesh2d.getExteriorFaces()) &
(sp.power(xfc,2) + sp.power(yfc,2) \
< (self.radius_domain -
self.radius_boundlayer)**2))
face_ex = (~face_in) & (self.mesh2d.getExteriorFaces())
So face_in are the interior faces in this case defined as the
exteriorfaces that satisfy being in a certain domain, and the other
faces is a negation of that. Then face_in and face_ex can be used to
set boundary conditions. xfc and yfc are the face centers.
Benny
Thanks a lot,
Bin
<grid.png>
On Oct 12, 2010, at 10:43 AM, Daniel Wheeler wrote:
On Tue, Oct 12, 2010 at 1:27 PM, BIN ZHANG <[email protected]> wrote:
Dear Daniel:
Thanks a lot for your reply. Greatly appreciated.
The trick you used to transform the equation is quite clever ;-).
I guess
now my question is how to represent $\partial_i V$ in the
convection term.
Try,
ConvectionTerm(V.getFaceGrad())
Though you confirmed to me that it's possible to have a convection
term with
a variable coefficient, like $\partial_i V$, I'm still not clear
how to do
that in Fipy. From the documentation on Convectionterm, I found the
following:
>>> cv = CellVariable(mesh = m)
>>> __ConvectionTerm(coeff = cv)
Traceback (most recent call last):
...
TypeError: The coefficient must be a vector value.
It looks to me that only, a scalar, or a constant vector is
allowed in this
term. Any suggestion?
The coeff for ConvectionTerm needs to be a vector. Try,
cv = CellVariable(mesh = m, rank=1)
__ConvectionTerm(coeff = cv)
"rank=1" makes the variable a vector value.
--
Daniel Wheeler
