[vtk-developers] SPH Interpolator

Sat Apr 23 16:24:29 EDT 2016

Will,

>So far we have implemented cubic, quartic and quintic SPH kernels.

There are more kernels you can easily adapt in the pv-meshless repo (they have been tested extensively over the years)
https://github.com/biddisco/pv-meshless/tree/master/sph

(the sph interpolation classes in there handle a variable kernel size per point)

JB

From: "vtk-developers-bounces at vtk.org<mailto:vtk-developers-bounces at vtk.org>" <vtk-developers-bounces at vtk.org<mailto:vtk-developers-bounces at vtk.org>> on behalf of Will Schroeder <will.schroeder at kitware.com<mailto:will.schroeder at kitware.com>>
Date: Saturday 23 April 2016 at 13:56
To: Favre Jean <jfavre at cscs.ch<mailto:jfavre at cscs.ch>>
Cc: VTK Developers <vtk-developers at vtk.org<mailto:vtk-developers at vtk.org>>
Subject: Re: [vtk-developers] SPH Interpolator

Hey Jean-

Here's a longer email than you'd probably like; hopefully it will give you a better idea of where this is heading. I want to emphasize that this a work in progress, and feedback from those more knowledgeable in this area is very welcome.

The recent efforts are partly due to the generous support of a collastomer (collaborator/customer), and the personal interest of myself and others at Kitware. For example, besides the SPH additions, I've implemented a bunch of point cloud algorithms in the vtkPointCloud remote module, and Ken Martin has put together a kick-butt point cloud renderer (also a VTK remote module vtkLODPointCloudMapper). We've also got some efforts going on with VeloView<http://www.paraview.org/Wiki/VeloView>, which is a derivative ParaView application for point cloud visualization and analysis. Finally, I've been reworking many algorithms using threaded approaches (vtkSMPTools) including the new vtkStaticPointLocator which is typically much faster than any other point locators and is critical to the speed of the SPH operations. (Caveat: I'm guessing that highly dynamic data where point densities vary by orders of magnitude e.g., astrophysics will require more work, probably a hierarchical static point locator much like an AMR grid.)

Okay with that preamble behind me, the current thinking is to support three workflows.

1. Fixed kernel size / specified smoothing length h (also called spatial step). The cutoff distance (sphere around an interpolated point) is a function of the SPH kernel. For example, a quintic kernel has cutoff distance 3*h. The volume is approximated as h^3. This has been tested with real simulation data incompressible CFD and is being validated/verified against simulation data so I have confidence that we are in good shape here. So far we have implemented cubic, quartic and quintic SPH kernels.

2. Volume defined on a pointwise basis using a pair of mass/density arrays. The API is there to support this, but I believe the last I left it there is a touch more work to do to get it working correctly, and of course testing is necessary. We could really use some data to test this properly. It's probably a few hours work and I plan to get back to finishing it off once I return from some travel next week.

3. Currently in #1 the kernel neighbors are found in the cutoff sphere, with the sphere radius constant across each point. However it would be possible to support a variable radius per point through a provided external array. This is not hard to do, but I really need data before proceeding.

My understanding from reading references like Price<http://arxiv.org/pdf/1012.1885.pdf> is that another, less desirable possibility for basis is to use N neighbors, rather than the neighbors contained in the cutoff sphere. I have done this in vtkPointInterpolator / vtkGeneralizedKernel classes, which are close cousin classes to vtkSPHInterpolator / vtkSPHKernel. So this could easily be extended to the SPH side as well.

The current implementation uses a gather method, and is done in parallel using vtkSMPTools. For each point to be interpolated, the basis neighbors around the point are retrieved using a provided abstract point locator (typically an instance of vtkStaticPointLocator but others are possible). The provided kernel is then invoked to perform the interpolation. So you can see that a vtkSMPTools::For(0,npts,functor) works really well to thread the operation (typically TBB is used because I like the load balancing that it provides).

So any suggestions you might have, or data you can provide would be more than welcome. Also note that more testing is desired, for example the cubic and quartic SPH kernels need more testing. Bill Lorensen has been helping a lot in that he is adding some unit tests for the generalized kernels and has already uncovered a couple of bugs (which have been fixed).

Anyway I leave for a week of travel early this coming week. If you want to continue the conversation when I return in early May that would be great.

Best,
W

On Fri, Apr 22, 2016 at 10:28 AM, Favre Jean <jfavre at cscs.ch<mailto:jfavre at cscs.ch>> wrote:
Hello

First of all, congrats on this new functionality. I am using it with results from an SPH simulation using the GADGET code[1]

Our scientists are asking some pretty specific questions about the SPH kernels:

Do you adopt a fixed kernel size for all particles or do you calculate it by adopting a certain number of neighbours? Do you use the scatter method (distributing each SPH particle on the grid) or the gather method (gathering the contribution of all SPH particles which are close enough to each grid point)?

They further comment that "to calculate physical quantities, such as column density etc. , these details become important. For instance, it is important to know if the SPH interpolation scheme which is used is conserving the total mass or not when we map the mass from SPH distribution onto a uniform grid, which depend on the details of the SPH interpolation."

If the Density and Mass arrays are set, [I use]:

interpolator = vtkSPHInterpolator()
interpolator.SetKernel(sphKernel)
interpolator.SetDensityArrayName("Density")
interpolator.SetMassArrayName("Mass")

can I assume the above to be true?

Thanks for clarifying these points. I anticipate some really good use of this new Interpolator.

[1] http://wwwmpa.mpa-garching.mpg.de/gadget/

-----------------
Jean/CSCS

_______________________________________________
Powered by www.kitware.com<http://www.kitware.com>

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Search the list archives at: http://markmail.org/search/?q=vtk-developers

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/vtk-developers

--
William J. Schroeder, PhD
Kitware, Inc. - Building the World's Technical Computing Software
28 Corporate Drive
Clifton Park, NY 12065
will.schroeder at kitware.com<mailto:will.schroeder at kitware.com>
http://www.kitware.com
(518) 881-4902
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtk-developers/attachments/20160423/7378d1e1/attachment.html>