David - The problem is here:
diffCoeff=Variable(value=1.) diffCoeff.constrain(0., mesh.facesLeft) It doesn't mean anything to constrain a Variable at particular faces. A Variable doesn't exist on a mesh; it's just a value. You can either do diffCoeff=CellVariable(mesh=mesh, value=1.) diffCoeff.constrain(0., mesh.facesLeft) or diffCoeff=Variable(value=1.) and leave off the constraint. Taking the divergence 51 faces to get 50 cells is not a problem. DiffusionTerm is declared right here on the page you linked: https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.terms.html#fipy.terms.diffusionTerm.DiffusionTerm DiffusionTerm is an implicit term, which is what you usually want. ExplicitDiffusionTerm is generally not what you want. > On Jun 4, 2018, at 11:13 PM, Ollodart, David Bernd <david...@illinois.edu> > wrote: > > Hello, > > I have a reaction-diffusion equation which has a flux boundary condition > dependent on the potential value at the boundary. I have seen in > examples.diffusion.mesh1d how to set the boundary conditions to a constant or > function of independent variables, but trying to insert the dependent > variable, or its face value (evaluated at some face(s)) directly into > .constrain does not work. I have seen at > https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.mail-archive.com%2Ffipy%40nist.gov%2Fmsg00389.html&data=02%7C01%7Cjonathan.guyer%40nist.gov%7Cb48df7e9f52d4b681f4608d5ca9277f3%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C636637652788127952&sdata=5f05NC56AzhaVB%2BGKg9gncSQTlDSsUQugUfNQzYNRDs%3D&reserved=0 > a question exactly on this topic but it uses FixedFlux which is now > deprecated: https://www.ctcms.nist.gov/fipy/documentation/USAGE.html. But I > have problems trying to implement the new methods suggested under "Applying > fixed flux boundary conditions". By the way,! there seems to be a typo; it should read ExplicitDiffusionTerm(diffCoeff) rather than DiffusionCoeff(diffCoeff), there being no DiffusionCoeff() in the fipy.terms package: https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.terms.html. > > I am first trying to use a constant external flux in the new method, which > could be achieved by potential.faceGrad.constrain(1., mesh.facesLeft), then > extend it to a potential-dependent flux. I have modified > examples.diffusion.mesh1d as below. I define the exterior flux on the whole > interior mesh, because the product with mesh.facesLeft makes it only non-zero > on the left face. The script is: > --------------------------------------------------------------------------------------- > from fipy import * > from scipy.special import erf > > #defining mesh of independent variables > nx = 50 > dx = 1. > mesh = Grid1D(nx=nx, dx=dx) > timeStepDuration = 0.9*dx**2/(2*1) > steps = 10 > x = mesh.cellCenters[0] > t = timeStepDuration * steps > > #setting parameter values > diffCoeff=Variable(value=1.) > diffCoeff.constrain(0., mesh.facesLeft) > > #defining solution variable and known analytical > phi = CellVariable(name="solution variable", > mesh=mesh, > value=0.) > > #defining boundary conditions > phi.constrain(1.,mesh.facesRight) > exteriorFlux = FaceVariable(name="fixed flux", > mesh=mesh, > value=1., > rank=1) > > #defining governing equations (might be matrix or vector valued for multiple > dependent variables) > eqn = TransientTerm() == ExplicitDiffusionTerm(coeff=diffCoeff) + > (mesh.facesLeft * exteriorFlux).divergence > > #defining the viewer > viewer = Viewer(vars=phi, > datamin=0.,datamax=1.) > > #iterating through solution > for step in range(steps): > eqn.solve(var=phi,dt=timeStepDuration) > viewer.plot() > > raw_input("Diffusion profile. Press <return> to proceed...") > -------------------------------------------------------------------------------------- > However, this throws an IndexError: too many indices for array. If i print > the product of exteriorFlux and mesh.facesLeft, it is an array of 51 > elements, the leftmost being 1, the others being 0, which is 1 greater than > the size of the array found from mesh.cellCenters which is what the > CellVariable is defined on. This error results if one specifies flux on one > or both sides of the mesh, the latter being the example provided. Since there > is 1 more face than there are cells, this result might be expected for > exteriorFlux, but then which mesh should be chosen? Trying to define a > different mesh throws an error for the product mesh.facesLeft * exteriorFlux. > > Thank you, > > David Ollodart > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
