Re: sources errors in advection/diffusion problems, and one solution
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 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 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
Re: sources errors in advection/diffusion problems, and one solution
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 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 diff
Re: sources errors in advection/diffusion problems, and one solution
Interesting. One addition to Zhekai's comment. When I read it, I thought Upwind might be doing better because it is adding additional numerical diffusion (as upwind usually does, at least in finite difference codes). However this is not the case. If you keep the advection u and diffusion D the 1D equation is D*psi_xx + u*psi_x =0 and the exponential term of the solution (the one that appears in the error) has the form exp((u/D)*x). An increase in diffusion D would make the error decay more slowly away from the boundary (you can test that this is so by increasing "DiffCoef" in the code). This increase in decay scale is not what appears when UpwindConvectionTerm is used for the convection term. Instead, the solution becomes perfect (Note -- in this context advection and convection are the same thing.) But it does give a hint as to the source of the error. The upwind schemes only use the gradient calculated in the direction the advection the current is coming from -- in this case u is from the positive x direction. This means an advection scheme in a finite-difference code (which is all that I know, I am new to finite element codes) would not usually use the gradient calculate on the x=0 boundary since the current is coming from +x). But the upwindConvection scheme is sensitive to the boundary condition on the x=0 boundary. (change the line "psi.faceGrad.constrain(0.0*faceNorm,where=xeq0bound)" to see this.) And it appears to be doing the right thing when the gradient is set to some other value than 0. So how does the upwind scheme implement the normal gradient boundary condition differently from the other schemes? At this point I should get off my lazy behind and dig through the code, but I have a number of deadlines for the rest of the week as I switch from being an incompetent applied mathematician to an incompetent writer and administrator. But if no one else has an idea, I might look next week. Jamie On Tue, Sep 20, 2016 at 6:52 PM, Zhekai Deng < zhekaideng2...@u.northwestern.edu> wrote: > Hi James, > > Thanks for providing this demo to illustrate the problem. I don't have any > particular ideas exactly why the initial value of psi helps to reduce the > error and why transient problem "advect" out. However, I have some > findings that may help on this. > > Finding # 1: I noticed you have tried ExponentialConvectionTerm, > PowerLawConvectionTerm, CentralDifferenceConvectionTerm, all of these > give exponential error. However, if you tried UpwindConvectionTerm, this > will give the right result on the steady state solution. Thus, maybe using > upwind method, the convection does not require the value from downstream, > consequently, the error from the downstream will not excite the error > toward the upstream. However, I am still very surprise with the magnitude > of the error from other methods, and how similar > > Finding # 2: In the transient state problem, if I increase the mesh size > from 50*50 to 100*100. The error actually grows larger for the > ExponentialConvectionTerm, PowerLawConvectionTerm, > CentralDifferenceConvectionTerm. To show this, if I change the mesh size > to 100*100, the max psi value I have is around 1.1 . However, if I change > the mesh size to 50*50, the max psi is 1.005366, which is several orders of > magnitude lower in terms of difference to the exact solution. This is also > the case for the UpwindConvectionTerm, however, the error for both mesh > size are very small (max(psi) = 1e-13 or 1e-14 + 1). So even in the > transient state, the mesh size appears to somehow amplify the error if we > use finer mesh. I am confused by this. > > To now, it seems UpWindConvectionTerm appears to the the work around > method to this issue. Of course, I will be willing to learn more about what > other people think on this. > > Best, > > Zhekai > > > > > On Fri, Sep 16, 2016 at 8:53 AM, James Pringle wrote: > >> No worries -- I had to do it to figure out the problem in my more complex >> domain and equation... The issue which surprised me was that the value the >> variable was initialized to had an effect on the steady solution. >> >> Jamie >> >> On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) < >> jonathan.gu...@nist.gov> wrote: >> >>> James - >>> >>> I think Daniel will have more insight into why this is happening and if >>> there is anything to be done about it besides artificial relaxation, but I >>> just want to say how much I appreciate your putting this together. This is >>> a very lucid illustration. >>> >>> - Jon >>> >>> > On Sep 15, 2016, at 5:13 PM, James Pringle 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. >>> > >>
Re: sources errors in advection/diffusion problems, and one solution
Hi James, Thanks for providing this demo to illustrate the problem. I don't have any particular ideas exactly why the initial value of psi helps to reduce the error and why transient problem "advect" out. However, I have some findings that may help on this. Finding # 1: I noticed you have tried ExponentialConvectionTerm, PowerLawConvectionTerm, CentralDifferenceConvectionTerm, all of these give exponential error. However, if you tried UpwindConvectionTerm, this will give the right result on the steady state solution. Thus, maybe using upwind method, the convection does not require the value from downstream, consequently, the error from the downstream will not excite the error toward the upstream. However, I am still very surprise with the magnitude of the error from other methods, and how similar Finding # 2: In the transient state problem, if I increase the mesh size from 50*50 to 100*100. The error actually grows larger for the ExponentialConvectionTerm, PowerLawConvectionTerm, CentralDifferenceConvectionTerm. To show this, if I change the mesh size to 100*100, the max psi value I have is around 1.1 . However, if I change the mesh size to 50*50, the max psi is 1.005366, which is several orders of magnitude lower in terms of difference to the exact solution. This is also the case for the UpwindConvectionTerm, however, the error for both mesh size are very small (max(psi) = 1e-13 or 1e-14 + 1). So even in the transient state, the mesh size appears to somehow amplify the error if we use finer mesh. I am confused by this. To now, it seems UpWindConvectionTerm appears to the the work around method to this issue. Of course, I will be willing to learn more about what other people think on this. Best, Zhekai On Fri, Sep 16, 2016 at 8:53 AM, James Pringle wrote: > No worries -- I had to do it to figure out the problem in my more complex > domain and equation... The issue which surprised me was that the value the > variable was initialized to had an effect on the steady solution. > > Jamie > > On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) < > jonathan.gu...@nist.gov> wrote: > >> James - >> >> I think Daniel will have more insight into why this is happening and if >> there is anything to be done about it besides artificial relaxation, but I >> just want to say how much I appreciate your putting this together. This is >> a very lucid illustration. >> >> - Jon >> >> > On Sep 15, 2016, at 5:13 PM, James Pringle 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 cha
Re: sources errors in advection/diffusion problems, and one solution
No worries -- I had to do it to figure out the problem in my more complex domain and equation... The issue which surprised me was that the value the variable was initialized to had an effect on the steady solution. Jamie On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) < jonathan.gu...@nist.gov> wrote: > James - > > I think Daniel will have more insight into why this is happening and if > there is anything to be done about it besides artificial relaxation, but I > just want to say how much I appreciate your putting this together. This is > a very lucid illustration. > > - Jon > > > On Sep 15, 2016, at 5:13 PM, James Pringle 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 > > > > transient.py>__ > _ > > 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=IaTL4nEetwjm4G4qfj8cwLzvztKzui > Fsw_Nhksv_oWQ&s=vG-nxTf76KxE_CqEHTjt2jIkoy0l9M6X8bm01ypXaBQ&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= > IaTL4nEetwjm4G4qfj8cwLzvztKzuiFsw_Nhksv_oWQ&s= > 7JGh0Zz82O69cJiLFKmkV3NfW2TLz6KB_ngkAGCrYGI&e= ] > > > ___ > 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=IaTL4nEetwjm4G4qfj8cwLzvztKzui > Fsw_Nhksv_oWQ&s=vG-nxTf76KxE_CqEHTjt2jIkoy0l9M6X8bm01ypXaBQ&e= > [ NIST internal ONLY: https://urldefense.proofpoint. > com/v2/url?u=https-3A__email.nist.gov_
Re: sources errors in advection/diffusion problems, and one solution
James - I think Daniel will have more insight into why this is happening and if there is anything to be done about it besides artificial relaxation, but I just want to say how much I appreciate your putting this together. This is a very lucid illustration. - Jon > On Sep 15, 2016, at 5:13 PM, James Pringle 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 > 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 ]
sources errors in advection/diffusion problems, and one solution
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 from pylab import * from numpy import * import fipy from collections import defaultdict import copy as Copy def ScipyDelaunay2mesh(points,simplices,debug=False): '''ScipyDelaunay2mesh(points,simplices): Creates a FiPy 2Dmesh object from the output of a scipy.spatial.Delaunay triangulation, Also returns a list of border faces for the application of boundary conditions. The list of border conditions is computed within this routine; it will handle the case where triangles/simplices have been removed from the output of the Delaunay() triangulation. Parameters -- points : array an (N,2) array of the location of N points that made the triangulation. If Delaunay is called with tri=scipy.spatial.Delaunay(), points is tri.points simplices : array an (M,3) array of the location of M triangles returned by Delaunay. If Delaunay is called with tri=scipy.spatial.Delaunay(), simplices is tri.simplices debug=False : Boolean if this is set to true, the mesh will be drawn in blue, boundaries will be drawn in red, boundary point numbers will be written at vertices, and numbers for triangles will be drawn in the triangle. This output is useful for finding where triangles have been removed incorrectly. Returns --- mesh: a FiPy mesh2D object ''' #Create the arrays that fipy.mesh.mesh2D needs. vertexCoords is #of shape (2,N), where N is the number of points, and contains the #coordinates of the vertices. It is equivalent to points in #content, but reshaped. print 'Entering ScipyDelaunay2Mesh; Making vertexCoords' vertexCoords=points.T #faceVertexID is of shape (2,M) where M is the number of faces/edges #(faces is the fipy terminology, edges is the scipy.spatial.Delaunay #terminology). faceVertexID contains the points (as indices into #vertexCoords) which make up each face. This data can be extracted #from s
