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