Re: [firedrake] Tetrahedral mesh
Dear David, I think you have confused what Lawrence was asking. Since you are working with Sebastian Reich, we assume you are solving atmospheric problems. In this context, the definite Helmholtz equation usually occurs: laplacian u - u = f while you appear to be solving the indefinite Helmholtz equation: laplacian u + u = f We were just attempting to ascertain whether you were deliberately solving the indefinite version of the problem. Is this the case? David On Thu, 10 Dec 2015 at 14:47 Angwenyi David <kipkoej@gmail.com> wrote:
Integration by parts is equivalent to using Green’s function; and so it turns out as follows:
laplacian u = - u + f Greens function: integ((laplacian u)v) + integ(dot(grad u, grad v)) = surfinteg((partial_n u)v)
So that the weak form is obtained by substituting laplacian u = -u + f from the Helmholtz equation in the Green’s function. The same is arrived at by integration by parts.
Sincerely,
David
On 10 Dec 2015, at 15:23, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 10/12/15 14:12, Angwenyi David wrote:
nx = 10 ny = 10 nz = 20 mesh = BoxMesh(nx, ny, nz, 1.0, 1.0, 2.0) V = FunctionSpace(mesh, "CG", 1) u = TrialFunction(V) v = TestFunction(V) f = Function(V) f.interpolate(Expression("(2*x[0])+(3*x[1])+(0.5*x[2])")) Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', tetrahedron, 1, dim=3), 3), FiniteElement('Lagrange', tetrahedron, 1)), 4) a = (dot(grad(v), grad(u)) - v * u) * dx
This is the indefinite form of the Helmholtz equation. Did you actually want the sign-positive Helmholtz equation?
The later has strong form:
-laplacian u + u = f
When integrating by parts the weak form becomes:
<grad u, grad v> + <u, v> = <f, v>
As to incorporation of boundary conditions, see
http://firedrakeproject.org/variational-problems.html#incorporating-boundary...
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Dear David, I was trying to solve the indefinite one. This is just an attempt to understand the working of Firedrake. We hope to solve PDEs on a sphere, for example the Shallow Water Equation on Icosahedral grids with time; and I am glad that such a geometry can be worked on using Firedrake. David
On 10 Dec 2015, at 15:53, David Ham <David.Ham@imperial.ac.uk> wrote:
Dear David,
I think you have confused what Lawrence was asking. Since you are working with Sebastian Reich, we assume you are solving atmospheric problems. In this context, the definite Helmholtz equation usually occurs:
laplacian u - u = f
while you appear to be solving the indefinite Helmholtz equation:
laplacian u + u = f
We were just attempting to ascertain whether you were deliberately solving the indefinite version of the problem. Is this the case?
David
On Thu, 10 Dec 2015 at 14:47 Angwenyi David <kipkoej@gmail.com <mailto:kipkoej@gmail.com>> wrote: Integration by parts is equivalent to using Green’s function; and so it turns out as follows:
laplacian u = - u + f Greens function: integ((laplacian u)v) + integ(dot(grad u, grad v)) = surfinteg((partial_n u)v)
So that the weak form is obtained by substituting laplacian u = -u + f from the Helmholtz equation in the Green’s function. The same is arrived at by integration by parts.
Sincerely,
David
On 10 Dec 2015, at 15:23, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk <mailto:lawrence.mitchell@imperial.ac.uk>> wrote:
On 10/12/15 14:12, Angwenyi David wrote:
nx = 10 ny = 10 nz = 20 mesh = BoxMesh(nx, ny, nz, 1.0, 1.0, 2.0) V = FunctionSpace(mesh, "CG", 1) u = TrialFunction(V) v = TestFunction(V) f = Function(V) f.interpolate(Expression("(2*x[0])+(3*x[1])+(0.5*x[2])")) Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', tetrahedron, 1, dim=3), 3), FiniteElement('Lagrange', tetrahedron, 1)), 4) a = (dot(grad(v), grad(u)) - v * u) * dx
This is the indefinite form of the Helmholtz equation. Did you actually want the sign-positive Helmholtz equation?
The later has strong form:
-laplacian u + u = f
When integrating by parts the weak form becomes:
<grad u, grad v> + <u, v> = <f, v>
As to incorporation of boundary conditions, see
http://firedrakeproject.org/variational-problems.html#incorporating-boundary... <http://firedrakeproject.org/variational-problems.html#incorporating-boundary-conditions>
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake> _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 10/12/15 15:06, Angwenyi David wrote:
Dear David,
I was trying to solve the indefinite one.
Ah, OK. A few things: 1. Your operator is indefinite, so a ksp_type of 'cg' won't work (conjugate gradient only works for SPD operators), try using 'minres' instead (which works for symmetric indefinite systems). 2. The indefinite Helmholtz problem is notoriously ill-conditioned, and moreover, very difficult to precondition well. If you're not trying to solve problems that are very large, try using a direct solver: {'pc_type': 'lu', 'ksp_type': 'preonly'} Cheers, Lawrence
Thank you!
On 10 Dec 2015, at 16:11, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 10/12/15 15:06, Angwenyi David wrote:
Dear David,
I was trying to solve the indefinite one.
Ah, OK. A few things:
1. Your operator is indefinite, so a ksp_type of 'cg' won't work (conjugate gradient only works for SPD operators), try using 'minres' instead (which works for symmetric indefinite systems).
2. The indefinite Helmholtz problem is notoriously ill-conditioned, and moreover, very difficult to precondition well. If you're not trying to solve problems that are very large, try using a direct solver:
{'pc_type': 'lu', 'ksp_type': 'preonly'}
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Lawrence, The indefinite one works very well with ‘minres’. Thanks again. David
On 10 Dec 2015, at 16:19, Angwenyi David <kipkoej@gmail.com> wrote:
Thank you!
On 10 Dec 2015, at 16:11, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 10/12/15 15:06, Angwenyi David wrote:
Dear David,
I was trying to solve the indefinite one.
Ah, OK. A few things:
1. Your operator is indefinite, so a ksp_type of 'cg' won't work (conjugate gradient only works for SPD operators), try using 'minres' instead (which works for symmetric indefinite systems).
2. The indefinite Helmholtz problem is notoriously ill-conditioned, and moreover, very difficult to precondition well. If you're not trying to solve problems that are very large, try using a direct solver:
{'pc_type': 'lu', 'ksp_type': 'preonly'}
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Angwenyi David
-
David Ham
-
Lawrence Mitchell