Re: sources errors in advection/diffusion problems, and one solution

[Zhekai Deng]

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
>> >
>> >

Re: Dynamic under-relaxation factors for FiPy sweep

[Guyer, Jonathan E. Dr. (Fed)]

There isn't enough here to tell what's going on. Can you show us the code?

> On Sep 19, 2016, at 6:53 PM, Gopalakrishnan, Krishnakumar 
>  wrote:
> 
> Hi Dan,
> 
> Thanks for your suggestion to use the ResidualTerm as per your gist posting 
> https://gist.github.com/guyer/f29c759fd7f0f01363b8483c7bc644cb  of the 
> Newton's method.
> 
> When I try to implement the Newton's method into my code, python interpreter 
> gives the following error message
> 
> TypeError: __init__() got an unexpected keyword argument 'var' 
> File " \fipy-3.1-py2.7.egg\fipy\terms\term.py", line 428, in __eq__
>return self - other
>  File "\fipy-3.1-py2.7.egg\fipy\terms\term.py", line 422, in __sub__
>return self + (-other)
>  File "\fipy-3.1-py2.7.egg\fipy\terms\abstractBinaryTerm.py", line 88, in 
> __neg__
>return (-self.term) + (-self.other)
>  File "\fipy-3.1-py2.7.egg\fipy\terms\nonDiffusionTerm.py", line 56, in 
> __neg__
>return self.__class__(coeff=-self.coeff, var=self.var)
> TypeError: __init__() got an unexpected keyword argument 'var'
> 
> I am following exactly the same steps given in the original gist posting. Any 
> idea what might be wrong here ?
> 
> 
> Krishna 
> 
> 
> 
> -Original Message-
> From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of 
> Daniel Wheeler
> Sent: Tuesday, September 6, 2016 1:46 PM
> To: Multiple recipients of list 
> Subject: Re: Dynamic under-relaxation factors for FiPy sweep
> 
> On Tue, Sep 6, 2016 at 6:58 AM, Krishna  wrote:
>> 
>> Since python to be a very distributed ecosystem, this question for 
>> some kind of a starter code,  may not fit well in  a 
>> general/computational math stackexchange post , nor in this mailing 
>> list. fipy's details are certainly required to implement an Aitken 
>> type dynamic under relaxation. , I.e. one needs access to the internal 
>> and residual matrices, in order to apply text book formulae, and then 
>> split the relaxation vectors into individual scalars for use in the 
>> 'underrelaxation' parameter for each sweep method. The first two 
>> sweeps must be static/initial 'underrelaxation' so that we can apply the 
>> formula.
> 
> I see. Here is an example of doing Newton iterations in FiPy
> 
>https://gist.github.com/guyer/f29c759fd7f0f01363b8483c7bc644cb
> 
> It uses the ResidualTerm. If you look at that code, it uses the 
> justResidualVector, which gives the residual vector. You can also get access 
> to the matrix and b vector separately. For the under relaxation, I don't 
> think it's possible to apply it as a vector that's different for each 
> equation. There is probably some way to do it akin to what's happening in the 
> ResidualTerm.
> 
> 
> 
> --
> Daniel Wheeler
> ___
> 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 ]


___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: sources errors in advection/diffusion problems, and one solution

[James Pringle]

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
>>> >