using sympy for representing the weakforms
Ondrej Certik
so we had a bonfire with Pavel today and as it usually goes, we made
some plans, so I am writing here to get some feedback. :)
We were discussing the possibility to define the weakforms in sympy
and we found that we need 4 things.
We'll start with the weak formulation (later we may even get some
automatic way to get from PDE to the weak formulation, but that's not
our current problem, as this could be done easily by the user). So
here are the steps:
1) The weak formulation (in 3D) contains volume integrals, boundary
(surface) integrals, all the weakforms will be represented as sympy
expressions
2) we prescribe a mesh, and at all boundaries we prescribe if we have
essential (we remove the dof from the system and prescribe a value for
it using the Dirichlet lift, e.g. this is the Dirichlet B.C.) or
natural (one calculates the boundary integral, and one can either
prescribe a value of the integrand, or leave it free) boundary
conditions. There needs to be a way to represent which boundary
integral belongs to which part of the mesh.
3) we prescribe the definitions of all material parameters in the
forms (e.g. either a number, or some function, either analytic or
using a table, or a previous solution), so in the Schroedinger
equation this means the potential, in the heat equation this means the
conductivity and sources.
4) in the last step one tells hermes what to do with the forms. E.g.
if it should construct matrices out of it and solve the eigenproblem,
or it should use a Newton method to solve the nonlinear problem etc.
The main idea of this approach is that sympy is used to represent the
equations and then the solver (be it hermes or libmesh or anything
else) uses the sympy api to poll any information it needs from it to
be able to solve it. And if it can't solve it -- let's say the user
prescribes some surface integral that the solver can't do yet -- it
will just raise an exception.
For example Robert in sfepy is using the input files to represent the
problem. I think if this is done in SymPy, one can then play with it,
for example one can even calculate things symbolically for small
number of base functions and one can even calculate things with non
local base functions, which can be useful in many cases (e.g. solving
the Schroedinger equation using atomic orbitals basis), not mentioning
that one can print the weak formulation in latex form, or convert to
Sage or to anything else that sympy can handle. Also it can be used to
automatically derive the Newton iterations, which is very tedious by
hand.
I'll start using this approach in hermes1D, this will be easy to
represent in sympy and also quite easy to solve. Then later as time
permits I'll go to 2D and 3D and to more complicated equations.
Ondrej
Ondrej Certik
> We were discussing the possibility to define the weakforms in sympy
> and we found that we need 4 things.
I have created an issue for it:
http://code.google.com/p/femhub/issues/detail?id=26
and quickly coded an example for a quantum mechanic harmonic
oscillator in 1D, you can browse both the sympy code and the pretty
printed output of it. Then I'll hook it into hermes1d and it will
serve as the standard (highest level) way of defining things.. It's
nice I like this approach a lot.
Ondrej
Jakub Cerveny
integration would be done. The integration loops simply cannot be done
in Python, which is too slow for this. I am not against doing some
high-level stuff in Python but the innermost cycles where most CPU
time is spent must be done in C.
Jakub
Ondrej Certik
>
> I don't see how the weak forms would be evaluated, i.e., how the
> integration would be done. The integration loops simply cannot be done
> in Python, which is too slow for this. I am not against doing some
> high-level stuff in Python but the innermost cycles where most CPU
> time is spent must be done in C.
Yes, this will be done in C. On the hermes part, almost nothing
changes, besides, the user is not required to use Python at all, if he
doesn't want to. As we discussed previously, ideally there should be
many levels on which the user can use hermes --- from very low level
one, up to a high level (Python) one.
Ondrej
David Andrs
Currently hermes1d is in python, so you can use this quite easily, I think.
But, 2D and 3D code are already in C++. If there is a predefined form in
hermes you just need to do some "match" to figure out what to call and
then call it. That is clear and I guess it is not a big deal in Python.
But what are you going to do if there is not a predefined form in
hermes? Or if you want to use some user-defined function in the
integral? Are you going to generate some Cython code along with
makefiles? Or something else?
I think this is what Jakub was asking about...
And let me also share my idea of this high level "layer". For now it is
just a concept. I'd like to have a meta language to describe the problem
and then use some stand-alone tool to prepare the front-end for hermes.
The output of such a tool would be Python/C++/whatever. And it would
also generate all necessary files. So in case of C++, it would be the
file with C++ code, but also a cmake script, or a Makefile, or whatever
would be needed.
--
David
Ondrej Certik
>
> I guess Jakub was asking *how* you are going to connect Python and C.
> Currently hermes1d is in python, so you can use this quite easily, I think.
Yes, but it will be written in C++ soon, so it's the same problem.
>
> But, 2D and 3D code are already in C++. If there is a predefined form in
> hermes you just need to do some "match" to figure out what to call and
> then call it. That is clear and I guess it is not a big deal in Python.
> But what are you going to do if there is not a predefined form in
> hermes? Or if you want to use some user-defined function in the
> integral? Are you going to generate some Cython code along with
> makefiles? Or something else?
Here is one way to do it:
In the first iteration, I'll create a special hermes2d form (in C++),
let's call it "python_form", that allows the user defined data to be
passed in. Then one will be able to control such data from Python.
This "python_form" is just C++ code, and as to hermes, it doesn't know
anything about Python.
I think Python is really confusing you guys. :) Maybe because you
think that if something is usable from Python, it has to be slow. It
is not true -- if done correctly, it will be as fast as in C++. If it
helps, forget about Python, and just imagine that you would like to
specify the forms on the commandline. E.g. you create a hermes2d
example and you will specify which forms to solve using options at the
runtime.
The crucial idea is that you don't have to compile anything. This is
really important.
In the second iteration, one can improve the above, we'll see. Maybe
there are also other ways too.
>
> I think this is what Jakub was asking about...
>
> And let me also share my idea of this high level "layer". For now it is
> just a concept. I'd like to have a meta language to describe the problem
> and then use some stand-alone tool to prepare the front-end for hermes.
> The output of such a tool would be Python/C++/whatever. And it would
> also generate all necessary files. So in case of C++, it would be the
> file with C++ code, but also a cmake script, or a Makefile, or whatever
> would be needed.
The big disadvantage of this approach is that you still have to
compile things whenever you want to solve something. I would like to
find a way so that this is not needed, because users will just install
a binary and that's it. Or use it over the notebook -- and while it is
feasible to compile stuff, I'd like to avoid it. It will be faster
too.
Ondrej
Pavel Solin
being able to define our weak formulations in such a way that they
can be further processed on symbolic level (I mean mainly symbolic
differentiation in order to construct the Jacobian matrix) has a huge
potential in it. We can make the Newton's method accessible to
everyone. Right now, most people are scared of differentiating the
weak formulation on paper. (And I understand why -- it took me seven
full pages to differentiate 2D compressible Euler equations and I the
results may be flawed. Honestly, I would not like to do Navier-Stokes
equations or MHD equations in 3D).
Jakub, are the following files in the src/ directory currently used?
mesh_lexer.cpp mesh_lexer.h mesh_lexer.l weakform_lexer.h
weakform_lexer.l
Best,
Pavel
--
Pavel Solin
University of Nevada, Reno
http://hpfem.math.unr.edu/people/pavel
Visit Hermes project at http://hpfem.org
Ondrej Certik
I didn't answer all your questions --- as to user defined functions,
in the first iteration, I will not allow it.
In the second iteration, I'll use the same approach. Many user defined
functions (e.g. material parameters) are given using a table, so I'll
create a C++ function that can handle any such functions.
In the third iteration, I'll add more functionality. And so on.
In the tenth iteration, one should be able to solve most of the common
equations, with varying parameters, etc. However, there will always be
cases, which it cannot handle -- then the user will just use a little
lower level of hermes API, and just write C++ code, the same way as
now. But all cases, that I have been using hermes2d for so far, I
think those can all be handled by the way I described.
Ondrej
David Andrs
example python that uses sympy notation/calls and the standalone tool is
some script, that can process that (so obviously also Python). But, with
such a tool you can do more than just generate the forms...
You would be able to generate for example tex commands (to include this
in tex documents) or html for web pages or whatever...
And yes, the purpose of such a tool is exactly the same what you are
writing about: to avoid 7 pages long derivation where you can easily
make a mistake...
--
David
Jakub Cerveny
> mesh_lexer.cpp mesh_lexer.h mesh_lexer.l
Yes, without them you could not load mesh files.
> weakform_lexer.h
> weakform_lexer.l, weakform_parser.cpp
These are not used.
This was supposed to be the weak form expression parser, but I never
got to implementing it. It would BTW solve your problem with the need
to compile stuff when the form changes.
Jakub
Robert Cimrman
I fear you are reinventing a wheel here a bit, checkout [1]. The guys at
the Simula lab have this for several years already (and yes, it's usable
from/written in Python). You seem to want to recreate their UFC
(universal form compiler) stuff but using sympy instead of ginac...
I suggest you read all the docs they have already published on this
topic before venturing further, and if possible, (re)use their ideas and
code. They are smart guys and have several (many) years of experience
with symbolic FE. They have some adaptivity too in the Fenics project.
Otherwise +1, I am tired of writing the C forms by hand, so count with
me (at least as the devil's advocate ;)). Also, as a start, we could
revive SymFE [*]? Or do you plan to develop within SymPy itself from the
very beginning?
r.
[1] http://www.fenics.org/wiki/SyFi
[*] SymFE (http://code.google.com/p/symfe/) was started some (quite a
long) time ago exactly for that purpose, but was stoved aside by me for
other more pressing things. The SyFi guys were involved too.
Ondrej Certik
>
> Hi guys,
>
> I fear you are reinventing a wheel here a bit, checkout [1]. The guys at
> the Simula lab have this for several years already (and yes, it's usable
> from/written in Python). You seem to want to recreate their UFC
> (universal form compiler) stuff but using sympy instead of ginac...
>
> I suggest you read all the docs they have already published on this
> topic before venturing further, and if possible, (re)use their ideas and
> code. They are smart guys and have several (many) years of experience
> with symbolic FE. They have some adaptivity too in the Fenics project.
Have you seen some more advanced problem solved by syfi? I can't find
any examples.
I think they use symbolics too much, but I may be wrong. I just want
to use sympy to represent the equations, everything else should be
done by hermes.
But otherwise you are right, we are definitely not alone with this
idea and we have to work very hard on femhub and all of this, so that
we can get wider acceptance.
>
> Otherwise +1, I am tired of writing the C forms by hand, so count with
> me (at least as the devil's advocate ;)). Also, as a start, we could
> revive SymFE [*]? Or do you plan to develop within SymPy itself from the
> very beginning?
Let's use symfe, no problem.
Ondrej
Ondrej Certik
> 2009/5/25 Robert Cimrman <cimr...@ntc.zcu.cz>:
>>
>> Hi guys,
>>
>> I fear you are reinventing a wheel here a bit, checkout [1]. The guys at
>> the Simula lab have this for several years already (and yes, it's usable
>> from/written in Python). You seem to want to recreate their UFC
>> (universal form compiler) stuff but using sympy instead of ginac...
>>
>> I suggest you read all the docs they have already published on this
>> topic before venturing further, and if possible, (re)use their ideas and
>> code. They are smart guys and have several (many) years of experience
>> with symbolic FE. They have some adaptivity too in the Fenics project.
>
> Have you seen some more advanced problem solved by syfi? I can't find
> any examples.
>
> I think they use symbolics too much, but I may be wrong. I just want
> to use sympy to represent the equations, everything else should be
> done by hermes.
Here is the manual:
http://www.fenics.org/pub/documents/syfi/syfi-user-manual/syfi-user-manual.pdf
and indeed, it contains many examples. Yes, that's pretty good. If it
works for all kinds of domains and boundary conditions, then I like
it.
Ondrej
Robert Cimrman
> 2009/5/25 Robert Cimrman <cimr...@ntc.zcu.cz>:
>> Hi guys,
>>
>> I fear you are reinventing a wheel here a bit, checkout [1]. The guys at
>> the Simula lab have this for several years already (and yes, it's usable
>> from/written in Python). You seem to want to recreate their UFC
>> (universal form compiler) stuff but using sympy instead of ginac...
>>
>> I suggest you read all the docs they have already published on this
>> topic before venturing further, and if possible, (re)use their ideas and
>> code. They are smart guys and have several (many) years of experience
>> with symbolic FE. They have some adaptivity too in the Fenics project.
>
> Have you seen some more advanced problem solved by syfi? I can't find
> any examples.
I follow the commit announcement list - they solve navier stokes,
hyperelasticity (large deformation), basically everything one needs. I
have not tried the code myself though. Anyway, I think it is worth
looking, when we set out on this path!
> I think they use symbolics too much, but I may be wrong. I just want
> to use sympy to represent the equations, everything else should be
> done by hermes.
Yes they start from the base functions, which leads to a very
complicated expressions later IMHO. I am thinking too about just
representing the equations, with some not so low-level mapping to the
"computing library". So one of the things we should look at is the
language/syntax they use to describe the forms [1].
[1] http://www.fenics.org/wiki/Documentation
As for the Newton, I have always thought that the automatic
differentiation stuff (AD) [2] is just the same as deriving the final
expression symbolically but at a compiler level. It is not so, instead
the gradient evaluation in AD is just as complex as the function itself
(I am over-simplifying here...) [3]. This is IMHO a topic worth looking too.
[2] http://en.wikipedia.org/wiki/Automatic_differentiation
[3]
http://justindomke.wordpress.com/2009/03/24/a-simple-explanation-of-reverse-mode-automatic-differentiation/
> But otherwise you are right, we are definitely not alone with this
> idea and we have to work very hard on femhub and all of this, so that
> we can get wider acceptance.
>
>> Otherwise +1, I am tired of writing the C forms by hand, so count with
>> me (at least as the devil's advocate ;)). Also, as a start, we could
>> revive SymFE [*]? Or do you plan to develop within SymPy itself from the
>> very beginning?
>
> Let's use symfe, no problem.
ok.
I have raised the issues above because those are very interesting topics
that I wanted to pursue (learning about fenics, AD), but never found
power to do it actually. You guys seem to have the drive to realize it :)
r.
Robert Cimrman
I like the input files. Playing with PDEs on a command line, or in a
one-shot script is nice, but once the equations settle, I just want to
describe the problem and solve it, I do not want to "program" again. So
I am saying that input files are generally useful, as an alternative.
In sfepy, we solve quite complex equations with different behaviour on
subdomains, various boundary conditions, material coefficients depending
on space/time, basically all you mentioned above, and our problem
description syntax handles it well (for us). Recently we have also made
some steps to make it less verbose, see [1]. Therefore let me offer you
this format or better "agreement", as a problem description file is just
a plain Python module with some "keywords" like "equations", "regions",
etc. SymPy could be easily made to read it and do the symbolic stuff
just as you describe.
I volunteer to tweak it so that it covers hp-fem related stuff too, in
case you are +1.
Pavel Solin
Ondrej Certik
> I like the input files. Playing with PDEs on a command line, or in a
> one-shot script is nice, but once the equations settle, I just want to
> describe the problem and solve it, I do not want to "program" again. So
Right. But you still want to tweak it a bit, like changing a mesh,
B.C. conditions, etc? If not, then I can also create a simple Python
function that does that.
> I am saying that input files are generally useful, as an alternative.
>
> In sfepy, we solve quite complex equations with different behaviour on
> subdomains, various boundary conditions, material coefficients depending
> on space/time, basically all you mentioned above, and our problem
> description syntax handles it well (for us). Recently we have also made
> some steps to make it less verbose, see [1]. Therefore let me offer you
> this format or better "agreement", as a problem description file is just
> a plain Python module with some "keywords" like "equations", "regions",
> etc. SymPy could be easily made to read it and do the symbolic stuff
> just as you describe.
I think those things from the input file can be easily passed as
arguments to some function. E.g. the solver pararmeters. The equation
itself will be represented in sympy and we can hook this into sfepy
too.
In the file you sent, you seem to be mixing both the weak formulation
and the type of elements and maybe other things. I think a cleaner
approach is what I described in my first email --- to just represent
the weak formulation somehow (sympy) and then specify other things,
like the spaces.
>
> I volunteer to tweak it so that it covers hp-fem related stuff too, in
> case you are +1.
As to SiFy, we have to check it out. Please talk about this with Pavel
in Roznov, or Plzen.
From what you said I think what they do is very similar to what we
want. They get very complicated expressions in the end, as you pointed
out, which honestly I don't like much, I'd like to know the pros and
cons of that approach. I prefer the way we seem to agree on, to just
use sympy to represent the equations, but otherwise use the regular
hp-FEM as we have now.
Ondrej
Pavel Solin
I hope that we are not losing the focus. It is not only representation, but
also (mainly) symbolic differentiation of the weak formulation for the
Newton's method. Does Fenics do this?
Pavel
Ondrej
Robert Cimrman
> Hi Robert,
>
>> I like the input files. Playing with PDEs on a command line, or in a
>> one-shot script is nice, but once the equations settle, I just want to
>> describe the problem and solve it, I do not want to "program" again. So
>
> Right. But you still want to tweak it a bit, like changing a mesh,
> B.C. conditions, etc? If not, then I can also create a simple Python
> function that does that.
Yes, I want this too. In sfepy, that's the purpose of a parametric hook
function.
>> I am saying that input files are generally useful, as an alternative.
>>
>> In sfepy, we solve quite complex equations with different behaviour on
>> subdomains, various boundary conditions, material coefficients depending
>> on space/time, basically all you mentioned above, and our problem
>> description syntax handles it well (for us). Recently we have also made
>> some steps to make it less verbose, see [1]. Therefore let me offer you
>> this format or better "agreement", as a problem description file is just
>> a plain Python module with some "keywords" like "equations", "regions",
>> etc. SymPy could be easily made to read it and do the symbolic stuff
>> just as you describe.
>
> I think those things from the input file can be easily passed as
> arguments to some function. E.g. the solver pararmeters. The equation
> itself will be represented in sympy and we can hook this into sfepy
> too.
Exactly.
> In the file you sent, you seem to be mixing both the weak formulation
> and the type of elements and maybe other things. I think a cleaner
> approach is what I described in my first email --- to just represent
> the weak formulation somehow (sympy) and then specify other things,
> like the spaces.
What do you mean by mixing?
In the file, "equations" are the weak formulation (but yes, the
integrals are specified there - there could be "auto" for hp-fem), and
"fields" are like the spaces in hermes.
>> I volunteer to tweak it so that it covers hp-fem related stuff too, in
>> case you are +1.
>
> As to SiFy, we have to check it out. Please talk about this with Pavel
> in Roznov, or Plzen.
We have lots of things to talk about, definitely.
r.
Robert Cimrman
>> As to SiFy, we have to check it out. Please talk about this with Pavel
>> in Roznov, or Plzen.
>>
>> From what you said I think what they do is very similar to what we
>> want. They get very complicated expressions in the end, as you pointed
>> out, which honestly I don't like much, I'd like to know the pros and
>> cons of that approach. I prefer the way we seem to agree on, to just
>> use sympy to represent the equations, but otherwise use the regular
>> hp-FEM as we have now.
>
>
> I hope that we are not losing the focus. It is not only representation, but
> also (mainly) symbolic differentiation of the weak formulation for the
> Newton's method. Does Fenics do this?
I think so. They have the weak form expressed symbolically.
r.
Pavel Solin
>> As to SiFy, we have to check it out.
+1
I looked briefly at this example at SyFi home page (below).
>The following example shows the computation of the element matrix for the Poisson problem in Python,
>triangle = ReferenceTriangle()This looks short and elegant but it only works for triangles with straight edges and constant material
> fe = LagrangeFE(triangle, 3)
>A = {}
>for i in range(0, fe.nbf()):
> for j in range(0, fe.nbf()):
> integrand = inner(grad(fe.N(i)), grad(fe.N(j)))
> A[(i,j)] = triangle.integrate(integrand)
properties. From its beginnings, Fenics was based on simplifications like this. Has something changed?
I would be very interested in seeing whether and how they handle curvilinear elements, nonconstant
PDE coefficients, etc. I have not found anything yet regarding symbolic expressions for the weak
formulation. If you mean this,
integrand = inner(grad(fe.N(i)), grad(fe.N(j)))it is not enough.
Pavel
Pavel Solin
I just wonder since for these elements it is not enough to
precompute a single local stiffness matrix on the reference
element and use it for all elements in the mesh. This was the
original concept of Fenics - they spent lots of time optimizing
some higher level stuff that was built on this naive foundation.
I am very curious whether they still stick to this original concept.
Pavel
Robert Cimrman
> Has Fenics nonsimplicial elements (quads, bricks, pyramids)?
> I just wonder since for these elements it is not enough to
> precompute a single local stiffness matrix on the reference
> element and use it for all elements in the mesh. This was the
> original concept of Fenics - they spent lots of time optimizing
> some higher level stuff that was built on this naive foundation.
> I am very curious whether they still stick to this original concept.
I cannot tell you such details as I do not know them :-) When I last
looked, the foundations was (IMHO) just as you say. I just wanted you to
look at Fenics, as it seems similar in goals, so there _might_ be ideas
useful for us - it is the only code that employs symbolic computations
for FE.
As for the matrices, I do not think they have to be precomputed, Fenics
can (IMHO) translate that to an optimized C(++?) code. The symbolic
stuff serves as a foactory for making dynamic C libraries of optimized
element matrices on the fly. But I can be wrong...
r.
Brian Granger
few thoughts...
... about specifiying weak forms symbolically using sympy:
* It should all be pure python. I don't like the idea of having to
deal with YAIFF (Yet Another Input File Format).
* The symbolics should be independent of the finite element methods.
Why? I would like to use the same symbolic expressions in my B-spline
solvers. Also, we want the symbolic API to be used by people who know
0.0 about the FEM.
My $0.02
Cheers,
Brian
Robert Cimrman
> I haven't been following this discussion to closely, but I do have a
> few thoughts...
>
> ... about specifiying weak forms symbolically using sympy:
>
> * It should all be pure python. I don't like the idea of having to
> deal with YAIFF (Yet Another Input File Format).
+1. Sfepy input files are just plain python modules. But e.g. the
equations are given in a string, and this is where some agreement have
to be reached (call it format or not), that this symbol denotes the
Laplacian operator, that denotes the linear elasticity term etc.
> * The symbolics should be independent of the finite element methods.
> Why? I would like to use the same symbolic expressions in my B-spline
> solvers. Also, we want the symbolic API to be used by people who know
> 0.0 about the FEM.
+1. FEM is a detail of a particular discretization.
cheers,
r.
Pavel Solin
+1. Sfepy input files are just plain python modules. But e.g. the
Brian Granger wrote:
> I haven't been following this discussion to closely, but I do have a
> few thoughts...
>
> ... about specifiying weak forms symbolically using sympy:
>
> * It should all be pure python. I don't like the idea of having to
> deal with YAIFF (Yet Another Input File Format).
equations are given in a string, and this is where some agreement have
to be reached (call it format or not), that this symbol denotes the
Laplacian operator, that denotes the linear elasticity term etc.
I think that the input format of the equations should be just "math".
This is not YAIFF and it is computer language independent.
I thought that SymPy would be a good tool to do this. If there is
a Python parser for the math expressions, this is fine with me.
Are we in disagreement here?
> * The symbolics should be independent of the finite element methods.
> Why? I would like to use the same symbolic expressions in my B-spline
> solvers. Also, we want the symbolic API to be used by people who know
> 0.0 about the FEM.
The only connection to FEM is that we need to define the weak
formulation (NOT a plain PDE). This is because several different
weak formulations can usually be obtained from a PDE,
and second, that the weak formulation has everything incorporated
in it (boundary conditions etc.), while they have to be presented
separately with a standard PDE. Third, the weak formulation
is the basis for FEM...
Pavel
+1. FEM is a detail of a particular discretization.
cheers,
r.
Ondrej Certik
Brian, is your splines code also working with the weak formulation?
Maybe it's not just called like that, but it seems to me that whenever
I discretize some operator in a basis, I need to multiply by basis
functions and integrate, which is a discretized weak formulation.
Ondrej
Ondrej Certik
> On Fri, May 29, 2009 at 4:57 AM, Pavel Solin <so...@unr.edu> wrote:
>>
>>
>> On Wed, May 27, 2009 at 11:35 PM, Robert Cimrman <cimr...@ntc.zcu.cz>
>> wrote:
>>>
>>> Brian Granger wrote:
>>> > I haven't been following this discussion to closely, but I do have a
>>> > few thoughts...
>>> >
>>> > ... about specifiying weak forms symbolically using sympy:
>>> >
>>> > * It should all be pure python. I don't like the idea of having to
>>> > deal with YAIFF (Yet Another Input File Format).
>>>
>>> +1. Sfepy input files are just plain python modules. But e.g. the
>>> equations are given in a string, and this is where some agreement have
>>> to be reached (call it format or not), that this symbol denotes the
>>> Laplacian operator, that denotes the linear elasticity term etc.
>>
>> I think that the input format of the equations should be just "math".
>> This is not YAIFF and it is computer language independent.
>> I thought that SymPy would be a good tool to do this. If there is
>> a Python parser for the math expressions, this is fine with me.
>> Are we in disagreement here?
I think we are not. Only that I don't want to write yet another parser
for mathematical expressions, since I already wrote one (sympy), and
there are lots of corner cases.
Ondrej
Brian Granger
>> equations are given in a string, and this is where some agreement have
>> to be reached (call it format or not), that this symbol denotes the
>> Laplacian operator, that denotes the linear elasticity term etc.
>
> I think that the input format of the equations should be just "math".
> This is not YAIFF and it is computer language independent.
> I thought that SymPy would be a good tool to do this. If there is
> a Python parser for the math expressions, this is fine with me.
> Are we in disagreement here?
Nope, I am in full agreement and I should have clarified this. I
think using sympy for this is a fantastic idea!
>> > * The symbolics should be independent of the finite element methods.
>> > Why? I would like to use the same symbolic expressions in my B-spline
>> > solvers. Also, we want the symbolic API to be used by people who know
>> > 0.0 about the FEM.
>
> The only connection to FEM is that we need to define the weak
> formulation (NOT a plain PDE). This is because several different
> weak formulations can usually be obtained from a PDE,
> and second, that the weak formulation has everything incorporated
> in it (boundary conditions etc.), while they have to be presented
> separately with a standard PDE. Third, the weak formulation
> is the basis for FEM...
The weak form works fine for my applications as well. But the great
thing about using sympy is that in principle we could have sympy
generate the weak form from the original PDE + boundary conditions.
In some contexts (like working with undergraduate students in classes)
that would be nice. But that is a detail - for now, using sympy to
specific the weak form itself is a great direction.
Cheers,
Brian
Brian Granger
you include the boundary condition information, they are equivalent.
Cheers,
Brian
Andy R. Terrel
Sorry to disturb your nice series of thoughts here, but I was curious
how Ondrej planned on using sympy to represent forms. Since I code
with FEniCS quite a bit I'll shoot some response to your questions.
On May 26, 4:59 am, Robert Cimrman <cimrm...@ntc.zcu.cz> wrote:
> Pavel Solin wrote:
> > Has Fenics nonsimplicial elements (quads, bricks, pyramids)?
The libraries are general enough that if someone wanted to add these
other elements it is possible, but there has been no push by the user/
developer community.
> > I just wonder since for these elements it is not enough to
> > precompute a single local stiffness matrix on the reference
> > element and use it for all elements in the mesh. This was the
> > original concept of Fenics - they spent lots of time optimizing
> > some higher level stuff that was built on this naive foundation.
> > I am very curious whether they still stick to this original concept.
(FFC). It has recently been revamped quite a bit to also generate
optimize quadrature code for the forms because there were far too many
forms that FFC tensor contractions were not optimal on. As an aside
FFC has been marked as a backend to another more general form compiler
that has just been release Unified Form Compiler (UFC) that does use a
richer symbolic engine for things like differentiation and what not.
>
> I cannot tell you such details as I do not know them :-) When I last
> looked, the foundations was (IMHO) just as you say. I just wanted you to
> look at Fenics, as it seems similar in goals, so there _might_ be ideas
> useful for us - it is the only code that employs symbolic computations
> for FE.
html/) which also support automatic differentiation and linearization
of forms, but then again only simplices. Or if you are talking about
*only* symbolic, that library would be very slow indeed.
>
> As for the matrices, I do not think they have to be precomputed, Fenics
> can (IMHO) translate that to an optimized C(++?) code. The symbolic
> stuff serves as a foactory for making dynamic C libraries of optimized
> element matrices on the fly. But I can be wrong...
generate C++ code and then use this code in your library (usually
DOLFIN but it conforms to a api called UFL (unified form language) )
or use pyDOLFIN and have the code generated and compiled on the fly
(in C++ called from python). UFC does have the ability to take a non-
linear form and produce the linearized left and right hand sides, FFC
does not since it turns to numeric methods very quickly in the
process.
-- Andy
Ondrej Certik
thanks for joining the discussion.
On Sat, May 30, 2009 at 4:53 PM, Andy R. Terrel <andy....@gmail.com> wrote:
>
>
> Sorry to disturb your nice series of thoughts here, but I was curious
> how Ondrej planned on using sympy to represent forms. Since I code
> with FEniCS quite a bit I'll shoot some response to your questions.
>
>
> On May 26, 4:59 am, Robert Cimrman <cimrm...@ntc.zcu.cz> wrote:
>> Pavel Solin wrote:
>> > Has Fenics nonsimplicial elements (quads, bricks, pyramids)?
>
> No FEniCS has remained simplicial, probably its biggest fault IMHO.
> The libraries are general enough that if someone wanted to add these
> other elements it is possible, but there has been no push by the user/
> developer community.
>
>> > I just wonder since for these elements it is not enough to
>> > precompute a single local stiffness matrix on the reference
>> > element and use it for all elements in the mesh. This was the
>> > original concept of Fenics - they spent lots of time optimizing
>> > some higher level stuff that was built on this naive foundation.
>> > I am very curious whether they still stick to this original concept.
>
> The precomputing was the grand challenge for the FEniCS Form Compiler
> (FFC). It has recently been revamped quite a bit to also generate
> optimize quadrature code for the forms because there were far too many
> forms that FFC tensor contractions were not optimal on. As an aside
> FFC has been marked as a backend to another more general form compiler
> that has just been release Unified Form Compiler (UFC) that does use a
> richer symbolic engine for things like differentiation and what not.
So do you think it's possible to implement quads and bricks into
fenics? The whole assembly would have to be modified, if I understand
correctly.
>
>>
>> I cannot tell you such details as I do not know them :-) When I last
>> looked, the foundations was (IMHO) just as you say. I just wanted you to
>> look at Fenics, as it seems similar in goals, so there _might_ be ideas
>> useful for us - it is the only code that employs symbolic computations
>> for FE.
>
> Not true. Look at Sundance (http://www.math.ttu.edu/~klong/Sundance/
> html/) which also support automatic differentiation and linearization
> of forms, but then again only simplices. Or if you are talking about
> *only* symbolic, that library would be very slow indeed.
>
>>
>> As for the matrices, I do not think they have to be precomputed, Fenics
>> can (IMHO) translate that to an optimized C(++?) code. The symbolic
>> stuff serves as a foactory for making dynamic C libraries of optimized
>> element matrices on the fly. But I can be wrong...
>
> There are two ways to use the system. Use UFC (or FFC or SyFi) to
> generate C++ code and then use this code in your library (usually
> DOLFIN but it conforms to a api called UFL (unified form language) )
> or use pyDOLFIN and have the code generated and compiled on the fly
> (in C++ called from python). UFC does have the ability to take a non-
> linear form and produce the linearized left and right hand sides, FFC
> does not since it turns to numeric methods very quickly in the
> process.
I really don't like that you have to compile it each time --- I
believe it must be possible to do design it in a way, so that you
don't have to compile anything, but still be able to change forms at
runtime.
Ondrej
Andy Ray Terrel
it is possible and yes the current software would require some better
abstractions, but the essential algorithms have been discovered. A
simple form compiler + element generator + topology library + linear
algebra/polynomial solver works really well. A colleague and I were
working with a system that basically tapped into some of FEniCS
middleware where you describe the topology and we provide the
numerical management system, but we only used Lagrange elements and
didn't polish much before we had to work on some other projects.
change the form, otherwise you have to make an interpreter and lose
out on the last 30 years of compiler theory*. pyDolfin gives the best
of both worlds, using just in time compiling and caching old forms.
* Although you could argue, as Marcus Püschel does, that compilers
don't help automatically generated code too much but the interpreter
has to be fast.
>
> Ondrej
>
> >
>
Pavel Solin
thanks a lot for contributing to the discussion.
I think that we met briefly in Bloomington.
See some comments/questions below.
When I saw FEniCS last time, it seems to me that it was built
really tightly on the Lagrange concept. Is this still the case, or
can you define shape functions in a more general way?
This reminds me of the BMW effort to re-use their original
regular gas engines for liquid hydrogen :) But honestly,
we hope to be able to figure this out, for us it is quite
important.
Best,
Pavel
pyDolfin gives the best
of both worlds, using just in time compiling
and caching old forms.
* Although you could argue, as Marcus Püschel does, that compilers
don't help automatically generated code too much but the interpreter
has to be fast.
>
> Ondrej
>
> >
>
--
Pavel Solin
University of Nevada, Reno
Hermes project: http://hpfem.org/
FemHub project: http://femhub.org/
Ondrej Certik
>> > believe it must be possible to do design it in a way, so that you
>> > don't have to compile anything, but still be able to change forms at
>> > runtime.
>>
>> So the compiled library will run faster. You only compile when you
>> change the form, otherwise you have to make an interpreter and lose
>> out on the last 30 years of compiler theory*.
>
>
> This reminds me of the BMW effort to re-use their original
> regular gas engines for liquid hydrogen :) But honestly,
> we hope to be able to figure this out, for us it is quite
> important.
I at least don't want to compile things for most of the frequent
equations. I have nothing against automatically generating some C++
stuff, like the shapefunctions in hermes, and in principle the forms
can be generated too.
My idea is just that this autogeneration should be done only at
compile time (e.g. when I build the FEM package on my computer and
upload a binary), not at runtime (when the user downloads the binary
and use it to solve his PDEs).
Ondrej
Ondrej Certik
I guess a solution is to just ship precompiled modules for the most
common equations, which are ready to be used in python, and provide
the user a way to compile his own advanced stuff, if needed.
Ondrej
Andy Ray Terrel
> Hi Andy,
> thanks a lot for contributing to the discussion.
> I think that we met briefly in Bloomington.
> See some comments/questions below.
having trouble with Taylor Hood hp methods. Did you ever get this
fixed? Another possibility that I would have thought unfruitful at
the time would be to use SUPG stabilization with equal order pressure
and velocity. Anywho that's a different conversation.
quite a few more elements such as Raviart Thomas, Nedelec, BDFM, BDM,
it is pretty easy to define anything based on some orthogonal
polynomial generators. Syfi supports a lot more but it does
everything symbolically so its quite a bit easier to define the
elements. Both generate code that can be used by the simulation
engine of your choice and then the correct mappings are required, but
DOLFIN supports most of these that can be supported in a reference
element way.
the first to jump on a system that could give me great interpreted
numbers. It would really help out the automatic differentiation
people as well since writing ad code in an interpreted environment is
much simplier. But I would be willing to bet you could stand the
penalty of compiling the form each time with some tricks to get the
compile time to some minimal time, such as using precompiled headers.
What's your use case for not recompiling the form?
-- Andy
Pavel Solin
Yes that was me. At the meeting you mentioned that your group was
On Sat, May 30, 2009 at 7:18 PM, Pavel Solin <so...@unr.edu> wrote:
> Hi Andy,
> thanks a lot for contributing to the discussion.
> I think that we met briefly in Bloomington.
> See some comments/questions below.
having trouble with Taylor Hood hp methods. Did you ever get this
fixed?
Probably, the code was flawed, it works pretty well now (I think).
But for some reason, hp-adaptivity is not so dramatically better
than usual adaptivity for time-dependent N-S equations. We may
need to invest more time into this.
Yes, Lagrange elements, but also other nodal elements (defined
via linear forms in the Ciarlet sense). Nedelec elements are also
nodal elements for me. I am more interested in elements where
you can define your shape functions directly, without going
through the linear forms. From what you say i understand that
you have this option.
We want Hermes to be usable Matlab-style 100% percent,
so no compilation at all if possible. Perhaps we will end up
with a combination, and allow the user to compile the
forms for best performance, after, say, testing them in
some slower format which does not require compilation.
Hey I would be
the first to jump on a system that could give me great interpreted
numbers. It would really help out the automatic differentiation
people as well since writing ad code in an interpreted environment is
much simplier. But I would be willing to bet you could stand the
penalty of compiling the form each time with some tricks to get the
compile time to some minimal time, such as using precompiled headers.
What's your use case for not recompiling the form?
Best,
Pavel
-- Andy
>
> Best,
>
> Pavel
>
>
>
>>
>> pyDolfin gives the best
>> of both worlds, using just in time compiling
>
>
>>
>> and caching old forms.
>>
>> * Although you could argue, as Marcus Püschel does, that compilers
>> don't help automatically generated code too much but the interpreter
>> has to be fast.
>>
>> >
>> > Ondrej
>> >
>> > >
>> >
>>
>>
>
>
>
> --
> Pavel Solin
> University of Nevada, Reno
> http://hpfem.math.unr.edu/people/pavel/
> Hermes project: http://hpfem.org/
> FemHub project: http://femhub.org/
>
> >
>
Ondrej Certik
[...]
>> I'm sorry which is important, not compiling the form?
>
> We want Hermes to be usable Matlab-style 100% percent,
> so no compilation at all if possible. Perhaps we will end up
> with a combination, and allow the user to compile the
> forms for best performance, after, say, testing them in
> some slower format which does not require compilation.
Exactly.
Ondrej
Brian Granger
>> so no compilation at all if possible. Perhaps we will end up
>> with a combination, and allow the user to compile the
>> forms for best performance, after, say, testing them in
>> some slower format which does not require compilation.
>
> Exactly.
How will you avoid compiling without things becoming horribly slow?
For the Schroedinger equation, there are basically infinitely many
potentials (this is no exageration!) with very complex algebraic
forms. These are the types of problems that I am particularly
interested in.
I don't see how you could cover all these with a pre-compiled
approach. The way I imagined using sympy to specify weak forms was
that you would take the sympy tree and then have it generate fast c++
code for everything that is compiled on the fly. The user wouldn't
know about the compilation to the experience is still 100% Matlab-ish.
Thoughts?
Brian
Ondrej Certik
I think we should just start coding it and we'll see. If it's
possible, let's not compile, if it's not possible, we'll have to
compile on the fly.
Ondrej
Robert Cimrman
The point is, apart having a library of precompiled modules (which never
encompasses everything a (power-)user might want) to allow also a
definition/compilation on the fly. This is a good thing IMHO. Of course,
the standard library should be rich enough for most regular users, so
that they do not need to bother.
r.
Robert Cimrman
Andy R. Terrel wrote:
>
> Sorry to disturb your nice series of thoughts here, but I was curious
> how Ondrej planned on using sympy to represent forms. Since I code
> with FEniCS quite a bit I'll shoot some response to your questions.
Great, more heads - more thoughts :) It good that someone actually
knowing Fenics chimed in.
> On May 26, 4:59 am, Robert Cimrman <cimrm...@ntc.zcu.cz> wrote:
>> Pavel Solin wrote:
>>> Has Fenics nonsimplicial elements (quads, bricks, pyramids)?
>
> No FEniCS has remained simplicial, probably its biggest fault IMHO.
> The libraries are general enough that if someone wanted to add these
> other elements it is possible, but there has been no push by the user/
> developer community.
>
>>> I just wonder since for these elements it is not enough to
>>> precompute a single local stiffness matrix on the reference
>>> element and use it for all elements in the mesh. This was the
>>> original concept of Fenics - they spent lots of time optimizing
>>> some higher level stuff that was built on this naive foundation.
>>> I am very curious whether they still stick to this original concept.
>
> The precomputing was the grand challenge for the FEniCS Form Compiler
> (FFC). It has recently been revamped quite a bit to also generate
> optimize quadrature code for the forms because there were far too many
> forms that FFC tensor contractions were not optimal on. As an aside
> FFC has been marked as a backend to another more general form compiler
> that has just been release Unified Form Compiler (UFC) that does use a
> richer symbolic engine for things like differentiation and what not.
So the UFC is a natural candidate to look at when thinking about
implementing similar technizues into SymPy, yes? I guess it is not tight
to the particular shape functions (simplices) on the abstract level of
weak forms, so another bases can be plugged in?
>> I cannot tell you such details as I do not know them :-) When I last
>> looked, the foundations was (IMHO) just as you say. I just wanted you to
>> look at Fenics, as it seems similar in goals, so there _might_ be ideas
>> useful for us - it is the only code that employs symbolic computations
>> for FE.
>
> Not true. Look at Sundance (http://www.math.ttu.edu/~klong/Sundance/
> html/) which also support automatic differentiation and linearization
> of forms, but then again only simplices. Or if you are talking about
> *only* symbolic, that library would be very slow indeed.
I am glad to be led out of mistake, Sundance seems very interesting,
too. Thank you for the link.
>> As for the matrices, I do not think they have to be precomputed, Fenics
>> can (IMHO) translate that to an optimized C(++?) code. The symbolic
>> stuff serves as a foactory for making dynamic C libraries of optimized
>> element matrices on the fly. But I can be wrong...
>
> There are two ways to use the system. Use UFC (or FFC or SyFi) to
> generate C++ code and then use this code in your library (usually
> DOLFIN but it conforms to a api called UFL (unified form language) )
> or use pyDOLFIN and have the code generated and compiled on the fly
> (in C++ called from python). UFC does have the ability to take a non-
> linear form and produce the linearized left and right hand sides, FFC
> does not since it turns to numeric methods very quickly in the
> process.
Nice food for thought, thanks!
r.
Robert Cimrman
Yeah, talk is cheap (especially mine...). Just once we have a C++ code
generator, there is _no_ theoretical/practical obstacle why it should
not be enabled (optionally) also for the on-the-fly compiling. Everybody
would be free to use it or not.
While I agree with the Matlab (or better SciPy ;)) analogy is good, it
seems to me that not using the (big IMHO) easy possibility of the
on-the-fly mode once the code generation works equals to "crippling" our
possibilities deliberately.
r.
Andy Ray Terrel
is the form language. But yes it is good candidate to look at,
http://fenics.org/wiki/UFL . My understanding is that it needs to
call a form compiler underneath it so if the form compiler doesn't
support the shape functions you like then you have to implement them,
but it handles all the form manipulations. One could conceivable do
all the backend support in SymPy even and then you would be able to
take any UFL form for the simulation.
The FEniCS core group has announced that they are trying to get
everything hardened up for a 1.0 release in October.
-- Andy
Vinzent Steinberg
On 31 Mai, 11:29, Pavel Solin <so...@unr.edu> wrote:
> We want Hermes to be usable Matlab-style 100% percent,
> so no compilation at all if possible. Perhaps we will end up
> with a combination, and allow the user to compile the
> forms for best performance, after, say, testing them in
> some slower format which does not require compilation.
usage. You could use for instance libtcc to compile C completely in
memory. I tried this (see sympy/utilities/compilef.py) and in some
cases it's faster for function evaluation than numpy, including the
overhead of compilation.
Vinzent
Pavel Solin
thank you for the valuable information. Hermes has
a powerful adaptive FEM engine that can be used for
arbitrary PDE problems, but many people warned us to
stay away from a very complex multi-functional GUI. So,
we would like to develop a tool ("super-GUI" ?) which
allows us (or the user eventually) to create low-level
problem-specific GUIs easily. The weak forms and/or other
necessary things would be compiled one time only, when
a low-level GUI is created/changed, not when it is used
for computing. Pavel Karban just started to develop
a GUI in Qt, but right now it is closer to a low-level GUI
than to the "super-GUI" described above. I am not sure
whether the "super GUI" can be written in Qt or if some
other library/language would be better for that. Since
this is among the highest priorities for us right now so
any suggestions or even help are extremely welcome!
With best regards,
Pavel
Pavel Solin
University of Nevada, Reno
http://hpfem.math.unr.edu/people/pavel/
Hermes project: http://hpfem.org/
Ondrej Certik
> Hello Vinzent,
> thank you for the valuable information. Hermes has
> a powerful adaptive FEM engine that can be used for
> arbitrary PDE problems, but many people warned us to
> stay away from a very complex multi-functional GUI. So,
> we would like to develop a tool ("super-GUI" ?) which
> allows us (or the user eventually) to create low-level
> problem-specific GUIs easily. The weak forms and/or other
> necessary things would be compiled one time only, when
> a low-level GUI is created/changed, not when it is used
> for computing. Pavel Karban just started to develop
> a GUI in Qt, but right now it is closer to a low-level GUI
> than to the "super-GUI" described above. I am not sure
> whether the "super GUI" can be written in Qt or if some
> other library/language would be better for that. Since
> this is among the highest priorities for us right now so
> any suggestions or even help are extremely welcome!
From what you say it seems to me that if there is some automatic
generator of the form, that has to be recompiled, it won't be a
problem, because this would be done once (per problem) anyway, then
one would use the GUI.
Ondrej
Pavel Solin
Pavel