Re: [firedrake] Vertical normal in operator assembly
Hi Eike, A shortcut is to use ast.FlatBlock("put your C code here") and then pass it to op2.Kernel. Let me know if this helps. -- Fabio 2014-12-05 12:03 GMT+00:00 Eike Mueller <e.mueller@bath.ac.uk>:
...looking through Lawrence’s emails on the 3d solver, I found the code below, which seems to build the zhat field. It looks pretty impenetrable to me, is there any easier way of doing this or do I only need parts of it?
Thanks,
Eike
### BEGIN NORMALS BLOCK ###
# Build new extruded coordinate function space normal_fs = VectorFunctionSpace(mesh, "DG", 0) zhat = Function(normal_fs)
base_coords = mesh._old_mesh.coordinates base_coords_dim = base_coords.function_space().cdim base_coords_arity = base_coords.cell_node_map().arity expr = Expression(('x[0]', 'x[1]', 'x[2]')) v0 = lambda x: ast.Symbol("v0", (x,)) v1 = lambda x: ast.Symbol("v1", (x,)) n = lambda x: ast.Symbol("n", (x,)) x = lambda x: ast.Symbol("x", (x,)) coords = lambda x, y: ast.Symbol("base_coords", (x, y)) normals = lambda x, y: ast.Symbol("normals", (x, y)) body = [] body += [ast.Decl("double", v(3)) for v in [v0, v1, n, x]] body.append(ast.Decl("double", "dot")) body.append(ast.Decl("double", "norm")) body.append(ast.Assign("norm", 0.0)) body.append(ast.Assign("dot", 0.0)) body.append(ast.Decl("int", "i")) body.append(ast.Decl("int", "c")) body.append(ast.Decl("int", "d")) body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Assign(v0("i"), ast.Sub(coords(1, "i"), coords(0, "i"))), ast.Assign(v1("i"), ast.Sub(coords(2, "i"), coords(0, "i"))), ast.Assign(x("i"), 0.0)])) # n = v0 x v1 body.append(ast.Assign(n(0), ast.Sub(ast.Prod(v0(1), v1(2)), ast.Prod(v0(2), v1(1))))) body.append(ast.Assign(n(1), ast.Sub(ast.Prod(v0(2), v1(0)), ast.Prod(v0(0), v1(2))))) body.append(ast.Assign(n(2), ast.Sub(ast.Prod(v0(0), v1(1)), ast.Prod(v0(1), v1(0)))))
body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Incr(x(j), coords("i", j)) for j in range(3)]))
body.extend([ast.FlatBlock("dot += (%(x)s) * n[%(i)d];\n" % {"x": x, "i": i}) for i, x in enumerate(expr.code)])
body.extend([ast.FlatBlock("norm += n[%(i)d] * n[%(i)d];\n" % {"i": i}) for i in range(3)])
body.append(ast.FlatBlock("norm = sqrt(norm);\n"))
# normals[d][c] = n[c] / norm; assign_block = [ast.Assign(normals("d", "c"), ast.Div(n("c"), ast.Symbol("norm")))] body.append(ast.IMul("norm", ast.Ternary(ast.Less("dot", 0), "-1", "1"))) body.append(ast.For(ast.Assign("d", 0), ast.Less("d", 1), ast.Incr("d", 1), [ast.For(ast.Assign("c", 0), ast.Less("c", base_coords_dim), ast.Incr("c", 1), assign_block)]))
# Looks complicated, actually just produces the code below # void normal_calc (double** base_coords , double** normals ) # { # double v0[3] ; # double v1[3] ; # double n[3] ; # double x[3] ; # double dot ; # double norm ; # norm = 0.0; # dot = 0.0; # int i ; # int c ; # int d ; # for (i = 0; i < 3; ++i) # { # v0[i] = base_coords[1][i] - base_coords[0][i]; # v1[i] = base_coords[2][i] - base_coords[0][i]; # x[i] = 0.0; # } # n[0] = v0[1] * v1[2] - v0[2] * v1[1]; # n[1] = v0[2] * v1[0] - v0[0] * v1[2]; # n[2] = v0[0] * v1[1] - v0[1] * v1[0]; # for (i = 0; i < 3; ++i) # { # x[0] += base_coords[i][0]; # x[1] += base_coords[i][1]; # x[2] += base_coords[i][2]; # } # dot += (x[0]) * n[0]; # dot += (x[1]) * n[1]; # dot += (x[2]) * n[2]; # norm += n[0] * n[0]; # norm += n[1] * n[1]; # norm += n[2] * n[2]; # norm = sqrt(norm); # norm *= (dot < 0 ? -1 : 1); # for (d = 0; d < 1; ++d) # { # for (c = 0; c < 3; ++c) # { # normals[d][c] = n[c] / norm; # } # } # }
kernel = op2.Kernel(ast.FunDecl("void", "normal_calc", [ast.Decl("double**", "base_coords"), ast.Decl("double**", "normals")], ast.Block(body)), "normal_calc")
base_coords = mesh._old_mesh.coordinates op2.par_loop(kernel, zhat.cell_set, base_coords.dat(op2.READ, base_coords.cell_node_map()), zhat.dat(op2.WRITE, zhat.cell_node_map()))
### END NORMALS BLOCK ###
On 5 Dec 2014, at 10:01, Eike Mueller <E.Mueller@bath.ac.uk> wrote:
Dear firedrakers,
I have the following function spaces:
W_2 = HDiv(U_2 x V_0) + HDiv(U_1 x V_1) [velocity] W_2^b = U_2 x V_0 [buoyancy]
where U_0 -> U_1 -> U_2 and V_0 -> V_1 are horizontal- and vertical de Rham complexes.
I need to be able to apply the operator Q: W_2^b -> W_2 defined by
Q_{ij} = <\vec{w}_i,\hat{z} \gamma_j>
where \vec{w}_i and \gamma_j are basis functions and \hat{z} is the unit normal in the vertical direction. To apply Q to a buoyancy field b, I want to write something like
u = assemble(dot(TestFunction(W_2),zhat*b)*dx)
but how do I get the zhat in there? Is that even possible?
Thanks a lot,
Eike
--
Dr Eike Hermann Mueller Research Associate (PostDoc)
Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom
+44 1225 38 5803 e.mueller@bath.ac.uk http://people.bath.ac.uk/em459/
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Fabio, thanks, I guess I could also just use the code that Lawrence commented out below… I figured out what the C code does, and that’s what I’m doing now. In a radially extruded mesh, are the normals in one column always the same? I think this is what the code assumes, since it only uses the coordinates in the base cell. Cheers, Eike -- Dr Eike Hermann Mueller Research Associate (PostDoc) Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5803 e.mueller@bath.ac.uk http://people.bath.ac.uk/em459/
On 5 Dec 2014, at 12:55, Fabio Luporini <f.luporini12@imperial.ac.uk> wrote:
Hi Eike,
A shortcut is to use ast.FlatBlock("put your C code here") and then pass it to op2.Kernel. Let me know if this helps.
-- Fabio
2014-12-05 12:03 GMT+00:00 Eike Mueller <e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk>>: ...looking through Lawrence’s emails on the 3d solver, I found the code below, which seems to build the zhat field. It looks pretty impenetrable to me, is there any easier way of doing this or do I only need parts of it?
Thanks,
Eike
### BEGIN NORMALS BLOCK ###
# Build new extruded coordinate function space normal_fs = VectorFunctionSpace(mesh, "DG", 0) zhat = Function(normal_fs)
base_coords = mesh._old_mesh.coordinates base_coords_dim = base_coords.function_space().cdim base_coords_arity = base_coords.cell_node_map().arity expr = Expression(('x[0]', 'x[1]', 'x[2]')) v0 = lambda x: ast.Symbol("v0", (x,)) v1 = lambda x: ast.Symbol("v1", (x,)) n = lambda x: ast.Symbol("n", (x,)) x = lambda x: ast.Symbol("x", (x,)) coords = lambda x, y: ast.Symbol("base_coords", (x, y)) normals = lambda x, y: ast.Symbol("normals", (x, y)) body = [] body += [ast.Decl("double", v(3)) for v in [v0, v1, n, x]] body.append(ast.Decl("double", "dot")) body.append(ast.Decl("double", "norm")) body.append(ast.Assign("norm", 0.0)) body.append(ast.Assign("dot", 0.0)) body.append(ast.Decl("int", "i")) body.append(ast.Decl("int", "c")) body.append(ast.Decl("int", "d")) body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Assign(v0("i"), ast.Sub(coords(1, "i"), coords(0, "i"))), ast.Assign(v1("i"), ast.Sub(coords(2, "i"), coords(0, "i"))), ast.Assign(x("i"), 0.0)])) # n = v0 x v1 body.append(ast.Assign(n(0), ast.Sub(ast.Prod(v0(1), v1(2)), ast.Prod(v0(2), v1(1))))) body.append(ast.Assign(n(1), ast.Sub(ast.Prod(v0(2), v1(0)), ast.Prod(v0(0), v1(2))))) body.append(ast.Assign(n(2), ast.Sub(ast.Prod(v0(0), v1(1)), ast.Prod(v0(1), v1(0)))))
body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Incr(x(j), coords("i", j)) for j in range(3)]))
body.extend([ast.FlatBlock("dot += (%(x)s) * n[%(i)d];\n" % {"x": x, "i": i}) for i, x in enumerate(expr.code)])
body.extend([ast.FlatBlock("norm += n[%(i)d] * n[%(i)d];\n" % {"i": i}) for i in range(3)])
body.append(ast.FlatBlock("norm = sqrt(norm);\n"))
# normals[d][c] = n[c] / norm; assign_block = [ast.Assign(normals("d", "c"), ast.Div(n("c"), ast.Symbol("norm")))] body.append(ast.IMul("norm", ast.Ternary(ast.Less("dot", 0), "-1", "1"))) body.append(ast.For(ast.Assign("d", 0), ast.Less("d", 1), ast.Incr("d", 1), [ast.For(ast.Assign("c", 0), ast.Less("c", base_coords_dim), ast.Incr("c", 1), assign_block)]))
# Looks complicated, actually just produces the code below # void normal_calc (double** base_coords , double** normals ) # { # double v0[3] ; # double v1[3] ; # double n[3] ; # double x[3] ; # double dot ; # double norm ; # norm = 0.0; # dot = 0.0; # int i ; # int c ; # int d ; # for (i = 0; i < 3; ++i) # { # v0[i] = base_coords[1][i] - base_coords[0][i]; # v1[i] = base_coords[2][i] - base_coords[0][i]; # x[i] = 0.0; # } # n[0] = v0[1] * v1[2] - v0[2] * v1[1]; # n[1] = v0[2] * v1[0] - v0[0] * v1[2]; # n[2] = v0[0] * v1[1] - v0[1] * v1[0]; # for (i = 0; i < 3; ++i) # { # x[0] += base_coords[i][0]; # x[1] += base_coords[i][1]; # x[2] += base_coords[i][2]; # } # dot += (x[0]) * n[0]; # dot += (x[1]) * n[1]; # dot += (x[2]) * n[2]; # norm += n[0] * n[0]; # norm += n[1] * n[1]; # norm += n[2] * n[2]; # norm = sqrt(norm); # norm *= (dot < 0 ? -1 : 1); # for (d = 0; d < 1; ++d) # { # for (c = 0; c < 3; ++c) # { # normals[d][c] = n[c] / norm; # } # } # }
kernel = op2.Kernel(ast.FunDecl("void", "normal_calc", [ast.Decl("double**", "base_coords"), ast.Decl("double**", "normals")], ast.Block(body)), "normal_calc")
base_coords = mesh._old_mesh.coordinates op2.par_loop(kernel, zhat.cell_set, base_coords.dat(op2.READ, base_coords.cell_node_map()), zhat.dat(op2.WRITE, zhat.cell_node_map()))
### END NORMALS BLOCK ###
On 5 Dec 2014, at 10:01, Eike Mueller <E.Mueller@bath.ac.uk <mailto:E.Mueller@bath.ac.uk>> wrote:
Dear firedrakers,
I have the following function spaces:
W_2 = HDiv(U_2 x V_0) + HDiv(U_1 x V_1) [velocity] W_2^b = U_2 x V_0 [buoyancy]
where U_0 -> U_1 -> U_2 and V_0 -> V_1 are horizontal- and vertical de Rham complexes.
I need to be able to apply the operator Q: W_2^b -> W_2 defined by
Q_{ij} = <\vec{w}_i,\hat{z} \gamma_j>
where \vec{w}_i and \gamma_j are basis functions and \hat{z} is the unit normal in the vertical direction. To apply Q to a buoyancy field b, I want to write something like
u = assemble(dot(TestFunction(W_2),zhat*b)*dx)
but how do I get the zhat in there? Is that even possible?
Thanks a lot,
Eike
--
Dr Eike Hermann Mueller Research Associate (PostDoc)
Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom
+44 1225 38 5803 e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/> _______________________________________________ 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 5 Dec 2014, at 17:14, Eike Mueller <e.mueller@bath.ac.uk> wrote:
thanks, I guess I could also just use the code that Lawrence commented out below… I figured out what the C code does, and that’s what I’m doing now. In a radially extruded mesh, are the normals in one column always the same? I think this is what the code assumes, since it only uses the coordinates in the base cell.
Yes, the normal to the horizontal facet is always in the same direction, at least with no orography, since the horizontal facets in a column are parallel with radial extrusion Lawrence
Thanks, at some point we have to think about orography, but I think then we can just use the coordinates of the extruded mesh instead of the host mesh. Cheers, Eike Sent from my iPod
On 5 Dec 2014, at 17:28, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 5 Dec 2014, at 17:14, Eike Mueller <e.mueller@bath.ac.uk> wrote:
thanks, I guess I could also just use the code that Lawrence commented out below… I figured out what the C code does, and that’s what I’m doing now. In a radially extruded mesh, are the normals in one column always the same? I think this is what the code assumes, since it only uses the coordinates in the base cell.
Yes, the normal to the horizontal facet is always in the same direction, at least with no orography, since the horizontal facets in a column are parallel with radial extrusion
Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
I don't think orography changes anything, but modelling Earth as an oblate spheroid would. (Had to teach my smart phone a bunch of new words there!) Cjc On 5 Dec 2014 17:33, "Eike Mueller" <e.mueller@bath.ac.uk> wrote:
Thanks, at some point we have to think about orography, but I think then we can just use the coordinates of the extruded mesh instead of the host mesh.
Cheers, Eike
Sent from my iPod
On 5 Dec 2014, at 17:28, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 5 Dec 2014, at 17:14, Eike Mueller <e.mueller@bath.ac.uk> wrote:
thanks, I guess I could also just use the code that Lawrence commented out below… I figured out what the C code does, and that’s what I’m doing now. In a radially extruded mesh, are the normals in one column always the same? I think this is what the code assumes, since it only uses the coordinates in the base cell.
Yes, the normal to the horizontal facet is always in the same direction, at least with no orography, since the horizontal facets in a column are parallel with radial extrusion
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
participants (4)
-
Colin Cotter
-
Eike Mueller
-
Fabio Luporini
-
Lawrence Mitchell