Lagrangian module enhancements

[paspro]

1. Implementation of the full turbulence dispersion scheme with rotation of the frame of reference.

2. Modification of CS flags to reflect these changes.

• modcpl = 1 now enables the full turbulence dispersion scheme and this is its default value. It takes a value of 0 when the user specifies that the particles are supposed to be fluid-like tracker particles in which case the scheme is disabled. The value of modcpl also affects the values of two other flags:

• modcpl = 1 => idistu = 1 (turbulent dispersion is enabled), idiffl = 1 (turbulent diffusion is disabled)

• modcpl = 0 => idistu = 1 (turbulent dispersion is disabled), idiffl = 1 (turbulent diffusion is enabled)

• irdla has been removed.

3. The Code_Saturne GUI has been modified to reflect the changes above.

4. Modification of the wall boundary conditions for particles.

[YvanFournier]

Hello,

As we use mostly a "rebase" logic, grouping commits in slightly larger blocs, reviewing the fine-grained commits is not easy, so I may have missed some things.

Comparing the reference master branch and the Renuda head, I have doubts about the update in cs_lagr_sde.c, even though the Renu branch was updated. Recent commits in that branch provided by C. Henry mention fixes in rotation formula, while this branch only mention updates, so I guess this is a merge issue, and kept the previous version, assuming that this part is not impacted by the Renuda enhancements beyond updates.

Also, the changes on modcpl conflict with fixes in recent days on this. The previous logic was incorrect, but I believe was fixed: since modcpl could only be active when modcpl >= idistnt, I fixed the tests (commit 2709c229) so that whatever the values, the model would only be active when the time step is above both constraints. But limiting this to "1" values is incorrect, as idistnt can be > 1 (and to my knowledge, this must be kept so as to allow updating statistics only once a number of initialization iterations have been run). Allowing modcpl > idistnt allows activating this option later than idisnt, though it might not be so useful (I need to check this with Martin and Mathieu). In any case, allowing this would be an avanced option, so setting modcpl to 1 (and letting the code modify that automatically to modcpl = idisnt) is good enough for GUI-defined setups.

I did not find the automatic changes of idistu and idiffl on the code, so also need to check this. If the requirement was for this to be fully automatic, then these keywords could be removed and the test on modcpl activation used instead, unless we want to keep them separate for advanced settings.

I rebased those changes (committing some additional style cleanup separately and pushing them first), and attach the resulting patch relative to the code_saturne master branch, sha1 4103d5e0.

This can be used for further review.

0003-Lagrangian-module-enhancements-based-on-Renuda-devel.patch.txt

[YvanFournier]

Also, since the test cases were placed in a separate repository (for easier merging in the future), should I obtain them from the history of this branch or from another repository ?

[paspro]

Also, since the test cases were placed in a separate repository (for easier merging in the future), should I obtain them from the history of this branch or from another repository ?

You can obtain them from the history of this branch.

[paspro]

Also, the changes on modcpl conflict with fixes in recent days on this. The previous logic was incorrect, but I believe was fixed: since modcpl could only be active when modcpl >= idistnt, I fixed the tests (commit 2709c22) so that whatever the values, the model would only be active when the time step is above both constraints. But limiting this to "1" values is incorrect, as idistnt can be > 1 (and to my knowledge, this must be kept so as to allow updating statistics only once a number of initialization iterations have been run). Allowing modcpl > idistnt allows activating this option later than idisnt, though it might not be so useful (I need to check this with Martin and Mathieu). In any case, allowing this would be an avanced option, so setting modcpl to 1 (and letting the code modify that automatically to modcpl = idisnt) is good enough for GUI-defined setups.

The project requirements regarding this issue are the following:

1. Make the complete turbulent dispersion model the default and only available choice (modcpl = 1)
2. There is no more a non-complete turbulent dispersion model available.
3. Statistics are required for (1) to work and therefore they should be enabled by default.
4. The user does not have the choice to enable the complete turbulent dispersion model after a specific number of time steps. It is enabled from step 1, always.
5. This model is disabled only if the user wants the particles to be considered as fluid-like tracer particles (modcpl = 0).

You may want to change these flags but the idea is that the user has the choice for:

1. Regular particles -> we set modcpl = 1, idistu = 1, idiffl = 0 (you can call these settings turbulent_dispersion = true)
2. Fluid-like tracer particles -> we set modcpl = 0, idistu = 0, idiffl = 1 (you can call these settings turbulent_dispersion = false)

I hope this clarifies things.

[YvanFournier]

Hello,

Yes, in this case I would probably go for a turbulent dispersion=true/false option rather than 3 separate options.

I am still waiting for confirmation from Mathieu that nothing was forgotten in the requirements: so far, (and in some test cases), when using the Lagrangian "steady flow" model, statistics are not reset after each iteration, but averaged over many iterations. We usually wait for the flow to be somewhat stabilized and particles to populate the whole domain before activating statistics. Otherwise, the statistics will be biased by the initialization steps. Using the "unsteady model", statistics are rest at each time step, so the issue disappears.

We might also be able to solve this removing idstnt and using only nstits, which determines the threshold to when we switch from unsteady to steady statistics. Which has the advantage of removing one more key word.

Also, if we now consider that particles are either "solid" or "tracer", adding an attribute to particles would allow mixing both options for different particle sets. I'll suggest this on our side, or check if it would be useful.

[YvanFournier]

I have further checked and edited (to ease diffs with the master branch) the patch, and now arrive at the attached patch.

I am not sure the rotation matrix coefficients are correct: in the master branch, commit 1795093a fixed some orientations, and the signs in the Renuda branch matched the old ones. I merged so as to have the fixed one in the stochastic and Jeffery models, and keep the Renuda ones in the sphere model with modcpl=1, but this seems strange.

Also, the sphere model with modcpl = 1 and the stochastic model seem very similar, except for the way the orient_loc vector is computed (the axis is named main_axis in the Renuda model, and renamed singularity_axis from the Frenglish axe_singularity name in the POPART project's spheroid model, but that is probably just a choice of name).

Given this, I wonder whether the modcpl and shape model options should be combined into a single keyword (at least based on the code in cs_lagr_sde.c, and assuming we sort out the rotation matrix orientation issues). I'll try to contact Jean-Pierre, as I believe he followed both developments.

I also need to automate the settings for idistu and idiffl, or base matching tests on modcpl, as these values do not seem to be modified automatically in the Renuda branch (or I can't find where).

0001-Lagrangian-module-enhancements-based-on-Renuda-devel.patch.txt

[paspro]

I have further checked and edited (to ease diffs with the master branch) the patch, and now arrive at the attached patch.

I am not sure the rotation matrix coefficients are correct: in the master branch, commit 1795093 fixed some orientations, and the signs in the Renuda branch matched the old ones. I merged so as to have the fixed one in the stochastic and Jeffery models, and keep the Renuda ones in the sphere model with modcpl=1, but this seems strange.

I had a look at this commit but I don't see any correction of the rotation matrix coefficients. Also, this rotation has been tested by two test cases and no errors have been found.

Also, the sphere model with modcpl = 1 and the stochastic model seem very similar, except for the way the orient_loc vector is computed (the axis is named main_axis in the Renuda model, and renamed singularity_axis from the Frenglish axe_singularity name in the POPART project's spheroid model, but that is probably just a choice of name).

This is correct. The only difference is the rotation direction. For modcpl = 1 the rotation is such so that the main axis for the calculations, i.e. the x-axis, becomes aligned with the average relative velocity vector. So, based on that requirement we compute the rotation angle 'rot_angle' and the rotation vector 'n_rot' (= cross product of the x-axis with the average relative velocity vector). The rotation matrix 'trans_m' is the same as with all types of rotations, the only difference are the values of 'rot_angle' and 'n_rot'. So, one could modify the code in order to compute these values for each case and then compute the rotation matrix with the same piece of code.

Given this, I wonder whether the modcpl and shape model options should be combined into a single keyword (at least based on the code in cs_lagr_sde.c, and assuming we sort out the rotation matrix orientation issues). I'll try to contact Jean-Pierre, as I believe he followed both developments.

I think that this is also possible.

I also need to automate the settings for idistu and idiffl, or base matching tests on modcpl, as these values do not seem to be modified automatically in the Renuda branch (or I can't find where).

Part of this pull request contains modifications to the Lagrangian part of the Code_Saturne GUI. The automatic setting of these flags happens at this level. In file gui/Pages/LagrangianView.py line 368 you can see the slot:

    @pyqtSlot()
    def slotMODCPL(self):
        """
        Input MODCPL.
        """
        if self.checkBoxMODCPL.isChecked():
            self.model.setRegularParticles(0)
            self.model.setTurbulentDispersion("off")
            self.model.setTurbulentDiffusion("on")
        else:
            self.model.setRegularParticles(1)
            self.model.setTurbulentDispersion("on")
            self.model.setTurbulentDiffusion("off")

So, the flags are adjusted automatically when the user creates the setup.xml file.

[YvanFournier]

The modified code breaks our "JET_2_PHI" test case.

Running the Renuda branch head, commit d4a8d7ad, (also under Valgrind), I see 2 errors, around line 1635 of cs_lagr_tracking.c:

• ustar is based on extra->ustar, but this points to NULL ustar (using k-epsilon, the "ustar" field is not always created).

• The update of particle_velocity_seen does not use a 3D vector expression (components are updated with different formulae).

[paspro]

The modified code breaks our "JET_2_PHI" test case.

Running the Renuda branch head, commit d4a8d7a, (also under Valgrind), I see 2 errors, around line 1635 of cs_lagr_tracking.c:

* `ustar` is based on `extra->ustar`, but this points to `NULL` ustar (using k-epsilon, the "ustar" field is not always created).

* The update of `particle_velocity_seen` does not use a 3D vector expression (components are updated with different formulae).

Regarding this piece of code, this is a modification of the wall boundary conditions for particles as specified by Jean-Pierre, If there are memory allocation issues I have not encountered before then we need to resolve them. The modifications involve the way the particle velocity seen is computed.

[YvanFournier]

Regarding "ustar", its presence may be forced (there is another example of this in cs_lagr_options), but Martin may have a more general solution, so I will wait for that (I discussed it briefly with him and hope to be able to solve this within a few days).

Regarding velocity_seen, Jean-Pierre's spec is incorrect, at least using general axes. The formula should be a true vector formula, not a projection a a specific system of axes (and in most places it is used, velocity seen seems to be in Cartesian, not local, coordinates). I sent an e-mail to Jean-Pierre, but did not receive any news so far, so I'm stuck here too.

I should be able to integrate the code once these 2 issues are resolved.

[YvanFournier]

Here are the latest patches relative to the current master branch.

The first is the same as before (rebase and edit of Renuda work), the second brings some improvements. But still waiting for input on the 2 issues described above (especially the one regarding the incorrect coordinate system).

0001-Lagrangian-module-enhancements-based-on-Renuda-devel.patch.txt 0002-Improvements-to-Lagrangian-settings-automation.patch.txt

[YvanFournier]

Hello,

Most of this work has been merged and rebased relative to the master branch in commit e53f1372.

I still suspect an error in cs_lagr_sde.c for the change of reference of the complete model. As mentioned before, the rotation expressions for CS_LAGR_SHAPE_SPHERE_MODEL were similar to those we had in the CS_LAGR_SHAPE_SPHEROID_STOC_MODEL case, prior to commit 1795093a which introduced a fix. Also, for the return transformation, the matrix was transposed twice : once when computing the coefficients, the second time when using the cs_math_33t_3_product function instead of cs_math_33_3_product as in the direct transformation.

I fixed the double transposition (along with code simplification) in follow-up commit d06d1042, but did not change the direct transformation, which should be checked.

Also, the new rebound model was not included:

• it depends on 0001-Lagrangian-wall-rebound-condition-based-on-Renuda-de.patch.txt u_star whether we are using a wall model or not (though this is a common assumption in the Lagrangian module)
• it is expressed in a local frame of reference while the "velocity" seen is used in the Cartesian from everywhere else in the code.
  • if the specs said this, the specs are wrong

[YvanFournier]

Though the previous comments are important to keep track of, since the essential part of this pull request has been integrated, we can close it now.