Hi Rose, I was trying to debug this and it seems like there is an inconsistency between the left hand boundary and the analytical solution.
By my calculations A~1 and B~0, that makes the gradient of the line about 1, yet the left hand boundary condition expects a gradient of about 0.0666 (if C is close to 0 on the left hand side, which the analytical solution wants). There may still be issues with FiPy, but I'd like to resolve that one first. I may also have gotten mixed up in the calculations, but I can't see where. Cheers, Daniel On Wed, Jun 25, 2014 at 4:55 PM, yuan wang <rose.w...@tufts.edu> wrote: > Hi, >> >> I tried to solve a problem with a Robin Boundary condition in fipy. The >> physical problem is very simple, but I could not get the numerical solution >> to match the analytical solution. >> >> [image: Inline image 2] >> >> Is there something wrong with the way I set up the boundary condition? >> I've copied my script below and bolded the boundary condition lines. Thanks >> for your help. >> >> Best regards, >> Rose >> ---------------------------------------------- >> >> from fipy import * >> >> L = 1. >> >> nx = 100 >> >> dx = L/nx >> >> mesh = Grid1D(nx=nx,dx=dx) >> >> zc = mesh.cellCenters[0] >> >> Dm = 8e-6 # moleculer diffusion in water column >> >> Dmp = 5e-6 # moleculer diffusion in porewater (corrected for tortuosity) >> >> Dbw = 1e-5 # bioturbation in porewater >> >> D = Dmp+Dbw >> >> >> Cw=1e-3 >> >> Cl=1. >> >> >> time = Variable() >> >> beta = 1e-3 >> >> #*abs(numerix.sin(w*time)) >> >> >> # steadystate solution >> >> A = beta*(Cl-Cw)/(D+beta*L) >> >> B = Cl-A*L >> >> Canal = CellVariable(name = "analytical solution steadystate", mesh=mesh, >> value = A*zc+B) >> >> >> C = CellVariable(name="concentration", mesh = mesh,value=0.) >> >> *C.faceGrad.constrain(beta*(C.faceValue-Cw)/D, mesh.facesLeft)* >> >> *C.constrain(Cl,mesh.facesRight)* >> >> eqI =TransientTerm()==ImplicitDiffusionTerm(coeff=D) >> >> >> viewer = Viewer(vars= (Canal,C)) >> >> dt = 864 >> >> while time()<86400: >> >> time.setValue(time() + dt) >> >> eqI.solve(var = C,dt=dt) >> >> viewer.plot() >> > > > > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > -- Daniel Wheeler
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
