It seems there's a bit of confusion here about initial condition vs. source term. Some comments in line. Hope it helps.
Ray On Sun, Apr 3, 2016 at 8:38 AM, Dario Panada <[email protected]> wrote: > Hi, > > I originally asked this on > http://stackoverflow.com/questions/36385367/solving-pde-with-sources-in-python > but am copying it to this list in case somebody has any advice :) I've > already read through > https://www.mail-archive.com/[email protected]/msg02439.html but didn't solve > my problem. > The list archives can be a good source, and I'd probably even more highly recommend the examples on the fipy webpage in terms of learning per time. > > I'm using FiPy to address a problem inspired by Biology. > > Essentially, I want to represent a 2D plane where at different points I > have sources and sinks. Sources emit substrate at a fixed rate (different > sources can have different rates) and sinks consume substrate at a fixed > rate (different sinks can have different rates). My code: > > import numpy.matlib > from fipy import CellVariable, Grid2D, Viewer, TransientTerm, > DiffusionTerm, ImplicitSourceTerm, ExplicitDiffusionTerm > from fipy.tools import numerix > from time import * > > nx = 10 > ny = nx > dx = 1. > dy = dx > L = dx*nx > mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny) > > arr_grid = numerix.array(( > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0,),'d') > > arr_source = numerix.array(( > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0.5,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0,),'d') > > arr_sink = numerix.array(( > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0, > 0,0,0,0,0,0,0,0,0,0.5,),'d') > > source = CellVariable(mesh=mesh, value = arr_source) > sink = CellVariable(mesh=mesh, value = arr_sink) > phi = CellVariable(name = "solution variable", mesh = mesh, value = > arr_grid) > As arr_grid is simply zeros, you could write phi = CellVariable(name="solution variable", mesh=mesh, value=0.) and avoid arr_grid entirely. > X,Y = mesh.cellCenters > phi.setValue(3.0, where=(X < 2.0) & (X > 1.0)) > phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0)) > You might consider setting up your source/sink CellVariables with this same sort of "where" logic so that you avoid writing full arrays out explicitly. What you've done with this line is establish an initial condition for the field variable (or an initial guess if you're trying to do a steady-state solve directly). Also, because this is an initial condition -- assuming phi is representing, e.g., concentration -- then, it's strange that phi is set to negative. I think this gets back to the issue that this is an initial condition and not a source/sink term. > D = 1. > eq = TransientTerm() == DiffusionTerm(coeff=D) > There is no source/sink in this equation, so the steady profile will be a perfectly uniform system with the average system concentration defined by the initial conditions. Actually, I'm not sure how FiPy treats the steady-state initial guess for Laplace's equation with no flux boundary conditions like yours here. The governing equation + BC's without an initial condition admits any uniform profile as a solution. Anyway, the approach you took below does have the source/sink included in the actual equations, so that seems like the direction to pursue. > > > > viewer = Viewer(vars=phi, datamin=0., datamax=1.) > > steadyState = False > > if(steadyState): > print("SteadyState") > DiffusionTerm().solve(var=phi) > viewer.plot() > sleep(20) > else: > print("ByStep") > timeStepDuration = 10 * 0.9 * dx**2 / (2 * D) > steps = 500 > for step in range(steps): > print(step) > eq.solve(var=phi, > dt=timeStepDuration) > if __name__ == '__main__': > viewer.plot() > > This works well, but FiPy treats the sources as "non-renewable" and > eventually I get an homogeneous concentration throughout the space as would > be expected. The alternative was to delete: > Here, I'm not sure what you mean by "non-renewable." As written above, there isn't any source/sink at all. > > X,Y = mesh.cellCenters > phi.setValue(3.0, where=(X < 2.0) & (X > 1.0)) > phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0)) > > And change the equation to: > > eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink > > Given that source and sink are never changed this offered "infinite" > sources and sinks. However, when I try solving for steady-state using > > eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink > I don't see any difference here. Did you forget to edit this line? To try to solve directly for steady state, you could write eq = DiffusionTerm(coeff=D) + source - sink then run eq.solve() noting that you call the eq solve method, not one from only the DiffusionTerm as written above. I'm still unsure about the treatment of the initial phi values in this sense as mentioned above. However, the direct-to-steady approach merits a few words of caution. First, there isn't really a good numerical way of directly computing steady state solutions for general systems. Often, your best bet is actually to solve the transient equation by time stepping from some initial condition until you reach steady state, as that's actually probably the most robust algorithm for solving for steady state profiles. Second, are you sure that your system admits a steady state? You have no-flux boundary conditions, but you have internal sources/sinks. Unless those sources/sinks produce/consume the field variable at exactly the same rate, you will have net accumulation in the system. Simple example with R = constant and uniform, non-zero value: dc/dt = R with no flux boundary conditions is an equation that doesn't admit any steady state. The addition of a diffusion term doesn't change that. If you were to add any dirichlet (specify value) boundary conditions, that would enable a steady solution, as the net production/consumption within the system could leave/enter via the system boundaries. > > I get: > > C:\Python27\python.exe > C:/Users/dario_000/PycharmProjects/mesa-test/mesher.py > SteadyState > > C:\Python27\lib\site-packages\fipy-3.1.dev134+g64f7866-py2.7.egg\fipy\solvers\scipy\linearLUSolver.py:71: > RuntimeWarning: invalid value encountered in double_scalars > if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= > self.tolerance: > > And the equation isn't solved. However if I solve it "by steps" using > again: > > eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink > > I get a nice picture similar to what I'd be expecting: > > http://i.stack.imgur.com/OsJda.png > > Any advice on how I can specify an initial setup with sources/sinks at > different spatial positions each with a different emission/consumption rate > in such a way that I may be able to obtain the steady-state solution? > > Thanks! > > > > > _______________________________________________ > 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 ]
