sources errors in advection/diffusion problems, and one solution

[James Pringle]

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 

Re: Solving advection, segregation and diffusion equation

[James Pringle]

Dear all -- in a bit I am going to post a little idealized
advection/diffusion problem to give an idea of how some of these errors in
steady solutions arise, and how iterating in time helps. It will be much
more basic then this work, but you might find it helpful. It was very
helpful for me in figuring out my more complex system to reproduce the
error in a simple way. It also reveals some surprising issues in the fipy
code -- for example, the error in a steady solution can depend on the
"initial" condition you feed the solver.

It will be under the title "sources errors in advection/diffusion problems,
and one solution." I hope to post it tonight.

Jamie

On Thu, Sep 15, 2016 at 1:27 PM, Zhekai Deng <
zhekaideng2...@u.northwestern.edu> wrote:

> Thanks for the reply. Those are very useful suggestions. I have some new
> findings which may help shedding some light on this problem.
>
> I have added TransientTerm and tried to see how the solution is developed
> over time. I also remove the datamin and datamax input on the viewer, so
> that I can see the full range of the solution variable data.
>
> When I see the temporal development of the solution (attached in the
> email), It looks stable. However, I noticed that the right outlet boundary
> has extreme high concentration and did not outflow the flow domain. I think
> It has due to the improper setup of my right outlet condition. Since the
> materials keep building up on the right outlet boundary, it appears to
> affect the solution in the actual flow domain(also attached in email) and
> drives the solution unstable. This could be the main reason that solving
> steady state equation directly never converges.
>
> I then begin to look into how to properly set up "open outlet" boundary,
> the closest I could find is following:
>
> https://www.mail-archive.com/fipy@nist.gov/msg02797.html
> 
> and in the section of "Applying fixed flux boundary conditions" of
> official website
> http://www.ctcms.nist.gov/fipy/documentation/USAGE.html
> 
>
> I tried to adapt the official method listed on the official website, but
> no success is achieved.  Thus, I wonder is it possible to point out some
> examples that used this method to define fixed flux outlet boundary
> condition ? You may find my current progress on this part in the code I
> attached.
>
> Another thing I noticed is that,  there is a reported issue on robin
> boundary condition:
>
> https://github.com/usnistgov/fipy/issues/426
> 
>
> Would it be possible to explain more about this ? I don't think the up &
> bottom boundary condition in my case causes the non convergent problem,
> because from the Transient equation, the segregation of the material starts
> to develop at the top boundary, which is consistent with my
> expectation.
>
> I have attached my current progress of the code, and the temporal
> development of the solution.
>
> Best,
>
> Zhekai
>
>
>
>
>
>
> On Thu, Sep 15, 2016 at 9:33 AM, Daniel Wheeler  > wrote:
>
>> On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
>>  wrote:
>> >
>> > I don't know offhand. With a Péclet number of 100, your problem is
>> almost completely hyperbolic, which FiPy (and cell-centered Finite Volume)
>> isn't very good at. Daniel knows more about this and may have some
>> suggestions.
>> >
>> > You might consider adding a TransientTerm for artificial relaxation and
>> trying the VanLeerConvectionTerm.
>>
>> Definitely use a TransientTerm and see the time step you can get away
>> with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm
>> isn't necessary as the accuracy is not important in a steady state
>> problem.
>>
>> --
>> Daniel Wheeler
>>
>> ___
>> fipy mailing list
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>> 
>>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
>> 

Re: Solving advection, segregation and diffusion equation

[Zhekai Deng]

Thanks for the reply. Those are very useful suggestions. I have some new
findings which may help shedding some light on this problem.

I have added TransientTerm and tried to see how the solution is developed
over time. I also remove the datamin and datamax input on the viewer, so
that I can see the full range of the solution variable data.

When I see the temporal development of the solution (attached in the
email), It looks stable. However, I noticed that the right outlet boundary
has extreme high concentration and did not outflow the flow domain. I think
It has due to the improper setup of my right outlet condition. Since the
materials keep building up on the right outlet boundary, it appears to
affect the solution in the actual flow domain(also attached in email) and
drives the solution unstable. This could be the main reason that solving
steady state equation directly never converges.

I then begin to look into how to properly set up "open outlet" boundary,
the closest I could find is following:

https://www.mail-archive.com/fipy@nist.gov/msg02797.html
and in the section of "Applying fixed flux boundary conditions" of official
website
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html

I tried to adapt the official method listed on the official website, but no
success is achieved.  Thus, I wonder is it possible to point out some
examples that used this method to define fixed flux outlet boundary
condition ? You may find my current progress on this part in the code I
attached.

Another thing I noticed is that,  there is a reported issue on robin
boundary condition:

https://github.com/usnistgov/fipy/issues/426

Would it be possible to explain more about this ? I don't think the up &
bottom boundary condition in my case causes the non convergent problem,
because from the Transient equation, the segregation of the material starts
to develop at the top boundary, which is consistent with my
expectation.

I have attached my current progress of the code, and the temporal
development of the solution.

Best,

Zhekai






On Thu, Sep 15, 2016 at 9:33 AM, Daniel Wheeler 
wrote:

> On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
>  wrote:
> >
> > I don't know offhand. With a Péclet number of 100, your problem is
> almost completely hyperbolic, which FiPy (and cell-centered Finite Volume)
> isn't very good at. Daniel knows more about this and may have some
> suggestions.
> >
> > You might consider adding a TransientTerm for artificial relaxation and
> trying the VanLeerConvectionTerm.
>
> Definitely use a TransientTerm and see the time step you can get away
> with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm
> isn't necessary as the accuracy is not important in a steady state
> problem.
>
> --
> Daniel Wheeler
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
from fipy import *

# mesh set up
nx = 100.
nz = 100.
dx = 1/nx
dz = 1/nz

mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)


# the follwing set up the velocity and shear rate in vector form.

# for those direction that I don't need it, I set its value as zero.

velocityX = mesh.faceCenters[1] # linear velocity profile
velocityY = FaceVariable(mesh=mesh, value = 0) # zero velocity in y direction
velocityVector = FaceVariable(mesh=mesh, rank=1)
velocityVector[0] = velocityX
velocityVector[1] = velocityY


local_shear_rate_X = FaceVariable(mesh=mesh, value = 0) # this term only acts 
on normal direction
local_shear_rate_Y = FaceVariable(mesh=mesh, value = 1) # linear velocity 
profile du/dz
local_shear_rate = FaceVariable(mesh=mesh, rank=1)
local_shear_rate[0] = local_shear_rate_X
local_shear_rate[1] = local_shear_rate_Y

Pe = 100. # some constants

#Setting up the diffusion part

D_Matrix = [[[0, 0],

 [0, 1/Pe]]]

#This is the new part to account for outflow fixed flux boundary (this does not 
set up correctly)
outflowflux = FaceVariable(mesh = mesh, value = phi.faceValue*velocityVector)
velocityVector.constrain(0.,mesh.facesRight)

# equation set up
eq = (TransientTerm() + UpwindConvectionTerm(coeff = velocityVector)\
+ UpwindConvectionTerm(coeff =local_shear_rate*(1-phi).faceValue)\
- (mesh.facesRight * outflowflux).divergence \
== DiffusionTerm (coeff=D_Matrix))


# set up the boundary conditions

phi.constrain(0.5, where = mesh.facesLeft) # fixed the inlet as 0.5 
concentration

phi.faceGrad.constrain(local_shear_rate*Pe*(1-phi.faceValue)*phi.faceValue, 
where = mesh.facesTop | mesh.facesBottom) #boundary conditions at the top and 
bottom


#phi.constrain(0, where = mesh.facesRight) #not working outlet boundary 
condition
#phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) # 
problematic outlet boundary 

JOR solver to help in accelerating sweeps?

[Gopalakrishnan, Krishnakumar]

Hello,

We have a specific question about how a solver helps in accelerating 
convergence of a loosely coupled system of PDEs.

Due to the fundamentally different properties/behaviour/geometry, we have a 
loosely coupled system of 8 PDEs defined along 5 different meshes, with data 
manually copied from certain regions of one mesh to suitable regions of another 
after sweeping each PDE.

Quite unsurprisingly, the system of equations needs a lot of sweeps (80 to 100) 
within each time-step to converge, and consequently the simulation model isn't 
useable. Although, the trilinos based parallel solvers are correctly set-up and 
working correctly for the example problems shipped with FiPy, the mpirun 
command complains when invoked on our (complicated script). We assume that it 
is due to the complicated setup of meshes and the classes/methods and objects 
that we have instantiated for our application. The serial code runs just fine.

So, in order to speed up our simulation, we looked to implement dynamic 
relaxation factors for the sweeps. This has only a limited success so far.

Next, I wonder whether including a Successive over-relaxation solver will help 
in this case ? There are 8 PDEs to be swept continuously until the residue of 
all of them gets lowered below a tolerance.  Does the Fipy's JOR solver 
(wrapper for Pysparse Jacobi/over-relaxation solvers)  help in this case ?  is 
there any other way to speed up our system ?


Best Regards

Krishna




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

Re: Solving advection, segregation and diffusion equation

[Daniel Wheeler]

On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
 wrote:
>
> I don't know offhand. With a Péclet number of 100, your problem is almost 
> completely hyperbolic, which FiPy (and cell-centered Finite Volume) isn't 
> very good at. Daniel knows more about this and may have some suggestions.
>
> You might consider adding a TransientTerm for artificial relaxation and 
> trying the VanLeerConvectionTerm.

Definitely use a TransientTerm and see the time step you can get away
with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm
isn't necessary as the accuracy is not important in a steady state
problem.

-- 
Daniel Wheeler

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