Particles "jumping"

[jhill1]

Hello, I've been running some tidal sims (SWE solver) in fluidity with particles turned on. These particles are simple ones, they don't interact with the fields, just drift along.

Analysing the paths, they seem to "jump" occasionally. An example is attached.

The particles timestep is 60 seconds and the labels in the figure are the time in seconds. You can see the particles at 97920 (right) and then 97980 (left) and the gap between.

The FLML has the following for particles:

  <particles>
    <particle_group name="Humans">
      <velocity_field name="Velocity"/>
        <particle_subgroup name="Site_3">
        <number_of_particles>
          <integer_value rank="0">100</integer_value>
        </number_of_particles>
        <initial_position>
          <python>
            <string_value type="code" language="python" lines="20">def val(t):
  det = []
  f = open("detectors_3.csv","r")
  for line in f:
    if line.startswith("X"):
      # skip header
      continue
    data = line.split(",")
    x = float(data[0])
    y = float(data[1])
    det.append([x,y])
  return det</string_value>
          </python>
        </initial_position>
      </particle_subgroup>
    </particle_group>
    <write_nan_outside_domain/>
    <lagrangian_timestepping>
      <subcycles>
        <integer_value rank="0">2</integer_value>
      </subcycles>
      <search_tolerance>
        <real_value rank="0">1e-7</real_value>
      </search_tolerance>
      <rk4_guided_search/>
    </lagrangian_timestepping>
  </particles>

and the h5part files (to generate the csv file to create the image above, for a single particle path) is:

import h5py

with h5py.File("Palaeo_redsea_dets.particles.Site_1.h5part", "r") as f:
        n_times = len(list(f.keys()))
        timestep = 60.
        for t in range(0,n_times):
            # Get the data; this is all 100 coords (x) at current time
            name = 1
            ts = 'Step#'+str(t)
            print(t*timestep, f[ts]['x'][1], f[ts]['y'][1])
            name += 1

The particles do not leave the domain (no nans, and look ok with respect to the domain boundaries). The above FLML snippet has 1 group, but there are three in total. Is this a bug, have I erred in the extraction script or have I got the wrong options?

thanks, Jon

[jhill1]

Got the same with an Euler directed search. Purple is Euler, blue RK. Jump goes from 12240 (bottom) to 12300 (top) on both.

[Patol75]

Hi @jhill1,

Thanks for the report, that is indeed an unexpected behaviour. That being said, it seems to me that you are not using the latest implementation of particles in Fluidity. For example, you do not seem to have a particle_io attribute in your flml, and you have a velocity_field attribute, which I am pretty sure has been removed. Regarding the lagrangian_timestepping, I have been using rk4_guided_search with subcycles and search_tolernace set to 3 and 1e-12, respectively.

Hope that helps.

[jhill1]

Thanks @Patol75 Looks like I had different versions on different computers. I had the latest main checkout on the HPC system where the above was run, but made the FLML on my desktop, which had an older version. Updated the flml and kicked it off. Will see what happens!

[jhill1]

Unfortunately same effect (at exactly the same time).

Will see if I can find time to debug this, but any help appreciated.

[Patol75]

Does your simulation include adaptivity? If yes, try turning it off and see if the behaviour persists. Similarly, is your simulation serial or parallel? If parallel, a serial run would also help.

Also, is there any fundamental difference between the type of simulation you are running and Stokes flow within Fluidity? Particles were designed using Stokes flow, so we could be missing something if the answer to the previous question is yes.

Regarding your h5part, if you load it in a recent version of ParaView, can you confirm that the proc_id and id of the particle shown above are the same at 12240 and 12300? You can use the Hover Points On functionality to interactively check.

[jhill1]

It's a non-adaptive, parallel (32), shallow water equation simulation. When visualising in paraview 5.9 (didn't know you could do that btw - very nice!); it all looks great. No jumping. Clearly it's my script to convert to shapefiles that's he issue here. So I think there is no bug and the error lies between keyboard and chair :-)

Thanks for the help!

J

[jhill1]

For the record, this now works to pull out particle 96:

import h5py

with h5py.File("../simulations/palaeo/simulation_888/Palaeo_redsea_dets.particles.Site_1.h5part", "r") as f:
        n_times = len(list(f.keys()))
        timestep = 60.
        for t in range(0,n_times):
            # Get the data; this is all 100 coords (x) at current time
            ts = 'Step#'+str(t)
            idx = list(f[ts]['id']).index(96)
            print(t*timestep, f[ts]['x'][idx], f[ts]['y'][idx])

Clearly the indices of the particles moves around, which then explains the "jumping" I saw. Doing the above keeps the right order.

[Patol75]

Yes, a particle is uniquely identified by its id and its proc_id, so this is what you should rely on when you post-process the h5part file.