advection-diffusion example using fipy
Ondrej Certik
so after today's debacle with hermes2d, I decided to give finite
volumes a shot, so instead of implementing it myself, I used fipy.
Took me about 3 hours to learn, but it's nicely documented, so general
it's pretty cool.
The script is attached. There is a problem, that one has to choose an
approximation for the convection term, as described here:
http://www.ctcms.nist.gov/fipy/documentation/numerical/scheme.html#sec-numericalschemes
e.g. in particular, one can play with:
CentralDifferenceConvectionTerm, ExponentialConvectionTerm,
HybridConvectionTerm, PowerLawConvectionTerm, UpwindConvectionTerm,
ExplicitUpwindConvectionTerm, or VanLeerConvectionTerm,
so I tried all of them, and I got 3 different results and one fails completely:
# (1) looks the same, maybe it's ok
# (2) looks a bit worse than (1)
# (3) different, maybe also ok
convection = CentralDifferenceConvectionTerm # fails
convection = ExponentialConvectionTerm # (1)
convection = HybridConvectionTerm # (1)
convection = PowerLawConvectionTerm # (1)
convection = UpwindConvectionTerm # (2)
convection = ExplicitUpwindConvectionTerm # (3)
convection = VanLeerConvectionTerm # (3)
The plot for (1), (2) and (3) is attached. The domain is the same, BC is:
facesTopRight = ((mesh.getFacesRight() & (y > L / 2))
| (mesh.getFacesTop() & (x > L / 2)))
facesBottomLeft = ((mesh.getFacesLeft() & (y < L / 2))
| (mesh.getFacesBottom() & (x < L / 2)))
BCs = (FixedValue(faces=facesTopRight, value=0),
FixedValue(faces=facesBottomLeft, value=1))
e.g. I am assigning 1 to bottomleft, up to half of the edge, and
similarly 0 to the top right. I took this example from the fipy
distribution, just adapted it a bit. I am using:
epsilon = 0.01
b = (1, 1)
and I use the grid 40x40 quads.
So that's a disaster. But it seems to me that the exponential
convection scheme is the best, so I sticked to that and changed the BC
to:
facesTopRight = ((mesh.getFacesRight()) | (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft()) | (mesh.getFacesBottom()))
which is what we have in hermes2d. And now it seems to calculate some
sensible results. In particular, it can do the 300x300 mesh quite
easily. It is pure python and what is interesting is that it easily
beats hermes in terms of speed on the same mesh. Of course I know that
finite volumes are cheap to compute (less dofs), but still.
All of this convinces me that we should have a robust and fast finite
volume solver, based on the weak formulation. fipy doesn't use the
weak formulation (at least from the user perspective), also it doesn't
use the riemann solver (well, not the sophisticated one).
Ondrej
Pavel Solin
thanks for doing this.
Hi guys,
so after today's debacle with hermes2d, I decided to give finite
volumes a shot, so instead of implementing it myself, I used fipy.
Took me about 3 hours to learn, but it's nicely documented, so general
it's pretty cool.
The script is attached. There is a problem, that one has to choose an
approximation for the convection term, as described here:
http://www.ctcms.nist.gov/fipy/documentation/numerical/scheme.html#sec-numericalschemes
I looked here for the equation that is solved, but did not find it. Could you
send me a link to it?
e.g. in particular, one can play with:
CentralDifferenceConvectionTerm, ExponentialConvectionTerm,
HybridConvectionTerm, PowerLawConvectionTerm, UpwindConvectionTerm,
ExplicitUpwindConvectionTerm, or VanLeerConvectionTerm,
so I tried all of them, and I got 3 different results and one fails completely:
Interesting.
I like to compare results on the same level of accuracy. The FVM
results are likely to be less accurate than an equally-sized
FEM discretization. I suggest that we postpone the comparison
to a point where we have working finite volumes as well.
All of this convinces me that we should have a robust and fast finite
volume solver, based on the weak formulation. fipy doesn't use the
weak formulation (at least from the user perspective), also it doesn't
use the riemann solver (well, not the sophisticated one).
If you mean the Riemann solver for EUler equations, it is not needed
here. This is a linear problem, so one uses instead some of the simpler
schemes that you mentioned above.
In Hermes2D, the finite volume method should be a special case of DG
with piecewise-constant approximation that has one integration point at
each interior edge.
Pavel
Ondrej
--
You received this message because you are subscribed to the Google Groups "femhub" group.
To post to this group, send email to fem...@googlegroups.com.
To unsubscribe from this group, send email to femhub+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/femhub?hl=en.
--
Pavel Solin
University of Nevada, Reno
http://hpfem.math.unr.edu/people/pavel/
Hermes project: http://hpfem.org/
FEMhub project: http://femhub.org/
Daniel Wheeler
interesting. We certainly need a more robust higher order accurate
convection term discretization. It's not a show-stopped for our own
work so we haven't been motivated, but maybe we should. I have looked
at Hermes with thoughts of incorporating some of the discretization
schemes, but haven't made much headway. Is there a way that we could
get Hermes 2D and fipy playing together? I'm thinking that fipy is a
relatively well documented simple front end and Hermes has better
schemes etc. Any thoughts or are the two approaches too far apart?
Cheers.
--
Daniel Wheeler
Pavel Solin
Hi, I work on the FiPy project and the analysis below is very
interesting. We certainly need a more robust higher order accurate
convection term discretization. It's not a show-stopped for our own
work so we haven't been motivated, but maybe we should.
It seems that the general trend is moving towards higher-order methods, so
probably this would not hurt.
I have looked
at Hermes with thoughts of incorporating some of the discretization
schemes, but haven't made much headway. Is there a way that we could
get Hermes 2D and fipy playing together?
That would be awesome. I think that we should start by learning more
about each other's code. This definitely holds for myself, I need to
learn more about FiPy.
I'm thinking that fipy is a
relatively well documented simple front end and Hermes has better
schemes etc. Any thoughts or are the two approaches too far apart?
Could you describe in a few words what your ambition with FiPy is --
I mean, what do you want FiPy to be known for, what are the equations
you want to solve, what functionality you would really like to do yourselves?
For Hermes some of this is described here:
http://hpfem.org/hermes2d/doc/src/introduction.html#what-is-hermes.
BTW, some time ago we included FiPy in FEMhub. The goals of FEMhub
are described here: http://femhub.org/. Do you guys also care for elliptic
problems, or you solve mostly transport? The most we can do ourselves
given our limited manpower is to wrap codes in Python (if they need that)
and add them to FEMhub, but we need help with unifying interfaces. Would
you be willing to modify existing FiPy interface or create a new simple
interface for FEMhub?
Best,
Pavel
P.S. The linear advection-diffusion equation that Ondrej
mentioned in his first email works now fine, we just forgot to use
integration weights in the definition of the weak forms. The example
is now in the git repo
http://hpfem.org/git/gitweb.cgi/hermes2d.git/tree/HEAD:/examples/linear-advection-diffusion
and I will create a Sphinx description today:
Ondrej Certik
On Mon, Mar 29, 2010 at 12:18 PM, Pavel Solin <so...@unr.edu> wrote:
> Hello Daniel,
[...]
>
> P.S. The linear advection-diffusion equation that Ondrej
> mentioned in his first email works now fine, we just forgot to use
> integration weights in the definition of the weak forms. The example
> is now in the git repo
> http://hpfem.org/git/gitweb.cgi/hermes2d.git/tree/HEAD:/examples/linear-advection-diffusion
> and I will create a Sphinx description today:
yes, it was a stupid mistake on our side, now that problem seems to be
working fine (including hp adaptivity).
As to fipy, here is the spkg package:
http://femhub.googlecode.com/files/fipy-2.0.2.spkg
that you can install into femhub by "femhub -i
http://femhub.googlecode.com/files/fipy-2.0.2.spkg", I need to check
that it builds on all platforms that we support (since it's pure
python, there should be no problems, but it might require some new
setuptools) and if all is fine, we'll just put it into our official
femhub.
As Pavel said, it'd be cool to integrate fipy well in femhub, so that
you can combine fipy and hermes, but we need some help with that.
Ondrej
Daniel Wheeler
> That would be awesome. I think that we should start by learning more
> about each other's code. This definitely holds for myself, I need to
> learn more about FiPy.
Sound reasonable. I'll try and carve out some time for Hermes and try
and solve a phase field problem, something that is important to our
group.
> Could you describe in a few words what your ambition with FiPy is --
> I mean, what do you want FiPy to be known for, what are the equations
> you want to solve, what functionality you would really like to do
> yourselves?
We have no overarching ambition for the code. FiPy was simply
implemented to answer our needs. We wanted to work in python and there
wasn't an easy to use code available that used the latest tools that
we could easily pass to others. The first priority is that it allows
us to do our work and publish results, the second is that it is
packaged in such a way that we can share our results, which is part of
NIST's mission.
Our primary interest id phase field problems and coupled PDEs that are
typical to phase file models. For example, I've just been working on a
"reactive wetting" problem that involves solving modified compressible
Navier-Stokes, a phase field equation and binary mass transport. Jon
is working on solving drift-diffusion equations for photo-voltaics,
which isn't a phase field method.
> For Hermes some of this is described here:
> http://hpfem.org/hermes2d/doc/src/introduction.html#what-is-hermes.
>
> BTW, some time ago we included FiPy in FEMhub. The goals of FEMhub
> are described here: http://femhub.org/. Do you guys also care for elliptic
> problems, or you solve mostly transport?
Some elliptic and some transport. The Allen-Cahn equation (phase field
equation) is an elliptic PDE.
> The most we can do ourselves
> given our limited manpower is to wrap codes in Python (if they need that)
> and add them to FEMhub, but we need help with unifying interfaces. Would
> you be willing to modify existing FiPy interface or create a new simple
> interface for FEMhub?
Yes, but this probably needs to work in two ways. Firstly, FiPy needs
to be accessible to FEMhub. Secondly, FiPy needs to have access to
better discretization routines.
The main problem is carving out time to take care of these things.
Probably trying to solve a phase field problem using Hermes may be the
best place for me to start and using FiPy in FEMhub to check how it
works.
--
Daniel Wheeler
Daniel Wheeler
> As to fipy, here is the spkg package:
>
> http://femhub.googlecode.com/files/fipy-2.0.2.spkg
>
> that you can install into femhub by "femhub -i
> http://femhub.googlecode.com/files/fipy-2.0.2.spkg", I need to check
> that it builds on all platforms that we support (since it's pure
> python, there should be no problems, but it might require some new
> setuptools) and if all is fine, we'll just put it into our official
> femhub.
>
> As Pavel said, it'd be cool to integrate fipy well in femhub, so that
> you can combine fipy and hermes, but we need some help with that.
I'm looking forward to getting my hands dirty with this some time soon. Cheers.
--
Daniel Wheeler
Pavel Solin
I'd be happy to help with the weak forms if you send me the equations.
Pavel
and using FiPy in FEMhub to check how it
works.
--
Daniel Wheeler
--
You received this message because you are subscribed to the Google Groups "femhub" group.
To post to this group, send email to fem...@googlegroups.com.
To unsubscribe from this group, send email to femhub+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/femhub?hl=en.