Make two point continuation more robust

[speth]

Changes proposed in this pull request

• Eliminate the time-dependent terms in the energy equation at the temperature control points. These terms were inherently inconsistent with the equation forcing the temperature to equal the control set point.
• Eliminate the time-dependent terms in the radial momentum equation when two point control. This allows the radial velocity gradient at the temperature control points to respond "instantaneously" to changes in the boundary mass flux and therefore allowing changes in the inlet mass flow rates to satisfy the equations fixing the temperature at the control points. Previously, the solver would almost always fail during time stepping.
• Ensure that internal arrays (such as Flow1D::m_rho) are up to date after a regridding operation that only removes points and where the flame is not re-solved. Previously, this could result in the wrong value of the density being stored in an output SolutionArray and incorrect values being calculated for other derived quantities like heat release rate.
• Correct order of operations when restoring a flame where two-point control was active so the value of Uo is set correctly
• Revise the two point continuation example to demonstrate a more efficient and fault tolerant algorithm for traversing the stable and unstable branches of the curve.

If applicable, fill in the issue number this pull request is fixing

Fixes #1747

If applicable, provide an example illustrating new features this pull request is introducing

Checklist

• [x] The pull request includes a clear description of this code change
• [x] Commit messages have short titles and reference relevant issues
• [x] Build passes (scons build & scons test) and unit tests address code coverage
• [x] Style & formatting of contributed code follows contributing guidelines
• [x] The pull request is ready for review

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 78.94737% with 4 lines in your changes missing coverage. Please review.

Project coverage is 73.21%. Comparing base (7aec821) to head (e20ec14). Report is 23 commits behind head on main.

Files	Patch %	Lines
src/oneD/Sim1D.cpp	70.00%	2 Missing and 1 partial :warning:
src/oneD/Flow1D.cpp	88.88%	1 Missing :warning:

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #1779 +/- ## ======================================== Coverage 73.20% 73.21% ======================================== Files 381 381 Lines 54371 54249 -122 Branches 9248 9239 -9 ======================================== - Hits 39803 39716 -87 + Misses 11590 11563 -27 + Partials 2978 2970 -8 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[wandadars]

This is a major improvement to the two-point control implementation. Thank you for the work on tracking down these issues @speth !

Have you looked at the heat release rate quantity between successive calls to the solve() with the two-point control method? I'm seeing the values jump around quite a bit, but when I look at the flow variables, they all seem fairly steady-ish between the calls to solve() with the exception of the temperature decrement. No large variations seem to be present, yet the heat release seems to change quite a bit. Below is two plots of the heat release for two solutions that differ by 5K at the control points.

At higher strain rates, the noisy-ness seems to go away. Or at least is hidden by the much higher heat release rates (starting at 10e8 and going to 10e14 as strain increases, for the H2/O2 case I'm looking at.)

The quantity is computed as:

- np.sum(self.partial_molar_enthalpies * self.net_production_rates, 0)

I'll print these out and look at which one is changing between the solutions later today.

[ischoegl]

@speth ... I am mostly waiting for a response from you to @wandadars' comment. Other than that, I am ok with approving.

[speth]

@wandadars, I had seen a few cases with funky looking heat release rate profiles early on, which I believe was related to incorrect density information being used in some cases and was fixed by d671f0faa6. With the current branch, looking at an H2/O2 flame at a few different pressures, the most negative heat release rate values I see are around -1e5 W/m^3, and I don't see any with the extreme spikes of your first graph. Can you share what conditions those are occurring under?

[wandadars]

@speth It might be something with my script that is testing the two-point method. For reference, this is the case I was looking at:

# Input parameters
p: 5.6e6  # pressure [Pascals]
fuel_inlet_temp: 800.0  # fuel inlet temperature [Kelvin]
oxidizer_inlet_temp: 711.0  # oxidizer inlet temperature [Kelvin]
fuel_mdot: 0.1 # kg/m^2/s
width: 50.0e-3 # Distance between inlets [m]

# Numerics
steady_tolerances:
    relative_tolerance: 1e-4
    absolute_tolerance: 1e-7

refine_criteria:
    ratio: 8
    slope: 0.1
    curve: 0.2
    prune: 0.03

#Mass Fraction composition
oxidizer_composition: 'O2:1.0'  # oxidizer composition (right)
fuel_composition: 'H2:1.0'  # fuel composition (left)

# Reaction mechanism
mechanism: 'h2o2.yaml'

# Two-Point Controls
spacing: 0.95
initial_temperature_increment: 5
maximum_temperature_increment: 50
maximum_change_in_max_temperature: 50  # This prevents the maximum flame temperature from dropping too fast

[speth]

I think this is a result of your very loose absolute tolerance solver interacting with reactions that are very nearly at equilibrium. I notice this mainly with the reaction H2 + OH <=> H + H2O, where the spikes in heat release come after a case where the solver regrids, adds a point for one of the species in this reaction, and then makes an immediately successful Newton solve. In such a case, it's possible to satisfy the Newton step criteria even though the species values at the new point could be closer to the true steady state, and this is all amplified by a reaction like this one which is nearly at equilibrium, with very high forward and reverse rates.

So, I don't think there's anything to address that's related to the two point solver, or that should be handled as part of this PR.

[wandadars]

I think this is a result of your very loose absolute tolerance solver interacting with reactions that are very nearly at equilibrium. I notice this mainly with the reaction H2 + OH <=> H + H2O, where the spikes in heat release come after a case where the solver regrids, adds a point for one of the species in this reaction, and then makes an immediately successful Newton solve. In such a case, it's possible to satisfy the Newton step criteria even though the species values at the new point could be closer to the true steady state, and this is all amplified by a reaction like this one which is nearly at equilibrium, with very high forward and reverse rates.

So, I don't think there's anything to address that's related to the two point solver, or that should be handled as part of this PR.

I ran with the default tolerances and for completeness, this is the profile as the temperature is marched down.

[speth]

That's a nice visualization. Can you share the full script for generating those results and plots? You could also consider posting that with a new issue. While I think it's separate from the issues being addressed here, I think there is room for improving the heat release profiles.

[wandadars]

Here's the script that generated the gif.

from pathlib import Path
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
import imageio
from pathlib import Path

p = 5.4e6 # pressure [Pascals]

fuel_inlet_temp= 800.0  # fuel inlet temperature [Kelvin]
oxidizer_inlet_temp= 711.0  # oxidizer inlet temperature [Kelvin]
fuel_mdot= 0.5 # kg/m^2/s
width= 50.0e-3 # Distance between inlets [m]

oxidizer_composition= 'O2:1.0'  # oxidizer composition (right)
fuel_composition= 'H2:1.0'  # fuel composition (left)

# Reaction mechanism
mechanism= 'h2o2.yaml'

# Define output locations
output_path = Path() / "testing_output"
output_path.mkdir(parents=True, exist_ok=True)

gas = ct.Solution(mechanism)

gas.TPY = fuel_inlet_temp, p, oxidizer_composition
density_f = gas.density

gas.TPY = oxidizer_inlet_temp, p, oxidizer_composition
density_o = gas.density

#Unity Lewis number testing
gas.transport_model = 'unity-Lewis-number'

f = ct.CounterflowDiffusionFlame(gas, width=width)

f.set_refine_criteria(ratio=15, slope= 0.15, curve= 0.08, prune= 0.04)
#Set the state of the two inlets
f.fuel_inlet.mdot = fuel_mdot
f.fuel_inlet.X = fuel_composition
f.fuel_inlet.T = fuel_inlet_temp

#Create a guestimate for what the oxidizer mdot would be
f.oxidizer_inlet.mdot = (fuel_mdot / density_f) * density_o*4
f.oxidizer_inlet.X = oxidizer_composition
f.oxidizer_inlet.T = oxidizer_inlet_temp

# Generate initial condition
f.solve(auto=True, loglevel=4)
f.save('backup.yaml', name="solution", overwrite=True)

print('mdot info:')
print('Fuel mdot: ' + str(f.fuel_inlet.mdot))
print('Oxidizer mdot: ' + str(f.oxidizer_inlet.mdot))
print('BC State:')
print('Left U: ' + str(f.velocity[0]))
print('Right U: ' + str(f.velocity[-1]))

print('Starting two-point control')
f.two_point_control_enabled = True

spacing = 0.75
temperature_increment = 10 # Initial temperature increment
maximum_temperature = []
a_max = []
filenames = []
for i in range(80):
    print('Iteration: ' + str(i+1))
    print('mdot info:')
    print('Fuel mdot: ' + str(f.fuel_inlet.mdot))
    print('Oxidizer mdot: ' + str(f.oxidizer_inlet.mdot))

    print('BC State:')
    print('Left U: ' + str(f.velocity[0]))
    print('Right U: ' + str(f.velocity[-1]))

    control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T))
    print(f'Iteration {i}: Control temperature = {control_temperature} K')
    print(f'Increment size: {temperature_increment} K')
    f.set_left_control_point(control_temperature)
    f.set_right_control_point(control_temperature)
    # This decrement is what drives the two-point control. If failure
    # occurs, try decreasing the decrement.
    f.left_control_point_temperature -= temperature_increment
    f.right_control_point_temperature -= temperature_increment

    try:
        f.solve(loglevel=0)
        #f.solve(loglevel=0,refine_grid=False)
    except ct.CanteraError as e:
        print('Solver did not converge.')
        print('Error:')
        print(e)

        if temperature_increment < 1:
            print('Smallest increment reached, stopping')
            break

        f.left_control_point_temperature += temperature_increment
        f.right_control_point_temperature += temperature_increment
        temperature_increment *= 0.5
        continue

    # Increase the increment if the solution worked out well
    if temperature_increment < 50:
        temperature_increment *= 1.1

    maximum_temperature.append(np.max(f.T))
    a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))

    # Plot the heat release rate as well as the locations of the control points using the
    # f.left_control_point_coordinate and f.right_control_point_coordinate attributes. Plot
    # The locations as vertical lines. Generate a gif of the heat release rate as a function
    # of position with the markers and save it as heat_release_rate.gif
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))  # 16x8 is a fairly large and wide size
    # Set the title as having the iteration number in it
    fig.suptitle(f"Iteration {i}")
    # Plot the temperature on the top subplot
    ax1.plot(f.grid, f.T, label='Temperature')
    ax1.axvline(f.left_control_point_coordinate, color='r', linestyle='--')
    ax1.axvline(f.right_control_point_coordinate, color='r', linestyle='--')
    ax1.set_xlabel('Position [m]')
    ax1.set_ylabel('Temperature [K]')
    ax1.legend()

    # Contrain the x axis to be between 0.02 and 0.03
    ax1.set_xlim(0.02, 0.03)

    # Plot the heat release rate on the bottom subplot
    ax2.plot(f.grid, f.heat_release_rate, label='Heat Release Rate', color='b')
    ax2.axvline(f.left_control_point_coordinate, color='r', linestyle='--')
    ax2.axvline(f.right_control_point_coordinate, color='r', linestyle='--')
    ax2.set_xlabel('Position [m]')
    ax2.set_ylabel('Heat Release Rate [W/m^3]')
    ax2.legend()
    ax2.set_xlim(0.02, 0.03)

    # Save each frame as a separate image file
    frame_filename = output_path / f"heat_release_rate_temp_{i}.png"
    plt.savefig(frame_filename)
    plt.close(fig)  # Close the figure to avoid memory leaks

    # Append the filename to the list
    filenames.append(frame_filename)

# Generate the gif(duration is in hundreths of a second per frame)
with imageio.get_writer(output_path / "heat_release_rate.gif", mode='I', duration=500, loop=0) as writer:
    for filename in filenames:
        image = imageio.imread(filename)
        writer.append_data(image)

# Plot the maximum temperature versus the maximum axial velocity gradient
plt.figure()
plt.plot(a_max, maximum_temperature)
plt.xlabel('Maximum Axial Velocity Gradient [1/s]')
plt.ylabel('Maximum Temperature [K]')
plt.savefig(output_path / "figure_max_temperature_vs_max_velocity_gradient.png")

# Plot maximum_temperature against number of iterations
plt.figure()
plt.plot(range(len(maximum_temperature)), maximum_temperature)
plt.xlabel('Number of Continuation Steps')
plt.ylabel('Maximum Temperature [K]')
plt.savefig(output_path / "figure_max_temperature_iterations.png")