Skip to content

Instantly share code, notes, and snippets.

@wd15
Created September 23, 2019 19:53
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save wd15/54f6d412a6f5e9f5ad5abd2f73549842 to your computer and use it in GitHub Desktop.
Save wd15/54f6d412a6f5e9f5ad5abd2f73549842 to your computer and use it in GitHub Desktop.
from fipy import Variable,FaceVariable,CellVariable,Grid1D,ImplicitSourceTerm,TransientTerm,DiffusionTerm,Viewer,ConvectionTerm
from fipy.tools import numerix
import math
D = 1
c0 = 1
ka = 1
r0 = 1
nx = 100
dx = 2*math.pi/100
mesh = Grid1D(nx=nx, dx=dx)
conc = CellVariable(name="concentration", mesh=mesh, value=0., hasOld=True) # This is the "phi" in the docs
valueLeft = c0
valueRight = 0
conc.constrain(valueRight, mesh.facesRight)
conc.constrain(valueLeft, mesh.facesLeft)
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100
show_per_steps = 50
A = 1 / (r0**2 * numerix.sin(mesh.x)[0])
dA = -(numerix.cos(mesh.x)[0])/(r0**2 * numerix.sin(mesh.x)[0]**2)
dsindA = (numerix.cos(mesh.x)[0])**3/(numerix.sin(mesh.x)[0])**2
eqX = TransientTerm() + ImplicitSourceTerm(ka) == DiffusionTerm(D*A*numerix.sin(mesh.x)) - ConvectionTerm(D*dA*numerix.cos(mesh.cellCenters))+ D*conc*dsindA
sweeps = 3
from builtins import range
for step in range(steps):
conc.updateOld()
for sweep in range(sweeps):
res = eqX.sweep(var=conc, dt=timeStepDuration)
print(res)
if __name__ == '__main__' and step % show_per_steps == 0:
viewer = Viewer(vars=(conc), datamin=0., datamax=c0)
viewer.plot()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment