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 <jprin...@unh.edu> 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 <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 >> > >> > <JMP_Make_Grids.py><SquareGrid_HomemadeDelaunay_transient. >> py><SquareGrid_HomemadeDelaunay.py>_________________________ >> ______________________ >> > fipy mailing list >> > fipy@nist.gov >> > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcm >> s.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3 >> EH6RqKWf16fbYIYKw&m=IaTL4nEetwjm4G4qfj8cwLzvztKzuiFsw_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=7JGh0Zz82O69 >> cJiLFKmkV3NfW2TLz6KB_ngkAGCrYGI&e= ] >> >> >> _______________________________________________ >> fipy mailing list >> fipy@nist.gov >> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcm >> s.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3 >> EH6RqKWf16fbYIYKw&m=IaTL4nEetwjm4G4qfj8cwLzvztKzuiFsw_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=7JGh0Zz82O69 >> cJiLFKmkV3NfW2TLz6KB_ngkAGCrYGI&e= ] >> > > > _______________________________________________ > 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 ]
