Dear Firedrake -- I have two questions related to mesh coordinates and boundaries: 1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors? 2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this? Thanks for your help with my questions so far! Ed -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
Hi Ed, I think that the answer to both of your questions is: x = SpatialCoordinate(mesh). x is then a symbolic representation of the coordinate vector. You can use x in any UFL expression, including for boundary conditions. In the case of your second question, you could use. It in either an interpolation or a projection: W = interpolate(w(x), fs) or W = project(w(x), fs) Where w(x) is a UFL expression in x (or a python function returning such an expression) and fs is the FunctionSpace you want to interpolate into. Regards, David Dr David Ham Department of Mathematics Imperial College London
On 10 May 2018, at 07:55, Ed Bueler <elbueler@alaska.edu> wrote:
Dear Firedrake --
I have two questions related to mesh coordinates and boundaries:
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
Thanks for your help with my questions so far!
Ed
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Here's an attempt to partially answer, and follow up on, my own questions. These were real questions a couple of days ago, but now I see that they are avoidable and don't really need answers.
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer. My much-reduced question has become: How to generate the indicator function of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g. the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec. My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way. In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine. It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier. Ed On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu> wrote:
Dear Firedrake --
I have two questions related to mesh coordinates and boundaries:
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
Thanks for your help with my questions so far!
Ed
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
Dear Ed, Here’s an attempt to lay out how boundary conditions work, which might enable you to do what you’re after. Let’s make some basic data structures: from firedrake import * m = UnitSquareMesh(2,2) fs = FunctionSpace(m, "CG", 2) bc = DirichletBC(fs, 1.0, 1) Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function. The (process-local) indices of the nodes in bc are given by: print(bc.nodes) array([ 0, 3, 4, 12, 14], dtype=int32) So let’s suppose we want that indicator function: f = Function(fs) bc.apply(f) If we actually look at the node values, we now find: print(f.dat.data_ro) array([1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) Note that the fact that the nodal values are 1 is a function of the node definition. It will be true for CG and DG spaces, because the node definition for these cases is point evaluation. For more esoteric finite elements such as those in H(div) and H(curl), the relationship between nodal values and Function value are more complex. If you want the physical locations of these nodes then you’ll need to interpolate the mesh coordinates into the VectorFunctionSpace corresponding to the FunctionSpace you’re working in. This process is described here: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter... In order to only perform the operation over the boundary nodes, you’d simply replace the last line there with: f.dat.data[bc.nodes] = mydata(X.dat.data_ro[bc.nodes]) Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of Ed Bueler <elbueler@alaska.edu> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Saturday, 12 May 2018 at 01:36 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] mesh coordinates and boundaries questions Here's an attempt to partially answer, and follow up on, my own questions. These were real questions a couple of days ago, but now I see that they are avoidable and don't really need answers.
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer. My much-reduced question has become: How to generate the indicator function of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g. the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec. My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way. In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine. It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier. Ed On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu<mailto:elbueler@alaska.edu>> wrote: Dear Firedrake -- I have two questions related to mesh coordinates and boundaries: 1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors? 2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this? Thanks for your help with my questions so far! Ed -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
David --
... Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function. The (process-local) indices of the nodes in bc are given by: print(bc.nodes) array([ 0, 3, 4, 12, 14], dtype=int32) So let’s suppose we want that indicator function: f = Function(fs) bc.apply(f) ...
Got it. It was the "bc.apply()" which I was missing. It is much more direct than the approach I eventually found in https://www.firedrakeprojec t.org/demos/linear_fluid_structure_interaction.py.html which allowed me to use par_loop() for the purpose.
I think my latest missive missed 2 below. You are right. The mechanism we provide will ask for data on each process. If your data is also distributed then that won’t necessarily be the same process you have data on. I think you just have to manage the communication in that case. It’s not clear to me how we could fix that one.
I don't know how to "fix" that in general either. But now I am not sure I need what I asked for ... Thanks! Ed On Sun, May 13, 2018 at 3:52 AM, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Ed,
Here’s an attempt to lay out how boundary conditions work, which might enable you to do what you’re after.
Let’s make some basic data structures:
*from* *firedrake* *import* *
m = UnitSquareMesh(2,2)
fs = FunctionSpace(m, "CG", 2)
bc = DirichletBC(fs, 1.0, 1)
Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function.
The (process-local) indices of the nodes in bc are given by:
print(bc.nodes)
array([ 0, 3, 4, 12, 14], dtype=int32)
So let’s suppose we want that indicator function:
f = Function(fs)
bc.apply(f)
If we actually look at the node values, we now find:
print(f.dat.data_ro)
array([1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0.])
Note that the fact that the nodal values are 1 is a function of the node definition. It will be true for CG and DG spaces, because the node definition for these cases is point evaluation. For more esoteric finite elements such as those in H(div) and H(curl), the relationship between nodal values and Function value are more complex.
If you want the physical locations of these nodes then you’ll need to interpolate the mesh coordinates into the VectorFunctionSpace corresponding to the FunctionSpace you’re working in. This process is described here: https://www.firedrakeproject.org/interpolation.html# interpolation-from-external-data
In order to only perform the operation over the boundary nodes, you’d simply replace the last line there with:
f.dat.data[bc.nodes] = mydata(X.dat.data_ro[bc.nodes])
Regards,
David
*From: *<firedrake-bounces@imperial.ac.uk> on behalf of Ed Bueler < elbueler@alaska.edu> *Reply-To: *firedrake <firedrake@imperial.ac.uk> *Date: *Saturday, 12 May 2018 at 01:36 *To: *firedrake <firedrake@imperial.ac.uk> *Subject: *Re: [firedrake] mesh coordinates and boundaries questions
Here's an attempt to partially answer, and follow up on, my own questions. These were real questions a couple of days ago, but now I see that they are avoidable and don't really need answers.
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id
number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form.
What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)"
(though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer.
My much-reduced question has become: How to generate the indicator function of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
2) Given a function along a boundary, how does one extend that to the whole mesh,
by holding constant in one variable? In the case I care about the domain is of the special form
D = {(x,y) | f(x) < y < g(x)}
and one has values w(x) along the y=g(x) boundary. How does one build a Function on the
whole mesh which is
W(x,y) = w(x)
(I understand that my domain is special, namely logically rectangular. And I can imagine an
answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested
in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g. the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec.
My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way.
In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine. It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier.
Ed
On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu> wrote:
Dear Firedrake --
I have two questions related to mesh coordinates and boundaries:
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form
D = {(x,y) | f(x) < y < g(x)}
and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is
W(x,y) = w(x)
(I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
Thanks for your help with my questions so far!
Ed
--
Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
--
Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
David -- I've made lots of progress and kargekt have the code I want ... except that it only runs in serial, which is why this follow-up. Context: I'm solving Stokes with a free surface z = h(x,t) which evolves according to a surface kinematical equation. I am (currently) doing this explicitly by stretching the mesh in the vertical direction. This now works as a conditionally-stable scheme in serial. The scheme depends on the following function: # return linear-interpolated surface profile function h(x) def getsurfaceprofile(mesh,top_id): from scipy import interpolate P1 = FunctionSpace(mesh, "CG", 1) bc = DirichletBC(P1, 1.0, top_id) # notes: 1) bc.nodes gives indices to mesh; is a 1D dtype=int32 numpy array # 2) f = Function(P1); bc.apply(f) gives indicator function x,z = SpatialCoordinate(mesh) xh = Function(P1).interpolate(x).dat.data[bc.nodes] # 1D numpy array zh = Function(P1).interpolate(z).dat.data[bc.nodes] return interpolate.interp1d(xh,zh,copy=False) # default is 'linear', which is what we want This function returns a piecewise-linear function which I then use on the x-coordinates in the rest of the mesh. But it will not work in parallel both because I have not (yet) handled the halo issue in the free surface *and* because the mesh is not partitioned in parallel according to x. For now I want to avoid extruded meshes, which could be partitioned only in the x direction, though I acknowledge that may be where I am forced to go. So basically the questions are: Any suggestions about redoing the above idea in parallel? Is there a firedrake tool I am missing? Would firedrake.projection.project() help? Thanks! I really appreciate what firedrake has done! Beautiful super-petsc layer. Ed On Mon, May 14, 2018 at 9:17 AM, Ed Bueler <elbueler@alaska.edu> wrote:
David --
... Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function. The (process-local) indices of the nodes in bc are given by: print(bc.nodes) array([ 0, 3, 4, 12, 14], dtype=int32) So let’s suppose we want that indicator function: f = Function(fs) bc.apply(f) ...
Got it. It was the "bc.apply()" which I was missing. It is much more direct than the approach I eventually found in https://www.firedrakeprojec t.org/demos/linear_fluid_structure_interaction.py.html which allowed me to use par_loop() for the purpose.
I think my latest missive missed 2 below. You are right. The mechanism we provide will ask for data on each process. If your data is also distributed then that won’t necessarily be the same process you have data on. I think you just have to manage the communication in that case. It’s not clear to me how we could fix that one.
I don't know how to "fix" that in general either. But now I am not sure I need what I asked for ...
Thanks!
Ed
On Sun, May 13, 2018 at 3:52 AM, Ham, David A <david.ham@imperial.ac.uk> wrote:
Dear Ed,
Here’s an attempt to lay out how boundary conditions work, which might enable you to do what you’re after.
Let’s make some basic data structures:
*from* *firedrake* *import* *
m = UnitSquareMesh(2,2)
fs = FunctionSpace(m, "CG", 2)
bc = DirichletBC(fs, 1.0, 1)
Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function.
The (process-local) indices of the nodes in bc are given by:
print(bc.nodes)
array([ 0, 3, 4, 12, 14], dtype=int32)
So let’s suppose we want that indicator function:
f = Function(fs)
bc.apply(f)
If we actually look at the node values, we now find:
print(f.dat.data_ro)
array([1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0.])
Note that the fact that the nodal values are 1 is a function of the node definition. It will be true for CG and DG spaces, because the node definition for these cases is point evaluation. For more esoteric finite elements such as those in H(div) and H(curl), the relationship between nodal values and Function value are more complex.
If you want the physical locations of these nodes then you’ll need to interpolate the mesh coordinates into the VectorFunctionSpace corresponding to the FunctionSpace you’re working in. This process is described here: https://www.firedrakeproject.org/interpolation.html#interpol ation-from-external-data
In order to only perform the operation over the boundary nodes, you’d simply replace the last line there with:
f.dat.data[bc.nodes] = mydata(X.dat.data_ro[bc.nodes])
Regards,
David
*From: *<firedrake-bounces@imperial.ac.uk> on behalf of Ed Bueler < elbueler@alaska.edu> *Reply-To: *firedrake <firedrake@imperial.ac.uk> *Date: *Saturday, 12 May 2018 at 01:36 *To: *firedrake <firedrake@imperial.ac.uk> *Subject: *Re: [firedrake] mesh coordinates and boundaries questions
Here's an attempt to partially answer, and follow up on, my own questions. These were real questions a couple of days ago, but now I see that they are avoidable and don't really need answers.
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id
number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form.
What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)"
(though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer.
My much-reduced question has become: How to generate the indicator function of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
2) Given a function along a boundary, how does one extend that to the whole mesh,
by holding constant in one variable? In the case I care about the domain is of the special form
D = {(x,y) | f(x) < y < g(x)}
and one has values w(x) along the y=g(x) boundary. How does one build a Function on the
whole mesh which is
W(x,y) = w(x)
(I understand that my domain is special, namely logically rectangular. And I can imagine an
answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested
in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g. the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec.
My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way.
In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine. It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier.
Ed
On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu> wrote:
Dear Firedrake --
I have two questions related to mesh coordinates and boundaries:
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form
D = {(x,y) | f(x) < y < g(x)}
and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is
W(x,y) = w(x)
(I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
Thanks for your help with my questions so far!
Ed
--
Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
--
Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
-- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
Hi Ed, If I understand correctly, you’re trying to solve a 1D PDE on a surface mesh. As I’m sure you’ve noticed, Firedrake doesn’t really support doing that (so you’ve got a hack which works in serial). This is a major wish-list item which hopefully the new Firedrake post-doc will work on once we get through the hiring process, but it’s not there now. As for the second part of what you want, which if I understand it correctly is to propagate the surface values through the domain, this can be done right now by solving the equation partial h/partial z = 0 with the surface value as boundary condition. Because that’s just a PDE, Firedrake’s usual parallel support works. However, all of the above comes with one huge caveat. Firedrake’s strong scaling performance only goes down to several tens of thousands of degrees of freedom per core. That’s significantly worse than world class but better than many codes out there. The practical upshot of this is that many 2D problems aren’t going to significantly benefit from parallelism. Unless you have hundreds of thousands of degrees of freedom, parallel is probably not worth it. (In 3D, of course, you hit hundreds of thousands of DoFs really fast). Regards, David From: Ed Bueler <elbueler@alaska.edu> Date: Saturday, 19 May 2018 at 01:34 To: "Ham, David A" <david.ham@imperial.ac.uk> Cc: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] mesh coordinates and boundaries questions David -- I've made lots of progress and kargekt have the code I want ... except that it only runs in serial, which is why this follow-up. Context: I'm solving Stokes with a free surface z = h(x,t) which evolves according to a surface kinematical equation. I am (currently) doing this explicitly by stretching the mesh in the vertical direction. This now works as a conditionally-stable scheme in serial. The scheme depends on the following function: # return linear-interpolated surface profile function h(x) def getsurfaceprofile(mesh,top_id): from scipy import interpolate P1 = FunctionSpace(mesh, "CG", 1) bc = DirichletBC(P1, 1.0, top_id) # notes: 1) bc.nodes gives indices to mesh; is a 1D dtype=int32 numpy array # 2) f = Function(P1); bc.apply(f) gives indicator function x,z = SpatialCoordinate(mesh) xh = Function(P1).interpolate(x).dat.data[bc.nodes] # 1D numpy array zh = Function(P1).interpolate(z).dat.data[bc.nodes] return interpolate.interp1d(xh,zh,copy=False) # default is 'linear', which is what we want This function returns a piecewise-linear function which I then use on the x-coordinates in the rest of the mesh. But it will not work in parallel both because I have not (yet) handled the halo issue in the free surface *and* because the mesh is not partitioned in parallel according to x. For now I want to avoid extruded meshes, which could be partitioned only in the x direction, though I acknowledge that may be where I am forced to go. So basically the questions are: Any suggestions about redoing the above idea in parallel? Is there a firedrake tool I am missing? Would firedrake.projection.project() help? Thanks! I really appreciate what firedrake has done! Beautiful super-petsc layer. Ed On Mon, May 14, 2018 at 9:17 AM, Ed Bueler <elbueler@alaska.edu<mailto:elbueler@alaska.edu>> wrote: David --
... Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function. The (process-local) indices of the nodes in bc are given by: print(bc.nodes) array([ 0, 3, 4, 12, 14], dtype=int32) So let’s suppose we want that indicator function: f = Function(fs) bc.apply(f) ...
Got it. It was the "bc.apply()" which I was missing. It is much more direct than the approach I eventually found in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... which allowed me to use par_loop() for the purpose.
I think my latest missive missed 2 below. You are right. The mechanism we provide will ask for data on each process. If your data is also distributed then that won’t necessarily be the same process you have data on. I think you just have to manage the communication in that case. It’s not clear to me how we could fix that one.
I don't know how to "fix" that in general either. But now I am not sure I need what I asked for ... Thanks! Ed On Sun, May 13, 2018 at 3:52 AM, Ham, David A <david.ham@imperial.ac.uk<mailto:david.ham@imperial.ac.uk>> wrote: Dear Ed, Here’s an attempt to lay out how boundary conditions work, which might enable you to do what you’re after. Let’s make some basic data structures: from firedrake import * m = UnitSquareMesh(2,2) fs = FunctionSpace(m, "CG", 2) bc = DirichletBC(fs, 1.0, 1) Notice I’ve set up a DirichletBC with the constant value 1.0. This can be used to create an indicator function. The (process-local) indices of the nodes in bc are given by: print(bc.nodes) array([ 0, 3, 4, 12, 14], dtype=int32) So let’s suppose we want that indicator function: f = Function(fs) bc.apply(f) If we actually look at the node values, we now find: print(f.dat.data_ro) array([1., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) Note that the fact that the nodal values are 1 is a function of the node definition. It will be true for CG and DG spaces, because the node definition for these cases is point evaluation. For more esoteric finite elements such as those in H(div) and H(curl), the relationship between nodal values and Function value are more complex. If you want the physical locations of these nodes then you’ll need to interpolate the mesh coordinates into the VectorFunctionSpace corresponding to the FunctionSpace you’re working in. This process is described here: https://www.firedrakeproject.org/interpolation.html#interpolation-from-exter... In order to only perform the operation over the boundary nodes, you’d simply replace the last line there with: f.dat.data[bc.nodes] = mydata(X.dat.data_ro[bc.nodes]) Regards, David From: <firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>> on behalf of Ed Bueler <elbueler@alaska.edu<mailto:elbueler@alaska.edu>> Reply-To: firedrake <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Date: Saturday, 12 May 2018 at 01:36 To: firedrake <firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk>> Subject: Re: [firedrake] mesh coordinates and boundaries questions Here's an attempt to partially answer, and follow up on, my own questions. These were real questions a couple of days ago, but now I see that they are avoidable and don't really need answers.
1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors?
Generating a Function or PETSc Vec defined only on the desired part of the boundary, and not defined elsewhere, would generally imply that some processes had no values, which can be awkward. In any case I don't need this any longer. My much-reduced question has become: How to generate the indicator function of part of the boundary? This can be locally-generated and is distributed the same way as a Function on the whole mesh. The FEniCS book shows an example of generating MeshFunctions which indicate boundary parts, but perhaps this is just Function in Firedrake? What is the idiom for using the subdomain id to generate the indicator function?
2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this?
I was missing something deeper here. In parallel each process only owns (knows about) its part of the mesh. Other than solve(), Firedrake operations are almost all local in the sense that each process works on its part of the mesh, or across a halo, but in any case without global knowledge. By contrast, solve() is a well-defined interface for spreading certain effects from some locales (e.g. the boundary) to all processes. One needs precise equations to solve ... but such equations give precise meaning to "spreading certain effects". Thus solve() is the key to global influence on the entries of a Function and/or its underlying mpi-type Vec. My question could therefore answered this way: Set up a weak form to solve a differential equation problem describing the field you want with the boundary values you want. Then call solve(). This approach can propagate information to all processes in a well-supported way. In fact it turns out that what I wanted could be done by setting up and solving a scalar, linear, 2nd-order PDE problem. It works fine. It is computationally-cheap compared to the "main" PDE I was solving, and scales better. In the longer term I almost certainly will use an extruded mesh, which may be even easier. Ed On Wed, May 9, 2018 at 10:54 PM, Ed Bueler <elbueler@alaska.edu<mailto:elbueler@alaska.edu>> wrote: Dear Firedrake -- I have two questions related to mesh coordinates and boundaries: 1) How does one get the coordinates of the points on a boundary? That is, suppose I know the id number of that part, for example "id" as in "ds(id)" used in forming boundary terms in the weak form. What is the corresponding idiom for getting the spatial coordinates? Something like "mesh.coordinates(id)" (though that does not work)? In parallel I guess this function will only be defined on a subset of the processors? 2) Given a function along a boundary, how does one extend that to the whole mesh, by holding constant in one variable? In the case I care about the domain is of the special form D = {(x,y) | f(x) < y < g(x)} and one has values w(x) along the y=g(x) boundary. How does one build a Function on the whole mesh which is W(x,y) = w(x) (I understand that my domain is special, namely logically rectangular. And I can imagine an answer specific to an extruded mesh ... but I would like to not go there now.) Again I am interested in a solution which works in parallel. Maybe I have to commit to solving a first-order PDE just to do this? Thanks for your help with my questions so far! Ed -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman -- Ed Bueler Dept of Math and Stat and Geophysical Institute University of Alaska Fairbanks Fairbanks, AK 99775-6660 301C Chapman
On 19/05/18 01:34, Ed Bueler wrote:
David --
I've made lots of progress and kargekt have the code I want ... except that it only runs in serial, which is why this follow-up. Context: I'm solving Stokes with a free surface z = h(x,t) which evolves according to a surface kinematical equation. I am (currently) doing this explicitly by stretching the mesh in the vertical direction. This now works as a conditionally-stable scheme in serial.
The scheme depends on the following function:
# return linear-interpolated surface profile function h(x) def getsurfaceprofile(mesh,top_id): from scipy import interpolate P1 = FunctionSpace(mesh, "CG", 1) bc = DirichletBC(P1, 1.0, top_id) # notes: 1) bc.nodes gives indices to mesh; is a 1D dtype=int32 numpy array # 2) f = Function(P1); bc.apply(f) gives indicator function x,z = SpatialCoordinate(mesh) xh = Function(P1).interpolate(x).dat.data[bc.nodes] # 1D numpy array zh = Function(P1).interpolate(z).dat.data[bc.nodes] return interpolate.interp1d(xh,zh,copy=False) # default is 'linear', which is what we want
This function returns a piecewise-linear function which I then use on the x-coordinates in the rest of the mesh. But it will not work in parallel both because I have not (yet) handled the halo issue in the free surface *and* because the mesh is not partitioned in parallel according to x.
I think (although you don't say) that one of this issues in parallel is likely to be that bc.nodes reads of the end of the data array. This is because we have not been careful with the naming: bc.nodes contains ghost nodes, where as `.dat.data` is a view only of the local dofs. So you need to say `.dat.data_with_halos`. But as David notes, there are other potential caveats. Cheers, Lawrence
participants (3)
-
Ed Bueler
-
Ham, David A
-
Lawrence Mitchell