On 2013-10-31 11:25, Mikael Mortensen wrote:
Den Oct 31, 2013 kl. 11:20 AM skrev Garth N. Wells:
On 2013-10-31 09:42, Mikael Mortensen wrote:
Dear all,
I have been inspired by the new LocalSolver class of Garth and I'm
currently trying out ideas to compute "projections" locally without
going through a global coefficient matrix and a global linear algebra
solve. I just wanted to ask if anyone has any experience with this?
If
so, I would be very interested to hear about it. I'm particularly
interested in speeding up the following:
V = FunctionSpace(mesh, "CG", 1)
p = Function(V)
dpdx = project(p.dx(0), V)
My first efforts on using a local averaging operator has been put in
branch
https://bitbucket.org/fenics-project/dolfin/branch/mikael/local-average-operator
[1]. For small problems I have found it to be 4-5 times faster than
the regular project shown above, and the accuracy is vey similar.
This can already be done with LocalSolver. For projection onto a
continuous basis, it is mathematically unsound because it depends on
the order in which cells are visited, but it is useful for
post-processing very large problems.
I have started this much because the LocalSolver is unsound (and
inaccurate) for a continuous basis.
I want a sound local solver that does not depend on the order of how
cells are visited. For a continuous basis, LocalSolver could possibly
be extended for dofs on the interface between cells. That would
require some local averaging, much like in the branch I've started.
Currently, LocalSolver is not accurate enough to be used for anything
but visualization of a continuous basis. If I want to project onto a
discontinuous basis, then sure it's great, but I want projection onto
P1!
I don't see how this is possible. An L2 projection onto a continuous
basis is fundamentally a global operation. Can you write down
mathematically what you want to compute?
For projection onto a discontinuous basis, LocalSolver gives the same
result as the global projection algorithm, but uses virtually no
memory and I would expect is more than 4-5 times faster.
It is about twice as fast in my tests for rather small problems.
For cell-wise averages, just use LocalSolver and a P0 basis.
Sure.
I don't think we should make simple wrappers for unsound operations,
even when the operation can be useful. A user should 'know' when
they're doing something that's fishy. We could consider
V = FunctionSpace(mesh, "CG", 1)
p = Function(V)
dpdx = local_project(p.dx(0)) # project onto P0 basis
But I want projection onto P1:-)
I need p.dx(0) onto P1 because in Navier-Stokes the velocity is updated
through
u_updated = u_old + dt*grad(p)
Regular solve is then something like
Solve(inner(u, v)*dx == inner(u_old, v)*dx + dt*inner(grad(p), v)*dx)
If I could project grad(p) onto the continuous velocity space I could
instead do:
u_update.vector()[:] = u_old.vector()[:]
u_updated.vector().axpy(dt, grad_p.vector())
which would incur a considerable speed-up.
These are errornorms of projecting p.dx(0) onto P1 on the UnitSquare
for regular project, LocalAverageOperator (my branch) and LocalSolver.
p = sin(x)*cos(y) and results are compared with exact solution.
Errornorm of projecting p.dx(0) onto P1
Dofs Regular project LAO LocalSolver
36 1.160083e-02 1.125136e-02 1.342422e-01
256 1.517463e-03 2.001789e-03 5.091767e-02
676 5.743440e-04 8.044087e-04 3.077231e-02
1296 3.013552e-04 4.353792e-04 2.200806e-02
2116 1.858203e-04 2.739904e-04 1.712131e-02
3136 1.261956e-04 1.888708e-04 1.400793e-02
4356 9.139609e-05 1.383886e-04 1.185168e-02
5776 6.930475e-05 1.059355e-04 1.027027e-02
7396 5.439590e-05 8.380729e-05 9.060989e-03
9216 4.385459e-05 6.802556e-05 8.106367e-03
So there is something to it, right?
Averaging can be ok, but I'd like to see a precise description of the
operation you want to perform.
Garth
Mikael
Garth
Best regards
Mikael
------------
Dr. Mikael Mortensen
Associate Professor
Department of Mathematics
Mechanics Division
University of Oslo
Phone: +47-22855866
Mob: +47-41407201
Email: mik...@math.uio.no
Links:
------
[1]
https://bitbucket.org/fenics-project/dolfin/branch/mikael/local-average-operator
_______________________________________________
fenics mailing list
fenics@fenicsproject.org
http://fenicsproject.org/mailman/listinfo/fenics
_______________________________________________
fenics mailing list
fenics@fenicsproject.org
http://fenicsproject.org/mailman/listinfo/fenics
