Hello FiPy users,
I am trying to implement Robin boundary condition for the below equation:
$$\mu \phi - \nabla (D \nabla \phi) = S on \Omega (interior)a \phi + b D
\partial (n . \phi) = 0 on \partial \Omega (boundary)$$
I read the previous messages and examples and found three ways to define Robin
boundary condition: 1- phi.faceGrad.constrain() (here)
| |
| | | | | |
| examples.convection.robin — FiPy 3.1 documentation>>> if __name__ ==
'__main__':... restol = 1e-5... anstol = 1e-3... else:... restol = 0.5...
anstol = 0.15 |
| |
| View on www.ctcms.nist.gov | Preview by Yahoo |
| |
| |
2- ImplicitSourceTerm() (here)
| |
| | | | | |
| Notebook on nbviewerCheck out this Jupyter notebook! |
| |
| View on nbviewer.ipython.org | Preview by Yahoo |
| |
| |
3- (a + b * D * mesh.faceNormals * mesh.exteriorFaces).divergence (here)
I think the (1) and (3) are not appropriate ways for my problem, because the
face value of phi doesn't update and I can not see the boundary on the solver
coefficient matrix. So I try to implement it through ImplicitSourceTerm, but I
don't know how exactly implement it. I wrote a simple code for 8*8 grid, But it
seems I miss something and I couldn't get the expected result. Please give some
hint to implement Robin boundary for my equation.
#=============================================================import
os,os.pathos.environ['FIPY_SOLVERS']='scipy'import fipynumOfAxNode = 8delta =
1.R_psi = 0.507683226223R_J = 0.365748652934806a = (1-R_psi)/4b = (1+R_J)/2mesh
= fipy.Grid2D(dx=delta, dy=delta, nx=numOfAxNode, ny=numOfAxNode)#Diffusion
CoefficientD = fipy.FaceVariable(name = "Diffusion Coefficient", mesh = mesh,
value = 0.7)#Absorption Coefficientmua = fipy.CellVariable(name = "Absorption
Coefficient", mesh = mesh, value = 0.2)mua.constrain(0., where=
mesh.exteriorFaces)#Source termS = fipy.CellVariable(name = "source variable",
mesh = mesh, value = 0.) S.constrain(1, where= ((mesh.x == 3.5) & (mesh.y ==
3.5)))#Field variable phi = fipy.CellVariable(name = "solution variable", mesh
= mesh, value = 0.) #Mask for boundary cells/facesmask_coeff =
(mesh.exteriorFaces * mesh.faceNormals).getDivergence()DB = 0.7 * a/b#equation
ieq = fipy.ImplicitSourceTerm(coeff=(mua+mask_coeff * DB)) -
fipy.DiffusionTerm(coeff=D) == S#coeff Matrix
eq.cacheMatrix()eq.solve(var=phi)coeMat = eq.matrixprint coeMat.matrix#boundary
face values boundaryValue =
phi.arithmeticFaceValue[mesh.exteriorFaces.value]print boundaryValue#plot field
valueviewer = fipy.MatplotlibViewer(vars=phi, title="Phi", limits={'datamin':
0., 'datamax': 1.})viewer.plot()viewer.axes.set_title("Solution variable
(steady state)")#=============================================================
Cheers,Vahid_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
