Dear Firedrakers, I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface. I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match. Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this? Cheers, Daniel -- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht@leeds.ac.uk Phone: +44 (0)113 343 22 01
Dear Daniel,
On 30 May 2016, at 13:29, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Dear Firedrakers,
I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface.
I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match.
This will work, but as you note, is very slow. I presume the identification of the degrees of freedom on the boundary that you need to move from one mesh to the other is static (i.e. you only need to do so once). If that is the case you could do this up front and then (assuming you solve in A, transfer to B, solve in B) you could do something like: solve in A boundary_function_in_B.dat.data[boundary_nodes_B] = solution_in_A.dat.data_ro[boundary_nodes_in_A] solve in B
Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this?
How are you describing the coupled problem. As two completely separate meshes, or as a single mesh in which you mask out the equations in one half or the other? Moreover, is the coupling of the form: solve in domain A transfer solution to domain B solve in domain B using the solution from A as a boundary condition Or do you solve the full coupled A-B problem in one go. The current approach in firedrake for solving in subdomains is currently a bit hacky, and what I describe now only works if you have a single mesh. What you can do is as follows. Mark the two regions of your domain (firedrake as of earlier last week correct reads cell subdomains from Gmsh files). Now you can define an indicator function in a DG0 functionspace I_A = Function(DG0) I_B = Function(DG0) which you can initialise using a par_loop to be 1 in the appropriate part of the domain. You can use these indicator functions in your forms to only select contributions from the appropriate domain in the appropriate parts. e.g. you can use: dx_A = I_A*dx dx_B = I_B*dx ds_A = I_A*ds ds_B = I_B*ds To apply the neumann boundary term you'll need to use an interior facet integral. This requires a little more thought as to how to restrict the indicator functions, but is not too bad. You'll then want to mask out the appropriate parts of the operators in the two domains with strong boundary conditions in the "other" half. You can do this by subclassing DirichletBC and overriding the "nodes" property to select all function space nodes that are in the wrong domain. A better approach to this latter problem is on our roadmap (although we do not currently have a timescale for it). See https://github.com/firedrakeproject/firedrake/issues/723 for details. Please comment if you have ideas on what should be available. Note that coupling two *unrelated* meshes is a harder problem https://github.com/firedrakeproject/firedrake/issues/714 effectively because the operator mapping between meshes (even if they are topologically conforming) is almost always not the identity. Hope this is helpful! Lawrence
Hi Lawrence, Thanks for the info - I do have two independent meshes (the first variant you describe), but the problem is not fully coupled: so yes, the sequence would be for every time step 1. Solve Problem A 2. Get Boundary Data of Problem A 3. Prescribe Neumann Data for Problem B based on Boundary Data from A (*) 4. Solve B (*) We do not need to actually couple the meshes in terms of identifying DOFs or so - it would be enough to some how create an expression from the boundary data of A and feed this into the boundary condition for B. So at the moment the idea is to derive from Expression and add an additional argument which can be a Function on mesh A. Then, in eval, I access values on points of mesh B through .at(X, dont_raise=True). I can't use a single mesh because the mesh of problem A will move in time - only in a simple way, so we can always identify where the two meshes would touch, but this still prevents us from doing everything on a single mesh. Cheers, Daniel On 30/05/16 15:13, Lawrence Mitchell wrote:
Dear Daniel,
On 30 May 2016, at 13:29, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Dear Firedrakers,
I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface.
I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match. This will work, but as you note, is very slow. I presume the identification of the degrees of freedom on the boundary that you need to move from one mesh to the other is static (i.e. you only need to do so once). If that is the case you could do this up front and then (assuming you solve in A, transfer to B, solve in B) you could do something like:
solve in A
boundary_function_in_B.dat.data[boundary_nodes_B] = solution_in_A.dat.data_ro[boundary_nodes_in_A]
solve in B
Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this? How are you describing the coupled problem. As two completely separate meshes, or as a single mesh in which you mask out the equations in one half or the other?
Moreover, is the coupling of the form:
solve in domain A
transfer solution to domain B
solve in domain B using the solution from A as a boundary condition
Or do you solve the full coupled A-B problem in one go.
The current approach in firedrake for solving in subdomains is currently a bit hacky, and what I describe now only works if you have a single mesh.
What you can do is as follows. Mark the two regions of your domain (firedrake as of earlier last week correct reads cell subdomains from Gmsh files).
Now you can define an indicator function in a DG0 functionspace
I_A = Function(DG0) I_B = Function(DG0)
which you can initialise using a par_loop to be 1 in the appropriate part of the domain.
You can use these indicator functions in your forms to only select contributions from the appropriate domain in the appropriate parts.
e.g. you can use:
dx_A = I_A*dx dx_B = I_B*dx ds_A = I_A*ds ds_B = I_B*ds
To apply the neumann boundary term you'll need to use an interior facet integral. This requires a little more thought as to how to restrict the indicator functions, but is not too bad.
You'll then want to mask out the appropriate parts of the operators in the two domains with strong boundary conditions in the "other" half. You can do this by subclassing DirichletBC and overriding the "nodes" property to select all function space nodes that are in the wrong domain.
A better approach to this latter problem is on our roadmap (although we do not currently have a timescale for it). See https://github.com/firedrakeproject/firedrake/issues/723 for details. Please comment if you have ideas on what should be available.
Note that coupling two *unrelated* meshes is a harder problem https://github.com/firedrakeproject/firedrake/issues/714 effectively because the operator mapping between meshes (even if they are topologically conforming) is almost always not the identity.
Hope this is helpful!
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht@leeds.ac.uk Phone: +44 (0)113 343 22 01
On 30 May 2016, at 15:24, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Hi Lawrence,
Thanks for the info - I do have two independent meshes (the first variant you describe), but the problem is not fully coupled: so yes, the sequence would be for every time step
1. Solve Problem A 2. Get Boundary Data of Problem A 3. Prescribe Neumann Data for Problem B based on Boundary Data from A (*) 4. Solve B
(*) We do not need to actually couple the meshes in terms of identifying DOFs or so - it would be enough to some how create an expression from the boundary data of A and feed this into the boundary condition for B.
So at the moment the idea is to derive from Expression and add an additional argument which can be a Function on mesh A. Then, in eval, I access values on points of mesh B through .at(X, dont_raise=True).
OK, so you effectively are evaluating (a function of) function_on_A at the nodes points of mesh B that live on the boundary? I would probably do this like so: meshA = ... meshB = ... V_a = FunctionSpace(meshA, element_a) V_b = FunctionSpace(meshB, "CG", 1) solution_A = Function(V_a) # Assume that meshB has boundary marker 1 on the appropriate boundary boundary_nodes = DirichletBC(V_b, 0, 1).nodes # Get the global coordinates of the boundary nodes by interpolating the # coordinate field into the vector-valued version of the space V_b c = Function(Vector_V_b) c.interpolate(SpatialCoordinate(meshB)) Xs = c.dat.data_ro[boundary_nodes, :] boundary_data_B = Function(V_b) # solve problem A ... # transfer to domain B. # You'll probably want to filter out the `None`s. boundary_data_B.dat.data[boundary_nodes] = solution_A.at(Xs, dont_raise=True) # solve problem B ... Note that .at is currently collective, so this won't work right in parallel. This is sort of what eval in the Expression is doing, but restricted to the boundary, and "vectorised" so it should be a little faster. Lawrence
I can't use a single mesh because the mesh of problem A will move in time - only in a simple way, so we can always identify where the two meshes would touch, but this still prevents us from doing everything on a single mesh.
Hi Lawrence, We set up the problem according to your info and it works initially. However, as we move meshA, at the some point all the results from solution_A.at(Xs, dont_raise=True) become None's and we have absolutely no idea why this is happening - it works for about 5 time steps and then it doesn't anymore. I attached a small example reproducing the problem and a short video showing what is happening. To illustrate, we print out solution_A.at((0.0, 2.2), dont_raise=True) and the point is marked in the video. For the first three steps, the result is the correct 5.0 (an arbitrary value we assign to the solution on meshA at the moment) but after that it becomes None, even though the point is still clearly ``in contact'' with meshA. Any ideas what is happening here? Cheers, Daniel On 30/05/16 16:04, Lawrence Mitchell wrote:
On 30 May 2016, at 15:24, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Hi Lawrence,
Thanks for the info - I do have two independent meshes (the first variant you describe), but the problem is not fully coupled: so yes, the sequence would be for every time step
1. Solve Problem A 2. Get Boundary Data of Problem A 3. Prescribe Neumann Data for Problem B based on Boundary Data from A (*) 4. Solve B
(*) We do not need to actually couple the meshes in terms of identifying DOFs or so - it would be enough to some how create an expression from the boundary data of A and feed this into the boundary condition for B.
So at the moment the idea is to derive from Expression and add an additional argument which can be a Function on mesh A. Then, in eval, I access values on points of mesh B through .at(X, dont_raise=True). OK, so you effectively are evaluating (a function of) function_on_A at the nodes points of mesh B that live on the boundary?
I would probably do this like so:
meshA = ... meshB = ...
V_a = FunctionSpace(meshA, element_a)
V_b = FunctionSpace(meshB, "CG", 1)
solution_A = Function(V_a)
# Assume that meshB has boundary marker 1 on the appropriate boundary
boundary_nodes = DirichletBC(V_b, 0, 1).nodes
# Get the global coordinates of the boundary nodes by interpolating the # coordinate field into the vector-valued version of the space V_b
c = Function(Vector_V_b) c.interpolate(SpatialCoordinate(meshB))
Xs = c.dat.data_ro[boundary_nodes, :]
boundary_data_B = Function(V_b)
# solve problem A
... # transfer to domain B. # You'll probably want to filter out the `None`s. boundary_data_B.dat.data[boundary_nodes] = solution_A.at(Xs, dont_raise=True)
# solve problem B ...
Note that .at is currently collective, so this won't work right in parallel.
This is sort of what eval in the Expression is doing, but restricted to the boundary, and "vectorised" so it should be a little faster.
Lawrence
I can't use a single mesh because the mesh of problem A will move in time - only in a simple way, so we can always identify where the two meshes would touch, but this still prevents us from doing everything on a single mesh.
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht@leeds.ac.uk Phone: +44 (0)113 343 22 01
On 31/05/16 10:08, Daniel Ruprecht wrote:
Hi Lawrence,
We set up the problem according to your info and it works initially. However, as we move meshA, at the some point all the results from
solution_A.at(Xs, dont_raise=True)
become None's and we have absolutely no idea why this is happening - it works for about 5 time steps and then it doesn't anymore. I attached a small example reproducing the problem and a short video showing what is happening.
To illustrate, we print out
solution_A.at((0.0, 2.2), dont_raise=True)
and the point is marked in the video. For the first three steps, the result is the correct 5.0 (an arbitrary value we assign to the solution on meshA at the moment) but after that it becomes None, even though the point is still clearly ``in contact'' with meshA.
Any ideas what is happening here?
Yes. The point evaluation stuff doesn't currently deal with moving meshes. We build an rtree for point queries, but only do that once. So after a while, once you've moved the mesh the bounding boxes are in the wrong place. To fix this, every time you move the mesh, you can do: del meshA.spatial_index If you report a bug, we will track it and add a public API for clearing the rtree. A few comments on your code. Please, where possible, use interpolation of UFL expressions (rather than C strings) as described here: http://firedrakeproject.org/interpolation.html#ufl-expressions This gets you much better type checking. Additionally, there's no need to do the WithGeometry dance for u_A, you're not changing the mesh's coordinate function space, only its values. Cheers, Lawrence
Hi Lawrence, That did the trick, thanks a lot! Also thanks for your comments on the code. I will open an issue on GitHub. Cheers, Daniel On 31/05/16 11:06, Lawrence Mitchell wrote:
On 31/05/16 10:08, Daniel Ruprecht wrote:
Hi Lawrence,
We set up the problem according to your info and it works initially. However, as we move meshA, at the some point all the results from
solution_A.at(Xs, dont_raise=True)
become None's and we have absolutely no idea why this is happening - it works for about 5 time steps and then it doesn't anymore. I attached a small example reproducing the problem and a short video showing what is happening.
To illustrate, we print out
solution_A.at((0.0, 2.2), dont_raise=True)
and the point is marked in the video. For the first three steps, the result is the correct 5.0 (an arbitrary value we assign to the solution on meshA at the moment) but after that it becomes None, even though the point is still clearly ``in contact'' with meshA.
Any ideas what is happening here? Yes. The point evaluation stuff doesn't currently deal with moving meshes. We build an rtree for point queries, but only do that once. So after a while, once you've moved the mesh the bounding boxes are in the wrong place.
To fix this, every time you move the mesh, you can do:
del meshA.spatial_index
If you report a bug, we will track it and add a public API for clearing the rtree.
A few comments on your code.
Please, where possible, use interpolation of UFL expressions (rather than C strings) as described here: http://firedrakeproject.org/interpolation.html#ufl-expressions
This gets you much better type checking.
Additionally, there's no need to do the WithGeometry dance for u_A, you're not changing the mesh's coordinate function space, only its values.
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht@leeds.ac.uk Phone: +44 (0)113 343 22 01
participants (2)
-
Daniel Ruprecht
-
Lawrence Mitchell