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