What was the result you were expecting? As discussed in the email thread you linked to, (1) has a known bug: https://github.com/usnistgov/fipy/issues/426 I also discovered that the typeset analytical expression in the robin.py example was also wrong, although the analytical solution in the code is correct.
On Sep 28, 2015, at 11:36 AM, vahid ettehadi <[email protected]> wrote: > 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 nbviewer > Check 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.path > os.environ['FIPY_SOLVERS']='scipy' > import fipy > > numOfAxNode = 8 > delta = 1. > R_psi = 0.507683226223 > R_J = 0.365748652934806 > > a = (1-R_psi)/4 > b = (1+R_J)/2 > > mesh = fipy.Grid2D(dx=delta, dy=delta, nx=numOfAxNode, ny=numOfAxNode) > > #Diffusion Coefficient > D = fipy.FaceVariable(name = "Diffusion Coefficient", > mesh = mesh, > value = 0.7) > > #Absorption Coefficient > mua = fipy.CellVariable(name = "Absorption Coefficient", > mesh = mesh, > value = 0.2) > mua.constrain(0., where= mesh.exteriorFaces) > > #Source term > S = 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/faces > mask_coeff = (mesh.exteriorFaces * mesh.faceNormals).getDivergence() > DB = 0.7 * a/b > > #equation i > eq = fipy.ImplicitSourceTerm(coeff=(mua+mask_coeff * DB)) - > fipy.DiffusionTerm(coeff=D) == S > > #coeff Matrix > eq.cacheMatrix() > eq.solve(var=phi) > coeMat = eq.matrix > print coeMat.matrix > > #boundary face values > boundaryValue = phi.arithmeticFaceValue[mesh.exteriorFaces.value] > print boundaryValue > > #plot field value > viewer = 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
