Adding temperature limits for incompressible solver

[Cristopher-Morales]

Proposed Changes

The aim of this pull request is to add temperature limits to prevent unphysical values for incompressible flows. Controlling temperature values is especially important for multicomponent and reacting flows, where low and high values of temperature can lead to unbounded density values, which might result in unphysical species mass fractions values (outside range 0.0-1.0).

Test H2-Air mixture adapted from test case rectangle passive transport validation:

top: develop
bottom: using temperature limits

Species mass fraction H2:

Density:

Related Work

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

[pcarruscag]

Have you considered turning on the automatic underrelaxation for the incompressible solver?

void CIncEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) {

  CompleteImplicitIteration_impl<false>(geometry, config);
}

I think we are not using it because the pressure and temperature can be negative in incompressible, but if there are physical limits it's easier to do.

I'm not saying "do this instead", maybe these custom limits can complement the auto relaxation (which feeds the CFL adaptation too).

[Cristopher-Morales]

Have you considered turning on the automatic underrelaxation for the incompressible solver?

void CIncEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) {

  CompleteImplicitIteration_impl<false>(geometry, config);
}

I think we are not using it because the pressure and temperature can be negative in incompressible, but if there are physical limits it's easier to do.

I'm not saying "do this instead", maybe these custom limits can complement the auto relaxation (which feeds the CFL adaptation too).

Hi!

I have tried the automatic underrelaxation for the incompressible solver and it helps. However, there might be some issues for future developments if the enthalpy for reacting flows (which range from positive to negative values) can be solved using the preconditioning.
This is because I observed that the Underrelaxation factor computation does not allow the energy variable to change sign, in particular, when the energy variable (enthalpy) is close to zero.
I added this ComputeUnderRelaxationFactor in other branch when I am solving for the enthalpy in the preconditioning:

void CIncEulerSolver::ComputeUnderRelaxationFactor(const CConfig* config) {

  /* Loop over the solution update given by relaxing the linear
   system for this nonlinear iteration. */

  const su2double allowableRatio = 0.2;

  SU2_OMP_FOR_STAT(omp_chunk_size)

  for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
    su2double localUnderRelaxation = 1.0;

    /*--- Energy ---*/

    const unsigned long index = iPoint * nVar + nVar - 1;

    su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, nVar - 1)) + EPS);

    if (ratio > allowableRatio) {
      localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
    }

    /* Threshold the relaxation factor in the event that there is
     a very small value. This helps avoid catastrophic crashes due
     to non-realizable states by canceling the update. */

    if (localUnderRelaxation < 1e-3) **localUnderRelaxation = 0.001**;

    /* Store the under-relaxation factor for this point. */

    nodes->SetUnderRelaxation(iPoint, localUnderRelaxation);
  }

  END_SU2_OMP_FOR
}

I had to add a minimum localUnderRelaxation = 0.001 to allow the enthalpy to change from positive to negative. This is because this quantity: allowableRatio / ratio can be quite small, and if I keep the threshold used by Euler and Nemo Solvers:

if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0;

The enthalpy does not change from positive to negative.

Please let me know what you think about this. There might be a better criterium that can help with this issue.

Thanks in advance!!!

[Cristopher-Morales]

Hi,

I have enabled the automatic under-relaxation for the incompressible solver, for now It is applied to the energy equation and it is skipped when the energy equation is off.
I am using the option GetMaxUpdateFractionFlow() which is used for the compressible solver. I am not sure if I should add one specific option for the incompressible solver or just add the option MAX_UPDATE_FLOW=1.0 to the existing test cases in SU2 in order to avoid breaking the current test cases
Please let me know what you think about this.
Thanks in advance!!!

[pcarruscag]

Nice. It is fine to reuse the option, are the limits still playing a big role after adding the relaxation?

[Cristopher-Morales]

Nice. It is fine to reuse the option, are the limits still playing a big role after adding the relaxation?

Hi!

Thanks for the comment.

Yes, the limits are still playing a role. The automatic under relaxation helps to delay that the solver diverges, but the temperature reaches values close to zero that causes an increase of density and finally a divergence in the species solver. I can add the est case that I am using for testing this feature.

Top : relaxation without temperature limits
Bottom: relaxation with temperature limits

Density:

Species:

Please let me know what you think about this.

Thank you so much in advance!!!

[bigfooted]

Hi Cristopher, good to add this in! I would not add the underrelaxation now if it is still unclear how to make it work in a robust way. But good to also keep that on the to-do list (give it to a student project maybe).
Why are so many regression tests affected, do they all hit the boundary of the temperature limits?
Does it make sense to choose as lower limit not zero but a positive value, maybe 250K or something?

[Cristopher-Morales]

Hi Cristopher, good to add this in! I would not add the underrelaxation now if it is still unclear how to make it work in a robust way. But good to also keep that on the to-do list (give it to a student project maybe). Why are so many regression tests affected, do they all hit the boundary of the temperature limits? Does it make sense to choose as lower limit not zero but a positive value, maybe 250K or something?

Hi Nijso,
Thanks for your comment!
Yes, that could also be an option to propose as a student project to investigate how to mae the automatic underrelation approach more robust, which could eventually be very useful for the species solver as well.
The regression test that are failing are due to enable the automatic underrelaxation where the option MAX_UPDATE_FLOW has a default value 0.2. I think if I add to the failing test cases the option MAX_UPDATE_FLOW=1.0, they should not fail. When I disable the automatic underrelaxation, the regression tests do not fail.
It does make sense to choose a lower limit 250K, in the figure of species and density that I obtained using relaxation with temperature limits (bottom figures), the temperature limits are 250,400.

[bigfooted]

Ok thanks, so with underrelaxation for the energy, is convergence improving as well for cases that were not problematic?

[pcarruscag]

We don't know how many people out there are using the incompressible solver for heat transfer, so I would advocate for default values that keep the default behavior.
Thanks for checking the underrelaxation, if most of the improvements still come from the limits I think it's fine to leave it for a project.

[Cristopher-Morales]

We don't know how many people out there are using the incompressible solver for heat transfer, so I would advocate for default values that keep the default behavior. Thanks for checking the underrelaxation, if most of the improvements still come from the limits I think it's fine to leave it for a project.

Thanks for the comment.

Yes, the most of the improvements come from the limits. I am going to disable the ComputeUnderRelaxationFactor in this pull request and let its implementation for a future project.

[bigfooted]

Hi Cristopher, can you add this testcase as well?

[Cristopher-Morales]

Hi Cristopher, can you add this testcase as well?

Hi Nijso,

Yes, I will add this test case as soon as possible.

Thanks in advance!!!