ENH: Allow for independent particle attribute values on initialisation and advection

[Patol75]

Currently, the particle implementation in Fluidity does not allow the user to set initial attribute values for particles initialised during the simulation. Instead, attributes are set to zero and immediately updated according to the provided constant, python or python_fields field in diamond. Such an approach does not allow to distinguish particles that have just been initialised during the current timestep from particles that were already present. In particular, this is an issue if particles carry some attributes that represent fields such as concentrations, which can require a non-zero initial value. Importantly, this behaviour only affects particles initialised during the simulation and not particles spawned as a result of control volume thresholds.

In this PR, a strategy to address the above problem is implemented. Within diamond, each particle attribute now contains a value_on_spawn and a value_on_advection fields. value_on_spawn is optional; if not provided, value_on_advection will be executed at particle initialisation. Before the timestep loop starts and whenever new particles are initialised during the simulation, value_on_spawn is executed. At any other time, value_on_advection is executed.

The coexistence of value_on_spawn and value_on_advection brings a potential conflict regarding the current implementation for storing old particle attribute values if python_fields are used. Here, the choice is made to only consider store_old_attribute under value_on_advection; if both value_on_spawn and value_on_advection are provided, turning on store_old_attribute under value_on_advection automatically stores old attribute values whether value_on_spawn or value_on_advection is executed. store_old_attribute is, nevertheless, available under value_on_spawn, but it is not being used, as mentioned in the updated schemas.

Within the timestep loop, initialise_particles_during_simulation is now executed after particle_cv_check and update_particle_attributes_and_fields. particle_cv_check deals with particle spawning, and it seems to me that simulations making use of spawning should not require the initialise during simulation option. update_particle_attributes_and_fields comes before initialise_particles_during_simulation as it updates attributes on particles that were present at the previous timestep. In turn, initialise_particles_during_simulation now both initialises new particles in the domain according to the initialise_during_simulation field in diamond and sets their attribute values according to the value_on_spawn field in diamond.

Such proposed changes are backwards-incompatible for two reasons:

• schema structure is different
• specific old attribute values logic for particles initialised during a given timestep is not available anymore (more details after the CI runs)

Additionally, considering the terminology in Fluidity, I am open to renaming value_on_spawn to value_on_initialisation, which seems more appropriate on second thought.

Finally, many thanks to @angus-g for some preliminary work on this PR.

---

Looking at the CI, tests suites have returned for Bionic and Focal. Groovy is not supported upstream anymore, which echoes to #338. Medium tests passed, so no worries there. For short tests, particle_attribute_array fails on Bionic, but not on Focal, because the default NumPy on Bionic does not have the newer random generator. Considering that Bionic support should drop pretty soon, this does not sound like an issue to me. lagrange_tracer_lowcourant failed because I missed updating its flml. particle_initialise_after_start fails, and the failures relate to my second backwards-incompatible reason above. Using the previous implementation, attributes of particles initialised during the simulation would be updated during the same timestep, whereas now there are not anymore. As a result, particles initialised at timestep 2, for example, have their attribute correctly equal to 0 and not 1, as it was before. Finally, diamond_validation fails because I missed updating some of the flml that include particles, as mentioned above.

---

Following the second CI, test_get_rotation_matrix reported one failure in Unit Focal; I have not investigated. particle_rayleigh_taylor_mu10, which I updated in the second commit, yields Sqr_entrainment_error_24 = 0.0020970997668790683 and Sqr_entrainment_error_48 = 0.0015043119771471875, and therefore falls short of the 1e-3 threshold value. When I run the same test on the supercomputer in Canberra, I get Sqr_entrainment_error_24 = 0.00014354568992100723 and Sqr_entrainment_error_48 = 0.022645421135390705. Something I am not sure about is if the mesh is supposed to be structured or not in these simulations, and I guess that would change the results. Remaining tests as before or solved.

[Patol75]

Thanks @angus-g for your time.

Yes, the logic is the one you implemented. Nevertheless, I tried to simplify a few things here and there where it made sense. The main thing that was missing was the per-processor counter for the number of particles initialised on each processor during a timestep, which I added. I also remember an issue with global not being initialised which led to some troubles that, for some reason, do not show up in the current main branch. The following code illustrates the problem:

program example
  implicit none

  logical :: global

  print *, global, global .eqv. .true., global .eqv. .false.
  call display_behaviour(global)
  print *, global, global .eqv. .true., global .eqv. .false.

contains

  subroutine display_behaviour(global)
    logical, intent(inout), optional :: global

    if (present(global)) then
       global = .not. global
    end if

  end subroutine display_behaviour

end program example

If run as is, it returns:

T T T
T T T

Instead, if we have logical :: global = .true. on line 4, then we get:

T T F
F F T

I think I have corrected that problem by setting global=.true. instead of global in what is now line 520 of femtools/Particles.F90, but I would not mind a second opinion.

No worries about the tests, and I agree that we should have a more careful look at particle_rayleigh_taylor_mu10.

Also, thanks for pointing out the trick about removing whitespaces! Atom always suppresses them automatically, which is helpful, but without that trick, diffs often look horrible.

As a side note, I noticed I forgot to remove ewrite(2, *) "Here_12" from line 331 of femtools/Node_Owner_Finder_Fortran.F90. I will correct in a later commit.

---

I forgot to mention it yesterday, but I have added a few ! FIX ME - Does not print in femtools/Detector_Tools.F90. I believe these ewrite calls, which should report Python failures, do not work properly. If I replace ewrite(-1, *) by ewrite(2, *), I see the messages in the log files, but the actual error does not show up.

[stephankramer]

Such an approach does not allow to distinguish particles that have just been initialised during the current timestep from particles that were already present. In particular, this is an issue if particles carry some attributes that represent fields such as concentrations, which can require a non-zero initial value.

Can you expand a little on your motivation. Why do we need to distinguish? Why does our inability to do so cause attributes to be initialised with a zero value? I can have a guess, but I think there are multiple assumptions here about your use case that are worth spelling out - a typical use case example for instance would help.

Am I right in assuming that for a lot of cases of particle attributes they are set using constant, or set from in python from a combination of other current attributes and interpolated field values - and it seems for these cases this issue doesn't apply? So I'm guessing the use case is when you're dealing with old and new values (so we're in effect trying to do some poor man's ODE time integration)?

If this is the case then I think I would much prefer a backwards-compatible solution with an optional initial_value level, where the existing top level is used during the simulation, rather than having a separate level for both the initial value and the value during the rest of the simulation. If I think of the use case of a user that merely wants to use particles to passively track some properties, it might be a little confusing they have to set something under "value_during_advection" ("what does that mean? a particle is always being advected?")

As a separate point, coming back to the whitespace issue: yes you can work around the whitespace in github diff issues - but that's just circumventing the problem: you should not be comitting arbitrary whitespace changes (even if it improves the indentation). The whitespace display in the diff is just one example of how it breaks development tools of other people that have to deal with your code, it also messes up the annotation of who changed which line in history which is an important tool in large projects like this, and even more importantly it increases the change of merge conflicts.

As a more general point you should try to be extremely conservative when changing code, and only change the code that is directly relevant to the thing you are trying to implement/fix. Otherwise you just create a lot of work for others and yourself in the future. Any code you change comes with the risk of introducing bugs that need to be dealt with. So "cleaning up code" and rewriting a bits of code purely because you see it has been done in a bit of a sub-optimal way, or there's some cleverer way to do it, is not really helpful.

Could you go through your PR with that in mind and have critical look at all the changes line by line (you should always do that before you request a PR), and remove anything that's not relevant - and if we're talking about reducing potential conflicts it might be worth squash-rebasing it, so the spurious changes disappear from history.

[Patol75]

Can you expand a little on your motivation. Why do we need to distinguish? Why does our inability to do so cause attributes to be initialised with a zero value? I can have a guess, but I think there are multiple assumptions here about your use case that are worth spelling out - a typical use case example for instance would help.

Am I right in assuming that for a lot of cases of particle attributes they are set using constant, or set from in python from a combination of other current attributes and interpolated field values - and it seems for these cases this issue doesn't apply? So I'm guessing the use case is when you're dealing with old and new values (so we're in effect trying to do some poor man's ODE time integration)?

The changes proposed here are indeed motivated by the integration of initial-value problems (IVP). In this regard, being able to distinguish between newly initialised particles and already present particles is crucial. In the current main branch, particles initialised during a timestep (this also applies before the timestep loop starts) have all their attributes hard-coded to 0. Immediately after, all particles (including the ones just initialised) have their attributes updated based on the diamond entry provided by the user. In writing this diamond entry, the user cannot distinguish between newly initialised particles and already present particles. As a result, it is impossible to concurrently set the initial state of the IVP on the newly initialised particles and have the integration of the IVP on the already present particles progress one step further. In other words, the issue this PR deals with occurs whenever a particle attribute needs to be initialised to a value different from the one it would be getting through the update routine.

Regarding how particle attributes are set, this very much depends on the problem that you are dealing with in your simulation. If particles are used to track materials, then a simple constant entry suffices. But if particles are used to solve equations, then a python_fields is likely required. I do not believe it is right to say that the former scenario is the one that occurs the most: I have been using particles for almost 4 years, and I have only dealt with the latter scenario. @Cmath2, who developed the particles, mainly used the former scenario when he was benchmarking the particle scheme, but later only used the latter scenario when he was tracking seismic anisotropy in his models. As a note, I do not believe it is currently possible to set particle attributes based on the current value of other attributes (at least that is what I understood from multiple conversations with @Cmath2).

I paste below a typical scalar_attribute_array that I use for my simulations. The first snippet corresponds to value_on_spawn and the second one to value_on_advection. These snippets make use of current finite-element field values and old attribute values (i.e. from the previous timestep). Before the changes proposed in this PR, both snippets would be combined such that the code executed in value_on_spawn would be subjected to an if t == 0:. This way, particles initialised before the timestep loop starts can be distinguished from already present particles, but indeed there are not any already present particles at the beginning of the simulation anyway. Later in the simulation, there is no equivalent to if t == 0: to distinguish newly initialised particles from already present particles.

def val(X, t, dt, fields, n):
    from numba import float64
    from numba.typed import Dict
    from numba.types import unicode_type
    from numpy import array, clip, ones, zeros, zeros_like
    from scipy.constants import g

    from Melt import Katz
    from MeltChemistryFunctions import non_zero_initial_melt, run_integrator

    from constants import (adiab_grad, cs_0, domain_dim, mant_temp,
                           melt_inputs, rho_continent, rho_crust, rho_mantle)

    depth = clip(domain_dim[2] - X[2], 0, domain_dim[2])
    lab_depth = fields['DepthLAB']
    crust_depth = min(4.1e4, lab_depth / 2)

    if depth < crust_depth:
        presGPa = rho_crust * g * depth / 1e9
    elif depth < lab_depth:
        presGPa = (rho_crust * g * crust_depth
                   + rho_continent * g * (depth - crust_depth)) / 1e9
    else:
        presGPa = (rho_crust * g * crust_depth
                   + rho_continent * g * (lab_depth - crust_depth)
                   + rho_mantle * g * (depth - lab_depth)) / 1e9

    temp = fields['Temperature'] + adiab_grad * depth

    F = Katz().KatzPT(presGPa, temp, inputConst=melt_inputs)
    if F == 0:
        return (F, F, 0, presGPa, temp, 0,
                *ones(cs_0.size, dtype=bool), 1e-6, F, presGPa, False,
                *zeros(6), *zeros(6), *cs_0, *zeros_like(cs_0))

    part_arr = Dict.empty(key_type=unicode_type, value_type=float64[:])
    part_arr["melt_fraction"] = array([F])
    part_arr["pressure"] = array([presGPa])
    part_arr["temperature"] = array([temp])

    lab_pres = (rho_crust * g * crust_depth
                + rho_continent * g * (lab_depth - crust_depth)) / 1e9
    dTdP_GPa = adiab_grad / rho_mantle / g * 1e9

    part_arr = non_zero_initial_melt(
        part_arr, fields['Temperature'], mant_temp, lab_pres, dTdP_GPa)

    mask_ode, ode_dX, Fn, pn, Dn, D_bulk, P_bulk, cs, cl = run_integrator(
        part_arr, ones(cs_0.size, dtype=bool), 1e-6,
        part_arr["melt_fraction"][0], part_arr["pressure"][0], zeros(6), cs_0)

    return (F, F, 0, presGPa, temp, 0,
            *mask_ode, ode_dX, F, presGPa, False, *Fn, *pn, *cs, *cl)

def val(X, t, dt, fields, n):
    from numba import float64
    from numba.typed import Dict
    from numba.types import unicode_type
    from numpy import array, clip, logical_xor
    from scipy.constants import g

    from Melt import Katz
    from MeltChemistryFunctions import run_integrator

    from constants import (adiab_grad, cs_0, domain_dim, melt_inputs,
                           rho_continent, rho_crust, rho_mantle)

    depth = clip(domain_dim[2] - X[2], 0, domain_dim[2])
    lab_depth = fields['DepthLAB']
    crust_depth = min(4.1e4, lab_depth / 2)

    if depth < crust_depth:
        presGPa = rho_crust * g * depth / 1e9
    elif depth < lab_depth:
        presGPa = (rho_crust * g * crust_depth
                   + rho_continent * g * (depth - crust_depth)) / 1e9
    else:
        presGPa = (rho_crust * g * crust_depth
                   + rho_continent * g * (lab_depth - crust_depth)
                   + rho_mantle * g * (depth - lab_depth)) / 1e9

    temp = fields['Temperature'] + adiab_grad * depth

    old_attrs = fields['old%katz_mckenzie_bdd21_']

    if abs(presGPa - old_attrs[3]) < 1e-9:
        return (old_attrs[0], old_attrs[1], 0, presGPa, temp, 0, *old_attrs[6:])

    temp_grad = (temp - old_attrs[4]) / (presGPa - old_attrs[3])
    sol = Katz().KatzPTF(old_attrs[3], presGPa, old_attrs[4], old_attrs[0],
                         temp_grad, inputConst=melt_inputs)
    T, F = sol(presGPa)
    F = 0 if F < 0 else F
    if F > old_attrs[1]:
        melt_rate = (F - old_attrs[1]) / dt * 8.64e4 * 365.25 * 1e6
        temp_src = (T - temp) / dt
    else:
        return (F, old_attrs[1], 0, presGPa, temp, 0, *old_attrs[6:])

    part_arr = Dict.empty(key_type=unicode_type, value_type=float64[:])
    part_arr["melt_fraction"] = array([old_attrs[1], F])
    part_arr["pressure"] = array([old_attrs[3], presGPa])
    part_arr["temperature"] = array([old_attrs[4], T])

    mask_ode = array(old_attrs[6:6 + cs_0.size], dtype=bool)
    if not mask_ode.any():
        return (F, F, melt_rate, presGPa, T, temp_src, *old_attrs[6:])

    ode_dX = old_attrs[6 + cs_0.size]

    Fn_old = array(old_attrs[10 + cs_0.size:16 + cs_0.size])
    if F / melt_rate < 4 and (old_attrs[2] > 0 or old_attrs[1] == 0):
        X_0, P_0 = old_attrs[1], old_attrs[3]
    elif old_attrs[9 + cs_0.size]:
        X_0, P_0 = old_attrs[1], old_attrs[3]
    else:
        X_0, P_0 = old_attrs[7 + cs_0.size:9 + cs_0.size]

    pn_old = array(old_attrs[16 + cs_0.size:22 + cs_0.size])
    cs_old = array(old_attrs[22 + cs_0.size:22 + 2 * cs_0.size])

    mask_ode, ode_dX, Fn, pn, Dn, D_bulk, P_bulk, cs, cl = run_integrator(
        part_arr, mask_ode, ode_dX, X_0, P_0, pn_old, cs_old)

    return (F, F, melt_rate, presGPa, T, temp_src,
            *mask_ode, ode_dX, X_0, P_0, logical_xor(Fn_old, Fn).any(),
            *Fn, *pn, *cs, *cl)

If this is the case then I think I would much prefer a backwards-compatible solution with an optional initial_value level, where the existing top level is used during the simulation, rather than having a separate level for both the initial value and the value during the rest of the simulation. If I think of the use case of a user that merely wants to use particles to passively track some properties, it might be a little confusing they have to set something under "value_during_advection" ("what does that mean? a particle is always being advected?")

What you are proposing here is what @angus-g had initially coded. However, I found that confusing. Why would two diamond entries that essentially achieve the same thing (set particle attributes values) be treated differently in the schema? I believe that they should not. As a result, I designed value_on_spawn and value_on_advection. As mentioned earlier, value_on_spawn would be better named value_on_initialisation. In any case, for simulations where spawning is disabled, particles can essentially do three things: be initialised, be advected, be deleted. A particle can be initialised before the timestep loop starts or during the simulation, as requested by the user. In both cases, the particle is only initialised, it is not advected. A particle is only advected if it was initialised at a previous timestep. During advection, if a particle moves outside of the domain, it is deleted. Therefore, I believe that value_on_initialisation and value_on_advection are meaningful names that reflect what is happening in the simulation. However, I agree that the current descriptions available in Diamond are not satisfactory, and I am happy to update them to make things clearer. Furthermore, value_on_initialisation is an optional field, and this is where the updated Diamond description would help by saying that (i) value_on_initialisation is only required if a particle attribute needs to be initialised to a value different from the one set in value_on_advection and (ii) value_on_advection is by default also executed before the timestep loop starts to determine the initial value of attributes. Additionally, I would like to emphasise that the behaviour displayed by the particle_initialise_after_start test in the current main branch is definitely undefined to the mere user. Indeed, the mere user has no idea that the attributes of the newly initialised particles are hard-coded to 0. Therefore, what should the mere user expect about this attribute that is incremented by 1 at each timestep? Making value_on_initialisation available, and clearly differentiating it from value_on_advection, solves this problem.

As a separate point, coming back to the whitespace issue: yes you can work around the whitespace in github diff issues - but that's just circumventing the problem: you should not be comitting arbitrary whitespace changes (even if it improves the indentation). The whitespace display in the diff is just one example of how it breaks development tools of other people that have to deal with your code, it also messes up the annotation of who changed which line in history which is an important tool in large projects like this, and even more importantly it increases the change of merge conflicts.

I agree that git blame is an important tool. However, it is still possible in git blame to see changes before the most recent change. Additionally, as is mentioned in this conversation, (i) GitHub provides a workaround for whitespace-only changes on its website and (ii) GitHub's own editor, Atom, automatically makes these whitespace-only changes by default. This, to me, means that GitHub encourages people to remove spurious whitespaces from the code they develop and to clean existing whitespaces from already-hosted codes. I am not aware of why such changes would break the development tools of other people.

As a more general point you should try to be extremely conservative when changing code, and only change the code that is directly relevant to the thing you are trying to implement/fix. Otherwise you just create a lot of work for others and yourself in the future. Any code you change comes with the risk of introducing bugs that need to be dealt with. So "cleaning up code" and rewriting a bits of code purely because you see it has been done in a bit of a sub-optimal way, or there's some cleverer way to do it, is not really helpful.

I do not believe that this is my approach. My approach, while implementing the main purpose of my PR, is to fix whatever seems to be broken or about to be broken (i.e. deprecated) and that also impacts my changes. I am not changing code for the fun of doing it, but rather to improve things where it is possible. It seems more logical to me to do it at the same time as the development of the main PR purpose as it avoids later work when one or another identified feature eventually breaks. Nonetheless, it is true that, in the case of Fluidity which has lots of spurious whitespaces, this leads to committing files that only contain whitespace changes. This can happen, for example, when tracking segmentation faults. I usually employ print statements, and yes, sometimes I put these print statements in source files that end up not requiring any changes, but the simple fact of saving the modified file in Atom automatically gets rid of all the whitespaces. Additionally, when I change indentation in source files, it is not just because it makes the code more beautiful-looking, but instead mainly because it makes the code readable. Sometimes, when I open these older Fortran or C or C++ source files, (I have to be honest) it is simply unreadable, and for example you can lose 5 minutes trying to determine if a block of code is included in an if condition or not.

Could you go through your PR with that in mind and have critical look at all the changes line by line (you should always do that before you request a PR), and remove anything that's not relevant - and if we're talking about reducing potential conflicts it might be worth squash-rebasing it, so the spurious changes disappear from history.

I do not believe there are any spurious changes in the present PR (although you could argue this for some other PR I have currently opened in this repository). I have already pointed out in my previous message that I have forgotten to remove a print statement, and I will remove it in a further commit.

[Patol75]

I have pushed a fix for an issue regarding adapt_at_first_timestep. This routine can be executed as part of a simulation restarted from a checkpoint. When particles are loaded from the checkpointed h5part file, value_on_advection should be used.

[Patol75]

I have gone through the changes again, and I have realised that I had forgotten to update some parts of the particle code, hence why particle_rayleigh_taylor_mu10 was failing. I have pushed a commit that should address what I had forgotten to include. CI now returns all green on Focal and Impish, but it looks like I have put the NumPy upgrade in the wrong Dockerfile for Bionic.

Something that I do not understand is why constant attributes are getting what looks like special treatment in multiple locations. For instance, why is there a routine called initialise_constant_particle_diagnostics in assemble/Particle_Diagnostics.F90, but there is no equivalent one for python/python_fields attributes? Also, in update_particle_subgroup_options within femtools/Particles.F90, I have commented out ~25 lines of code associated with checkpointing of constant attributes, and it seems it did not cause any test to break. @angus-g Do you have some insight about how constant attributes are treated here?

[Patol75]

Final update on my side, or at least it should be:

• Apart from the Impish failures due to the latter's end of life, CI now passes successfully.
• I have renamed value_on_initialisation and value_on_advection to initial_attribute_value and attribute_value, while also clarifying their descriptions in the schema.
• I have fixed the behaviour when restarting from a checkpoint: whether adapt_at_first_timestep is called or not, attribute_value should always be used. Consequently, initial_attribute_value should only be called before the timestep loop starts when particles are initialised through Python and whenever initialise_during_simulation is called. The implementation for the latter case implicitly assumes that any simulation using particle spawning should not use the initialise_during_simulation feature, which sounds reasonable to me.
• I have removed entirely checkpoint_from_file from all particle-related code. It was simply not used at all. Also, I now better understand why constant attributes get special treatment: the particle code implicitly assumes that all material volume fraction fields calculated from particles should depend on constant attributes. As these attributes need to be calculated before MVF fields are initialised, constant attributes are dealt with separately in a few instances. This is very likely fine. However, there may be an issue with adapt_at_first_timestep. It seems to me that MVF fields derived from particles are not calculated within adapt_at_first_timestep and, as a result, I do not see how such fields could contribute to the adaptivity iterations. I was expecting to see calculate_particle_material_fields being called as is the case in Fluids.F90, but for some reason, it is not.
• I have updated a couple of XML files to reflect some deprecations within NumPy and SciPy that generated errors when I updated NumPy and later on SciPy within the Bionic Actions Dockerfile.

Happy to hear your thoughts @stephankramer and @angus-g.

[Patol75]

Hi @stephankramer and @angus-g,

When the two of you have some time, would you mind giving another look at the present PR? We have a few group members here at ANU who have been using this branch for quite a while now without issues; it would be great if it could make its way into main, given that it extends particle capabilities.

I merged main into the branch last month, which triggered the CI. It returned all green, which I am sure you will agree is great. Please let me know if there are outstanding issues.