Fixing the turbulence advection order bug

[anilyil]

Purpose

This PR fixes the issue explained in #58. I do have a much deeper explanation of the issue and the fix below, so feel free to just read this PR. Furthermore, I renamed a number of turbulence solver routines, which I detail below. This renaming does not change any internal functionality, or user interface.

Type of change

• Bugfix (non-breaking change which fixes an issue)
• Code style update (formatting, renaming)

Summary of the issue

In ADflow, we use a first order discretization in space for the turbulence models by default. Optionally, users can select to use a second order accurate discretization for the turbulence model advection terms. This is done by setting 'turbulence order' to 'second order', where the default is 'first order' in the Python layer.

The option on the Python level is used to set the Fortran level integer-type variable orderturb found in inputDiscretization module in file https://github.com/mdolab/adflow/blob/master/src/modules/inputParam.F90#L15. Then, when the advection terms are actually computed, there is another variable that is used to determine if second order terms should be active. This logical-type variable is called secondOrd, and is found in https://github.com/mdolab/adflow/blob/master/src/turbulence/turbMod.F90#L10.

The idea with these two variables is that; orderturb is used to set the spatial discretization order on the final solution for the turbulence advection terms. However, because we also may use a full multigrid startup procedure, we want to only use these second order terms when the "ground level" is 1, i.e. we do not want to bother computing the second order accurate terms when we are going through the full multigrid startup procedure. I did not make this decision, and this is just my interpretation of what is happening in the code. As always, I can be wrong, and it is safest to assume I am wrong and verify this yourselves.

The bug fixed in this PR stems from this approach. The secondOrd variable was only set in two places in the code: https://github.com/mdolab/adflow/blob/master/src/turbulence/turbAPI.F90#L33-L38 and https://github.com/mdolab/adflow/blob/master/src/turbulence/turbAPI.F90#L126-L131. These lines are only run if we call these methods, and depending on the solver algorithms selected, they may never be run, so the value for the variable secondOrd may never be set, and we use the uninitialized value throughout the simulation and adjoint evaluation.

I have checked that the uninitialized value of this variable is .false. with both GNU and Intel compilers. Therefore, we should not have encountered a "consistency issue", where this variable would be determined at runtime, and possibly got different values on different processors. In both my docker and Pleiades installations, this variable has the value .false. whenever it is actually never set by our code. Furthermore, the same value of secondOrd is used for the adjoint solver, and therefore, regardless of the value, the gradients should have been accurate.

This PR fixes this issue, and with this PR, the variable secondOrd is always explicitly set. Below, I list the okay and not-okay scenarios, where with the former ones, the variable is set correctly, and with the latter, the variable is not set correctly. After this PR, both scenarios will have the variable set correctly.

Okay Cases

This variable is set in the two places I mentioned above. If the simulation starts with a smoother (i.e. RK or DADI), or with decoupled ANK and the turbDADI solver in ADflow (i.e. not the turbKSP in ANK), the solver does call the method for solving the turbulence model equations using a DDADI approach. This call sets this value, and the correct value is propagated throughout in most cases.

Not Okay Cases

If the solver doest not call the DDADI based turbulence solver, the value of secondOrd will never be set explicitly. This would happen if we start with:

1. Coupled ANK
2. turbKSP solver that implements a decoupled ANK method to the turbulence
3. NK

If any of these combinations are used to start the solver (NK start might happen after a restart from a converged state), the secondOrd variable is never set, and the solver defaults to the uninitialized value provided by the compiler. Like I said, this was .false. for GNU and Intel. Furthermore, all of these are scenarios require non-default options, except for restarting straight into the NK solver. Even then, almost nobody uses second order turbulence, so the correct value is used.

On a similar note, because this variable is not explicitly set with these scenarios, if the user wanted to use second order accurate advection terms by setting 'turbulence order' to 'second order', the code would not change anything and they would still get first order turbulence. This is how I originally discovered this bug.

Fix

This PR fixes the issue above by changing where the secondOrd variable is set. This variable is a very low-level variable that is only used inside the nested loops to compute the advection terms. These loops exist in two places in the code:

1. Block residuals: https://github.com/mdolab/adflow/blob/master/src/turbulence/turbUtils.F90#L826
2. Blockette residuals: https://github.com/mdolab/adflow/blob/master/src/NKSolver/blockette.F90#L1391

Because this variable is only used in these routines, I thought the best place to set the variable would be inside these routines, so that we do not miss any edge cases. I have made these modifications in these two places: https://github.com/anilyil/adflow/blob/turb_fix/src/turbulence/turbUtils.F90#L871-L875 and https://github.com/anilyil/adflow/blob/turb_fix/src/NKSolver/blockette.F90#L1407-L1411. Because of this modification, I re-ran tapenade and the all three AD versions of turbAdvection function from turbUtils are updated accordingly.

Another required change as a result of this approach: We used to save the value in secondOrd and set it to .false. explicitly whenever we wanted to get first order advection routines. These are with first order ANK routines, and with preconditioner creation routines. Now, instead of saving and setting secondOrd, ADflow saves and sets the variable orderturb directly, which is then used in the two low-level routines I listed above. After we are done with the first order routines, we set back the saved value of orderturb. These changes can be found by looking at the diff for the following files in this PR:

1. https://github.com/mdolab/adflow/blob/master/src/NKSolver/NKSolvers.F90
2. https://github.com/mdolab/adflow/blob/master/src/adjoint/adjointUtils.F90
3. https://github.com/mdolab/adflow/blob/master/src/adjoint/fortranPC.F90

As a result of these modifications, this PR closes #58

Separate change of turbulence solver naming

The DDADI based turbulence solver in ADflow has always been referred to as a "segregated solver" in the code. I disliked this name from the beginning, and in this PR, I changed all instances of the word "segregation" to the more accurate and appropriate versions like decoupled, or the actual algorithm, which is DDADI.

Testing

All regression tests pass. I also tested if the second order advection terms are activated properly with all the solver options listed above.

Checklist

• I have run unit and regression tests which pass locally with my changes

[sseraj]

First off, thanks for the detailed description. It explained the changes really well.

I took a look through the code and the changes make sense:

1. secondOrd is only used in the block and blockette residual routines
2. setting secondOrd explicitly inside these routines addresses the bug
3. by saving and setting turborder instead of (the unset) secondOrd outside of these routines, subsequent calls to the residual routines have the intended order

[joanibal]

Nice work