Thank you. I am busy for the week, but then will get a chance to think about this :-)
Cheers, Jamie On Tue, Oct 4, 2016 at 1:16 PM, Daniel Wheeler <daniel.wheel...@gmail.com> wrote: > Hi James, > > Thanks for this example. I think I understand what is happening. Let's > start with the code for solving the problem, > > ~~~~ > import fipy > import numpy > > mesh = fipy.Grid1D(nx=10) > > var = fipy.CellVariable(mesh=mesh) > var[:] = 2 * numpy.random.random(mesh.numberOfCells) > var.constrain(1., where=mesh.facesRight) > var.faceGrad.constrain(0, where=mesh.facesLeft) > > eqn = fipy.ConvectionTerm([-1]) == fipy.DiffusionTerm(1.0) > > for _ in range(100): > eqn.sweep(var) > print var > ~~~ > > Note that the "var.faceGrad.constrain" is required with convection > terms to make sure that there is an outlet flow. This is a weird > decision in many ways, but was made for backwards compatibility with > previous FiPy versions, but this is not the main issue here. > > Notice that when the diffusion coefficient is small or large then this > system solves in one sweep or close to one sweep (set it to 1000.0 or > 0.0001 for example). The number of sweeps required to solve the system > is at its worst when the diffiusion and convection coefficients are > close in magnitude. So what is going on? Basically, the system is not > fully implicit. > > The left hand boundary condition is dependent on the left hand cell > value but gets added to the b vector in the case of non upwind > convection terms. In the case of fully upwind convection terms this > contribution is zero so as the diffusion coefficient goes to zero the > CovectionTerm becomes fully upwind (or if a UpwindConvectionTerm is > selected). As the diffusion coefficient gets large the contribution of > the convection term decreases and thus the contribution from the > explicit boundary condition is reduced. > > The discretization for the left most cell for a central difference > scheme (we're using power law above, but it has some central > difference contribution at D=1 and u=-1) is > > D * phi0 - (u / 2 + D) * phi1 = -u * phiB / 2 > > where phi0 is the value at the left most cell and phi1 is the cell to > the right. In FiPy, phiB (the boundary value of phiB) is not > calculated implicitly, but is explicitly set for the convection term > boundary condition. Basically phiB is the previous value of phi0 in > the FiPy code. That may not be the best way to deal with this. It may > be that it needs to be explicit for stability reasons. > > A possible way to deal with this issue is to use terms for the outlet > condition that are upwind at the outflow, but power law everywhere > else. > > Anyway, I hope that helps with understanding what's going on even if > the implementation is sub-optimal. > > Cheers, > > Daniel > > On Thu, Sep 15, 2016 at 5:13 PM, James Pringle <jprin...@unh.edu> wrote: > > Dear FiPy users -- > > > > This is a simple example of how and why fipy may fail to solve a > > steady advection diffusion problem, and how solving the transient > > problem can reduce the error. I also found something that was a > > surprise -- the "initial" condition of a steady problem can affect > > the solution for some solvers. > > > > At the end are two interesting questions for those who want to > > understand what FiPy is actually doing.... I admit to being a bit > > lost > > > > The equation I am solving is > > > > \Del\dot (D\del psi + u*psi) =0 > > > > Where the diffusion D is 1, and u is a vector (1,0) -- so there is > > only a flow of speed -1 in the x direction. I solve this equation > > on a 10 by 10 grid. The boundary conditions are no normal gradient > > on the y=0 and y=10 boundary: > > > > psi_y =0 at y=0 and y=10 > > > > For the x boundary, I impose a value of x=1 on the inflow boundary at > > x=10 > > (this is a little tricky -- the way the equation is written, u is > > the negative of velocity). > > > > psi=1 at x=10 > > > > and a no-normal-gradient condition at x=0. > > > > psi_x=0 at x=0 > > > > since all of the domain and boundary is symmetrical with respect to > > y, we can assume psi_y=0 is zero everywhere. This reduces the > equation to > > > > psi_xx + psi_x =0 > > > > The general solution to this equation is > > > > psi=C1+C2*exp(-x) > > > > Where C1 and C2 are constants. For these boundary conditions, C1=1 > > and C2=0, so psi=1 everywhere. > > > > Now run the code SquareGrid_HomemadeDelaunay and look at figure(3) > > -- this is the plot of psi versus x, and you can see that it does > > not match the true solution of psi=1 everywhere! Instead, it > > appears to be decaying exponential. The blue line is actually just > > (1+exp(-x)). What is going on? It appears that small errors in the > > boundary condition at x=0 are allowing C2 to not be exactly 0, and > > this error is this exponential mode. The error is the artificial > > exiting of a correct mode of the interior equation, albeit one that > > should not be excited by these BC's. > > > > Interestingly, the size of this error depends on the value the psi > > is initially set to. If the line > > > > psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0) > > > > is changed so psi is initially 1, the error goes away entirely; if > > it is set to some other value, you get different errors. I do not > > know if this is a bug, or just numerical error in a well programmed > > solver. > > > > Now if you run SquareGrid_HomemadeDelaunay_transient which > implements > > > > psi_t = \Del\dot (D\del psi + u*psi) > > > > you can see that the error in the numerical solution is advected > > out of the x=0 boundary, and the solution approaches the true > > solution of psi=1 rapidly. > > > > The interesting question is if the formulation of the boundary > > condition at x=0 could be altered to less excite the spurious mode? > > > > Also, why does the "initial" condition have any effect on the > > steady equation? > > > > Cheers, > > Jamie > > > > > > _______________________________________________ > > fipy mailing list > > fipy@nist.gov > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www. > ctcms.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r= > 7HJI3EH6RqKWf16fbYIYKw&m=tZaHEx9flGw_0A9dbKvkhX2A2ZKYl0jHC2CQn6iZFO > 4&s=lIeSkvrgSp7C56l_JXBrrtJGnWXSYmtREoARHhZ8gYI&e= > > [ NIST internal ONLY: https://urldefense.proofpoint. > com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_ > fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m= > tZaHEx9flGw_0A9dbKvkhX2A2ZKYl0jHC2CQn6iZFO4&s=nvZbf- > 2vJlO4vBbZz4Vwi2RrSieAIr7clFYllQqV4ng&e= ] > > > > > > -- > Daniel Wheeler > _______________________________________________ > fipy mailing list > fipy@nist.gov > https://urldefense.proofpoint.com/v2/url?u=http-3A__www. > ctcms.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r= > 7HJI3EH6RqKWf16fbYIYKw&m=tZaHEx9flGw_0A9dbKvkhX2A2ZKYl0jHC2CQn6iZFO > 4&s=lIeSkvrgSp7C56l_JXBrrtJGnWXSYmtREoARHhZ8gYI&e= > [ NIST internal ONLY: https://urldefense.proofpoint. > com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_ > fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m= > tZaHEx9flGw_0A9dbKvkhX2A2ZKYl0jHC2CQn6iZFO4&s=nvZbf- > 2vJlO4vBbZz4Vwi2RrSieAIr7clFYllQqV4ng&e= ] >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
