Parallel issue with Subroutine: Vertical_Extrapolation.F90

[fangf]

There is an issue with parallel computing when using the subroutine: Vertical_Extrapolation.F90

Let take the test case (water_collapse_2d) as an example. The error message is: FLUIDITY ERROR Source location: (Vertical_Extrapolation.F90, 431) Error message: Something wrong with the geometry. Backtrace will follow if it is available: Initialising picker for field DistanceToTopHorizontalCoordinate New picker ID: 1 Updating halo CoordinateMeshLevel2Halo for field DistanceToTop

in code: do i=1, size(seles) if (seles(i)<0) then ewrite(-1,) "For node with coordinate", node_val(mesh_positions, i) ewrite(-1,) "no top surface node was found." FLAbort("Something wrong with the geometry.") else seles(i) = horizontal_mesh_list(seles(i)) end if end do

[stephankramer]

That's a very timely issue!

Yes, vertical extrapolation has never worked in parallel for meshes that are not extruded. This is any setup with /geometry/ocean_boundaries/ which is required e.g. for cases with both free surface and mesh movement. In the case of water_collapse_2d I'm not actually sure why it has /geometry/ocean_boundaries, so I think you can just remove it.

The issue is that vertical extrapolation requires interpolation from the top and bottom surface mesh to the interior mesh, but in parallel the surface mesh will be decomposed completely differently from the interior mesh, i.e. a local domain somewhere in the interior will not even be aware where the surface mesh (above and below it) is or which other processor owns it. So it's not just a case of exchanging halo data with neighbouring processes. If you use an extruded mesh, then the decomposition is in the horizontal, so each local process still has the relevant bit of top and bottom surface mesh on top and beneath its local part of the mesh - of course this means that you can't use general adaptivity (only 2+1 adaptivity) either.

If you do need to use vertical extrapolation in parallel on a non-extruded mesh, you're in luck, because I have just yesterday put in a PR (#334 ) that implements it. Just to warn, it requires quite a lot of communication (between owners of the domain decomposed surface mesh and all other processes) so might not scale very well on a large number of cores, but it's probably fine for medium sized problems.

[jrper]

I've got a feeling that that particular example uses (or used) the DisfanceToTop field to calculate and subtract hydrostatic components of pressure, so may not function as intended with the /geometry/ocean_boundaries/ option removed, but it's a long time since I last looked at it.

[stephankramer]

Fixed in #334