Created
September 23, 2019 19:53
-
-
Save wd15/54f6d412a6f5e9f5ad5abd2f73549842 to your computer and use it in GitHub Desktop.
Stackoverflow answer, https://stackoverflow.com/questions/57937418/how-to-write-equation-terms-with-independent-variable-in-fipy
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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