Counterflow Flame control

[wandadars]

This pull request supersedes #1510 and is based on the newer coding in the StFlow class. This feature aims to provide a capability for users to generate 1D flame solutions that include solutions along the unstable burning branch of the typical s-curve that is used for flamelet progress variable methods.

• Introduce two-point flame control method that adjusts the Lambda equation and adds a new equation, Uo as a dummy equation.
• C++ adjustments to the StFlow class for axisymmetric stagnation flow configurations (only valid case for the two-point control currently)
• Python interface to allow users to control the two-point control.

Checklist

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

[wandadars]

For reference, this is the control file that I've been using for testing the overall execution of the two-point flame control method. It's still not very stable, which I believe it should at least be somewhat stable based on the lack of any caveats present in the 1994 paper from Nishioka, Law, and Takeno (the paper that serves as the basis for this implementation).

Also attached is the mechanism file I am using to debug the implementation. I had to rename the YAML file to have a .txt suffix to make Github happy.
burke.txt

import os
import time

import matplotlib.pyplot as plt
import numpy as np

import cantera as ct


def plot_solution(f, solution_num, control_points=False):
    variables = [('temperature', f.T, 'Temperature (K)'), ('Uo', f.Uo, 'Oxidizer Velocity (m/s)'),
                 ('L', f.L, 'Lambda (m/s)'), ('u', f.velocity, 'Axial Velocity (m/s)'),
                 ('V', f.spread_rate, 'Spread Rate (1/s)')
                ]

    for folder, data, ylabel in variables:
        fig, ax = plt.subplots()
        ax.plot(f.grid, data, 'o-')
        if control_points and folder == 'temperature':
            plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
            plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
            ax.legend()
        ax.set_ylabel(ylabel)
        plt.savefig(f'output/{folder}/flamelet_{solution_num}.png')
        plt.close()

    # Combined plot
    fig, axs = plt.subplots(3, 2, figsize=(10,10))
    for i, (folder, data, ylabel) in enumerate(variables + [('temperature_zoom', f.T, 'Temperature (K) Zoomed')]):
        ax = axs.flat[i]
        ax.plot(f.grid, data, 'o-')
        ax.set_xlabel('Distance, m')
        ax.set_ylabel(ylabel)
        if folder == 'temperature_zoom':
            high_temp_indices = np.where(data > 310)[0]
            if len(high_temp_indices) > 0:
                zoom_range_min = f.grid[high_temp_indices[0]]
                zoom_range_max = f.grid[high_temp_indices[-1]]
                ax.set_xlim([zoom_range_min, zoom_range_max])
                plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
                plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
                ax.legend()

    # Add a title to the figure
    fig.suptitle(f'Flamelet Number: {solution_num}', fontsize=16)
    plt.subplots_adjust(wspace=0.4, hspace=0.6)
    plt.savefig(f'output/combined/flamelet_{solution_num}.png')
    plt.close()

def plot_mixture_fraction_solution(f, solution_num, control_points=False):
    variables = [('temperature_z', f.T, 'Temperature (K)'), ('Uo_z', f.Uo, 'Oxidizer Velocity (m/s)'),
                 ('L_z', f.L, 'Lambda (m/s)'), ('u_z', f.velocity, 'Axial Velocity (m/s)'),
                 ('V_z', f.spread_rate, 'Spread Rate (1/s)')
                ]

    # Get the practical range for the mixture fraction for plotting only valuable data
    epsilon = 1e-3
    grid = f.mixture_fraction('Bilger')
     # Find indices where grid crosses epsilon threshold
    greater_than_epsilon_indices = np.where(grid > epsilon)[0]
    less_than_epsilon_indices = np.where(grid < 1 - epsilon)[0]
    
    if greater_than_epsilon_indices.size > 0 and less_than_epsilon_indices.size > 0:
        first_index = greater_than_epsilon_indices[0]
        last_index = less_than_epsilon_indices[-1]
        subset_indices = range(first_index, last_index + 1)
    else:
        print("No indices found for the specified epsilon threshold")
        return

    for folder, data, ylabel in variables:
        fig, ax = plt.subplots()
        ax.plot(grid[subset_indices], data[subset_indices], 'o-')
        if control_points and folder == 'temperature':
            plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
            plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
            ax.legend()
        ax.set_ylabel(ylabel)
        ax.set_xlabel('Mixture Fraction(Bilger), Z')
        plt.savefig(f'output/{folder}/flamelet_{solution_num}.png')
        plt.close()

    # Combined plot
    fig, axs = plt.subplots(3, 2, figsize=(10,10))
    for i, (folder, data, ylabel) in enumerate(variables + [('temperature_control', f.T, 'Temperature (K)')]):
        ax = axs.flat[i]
        ax.plot(grid[subset_indices], data[subset_indices], 'o-')
        ax.set_xlabel('Mixture Fraction(Bilger), Z')
        ax.set_ylabel(ylabel)
        if folder == 'temperature_control':
            plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
            plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
            ax.legend()

     # Add a title to the figure
    fig.suptitle(f'Flamelet Number: {solution_num}', fontsize=16)
    plt.subplots_adjust(wspace=0.4, hspace=0.6)
    plt.savefig(f'output/combined_z/flamelet_{solution_num}.png')
    plt.close()

def create_output_folders():
    folders = ['temperature', 'Uo', 'L', 'u', 'V', 'combined','temperature_z',
               'Uo_z', 'L_z', 'u_z', 'V_z', 'combined_z', 'flamelets']
    for folder in folders:
        path = f'output/{folder}'
        if not os.path.isdir(path):
            os.makedirs(path, exist_ok=True)  # exist_ok=True will prevent any error if the directory exists.
        
#Callback testing
def plot_variable(grid, data, title, xlabel, ylabel, filename):
    plt.figure()
    plt.plot(grid, data, 'o-', lw=2, label=ylabel)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.legend()
    plt.savefig(filename)
    plt.close()
    
def steady_callback(t, f):
    variables = [
        (f.T, 'Temperature (K)', 'temperature_plot_at_'),
        (f.velocity, 'Axial Velocity (m/s)', 'axial_velocity_plot_at_'),
        (f.L, 'Lambda', 'lambda_plot_at_'),
        (f.Uo, 'Oxidizer Velocity Uo (m/s)', 'oxidizer_velocity_plot_at_')
    ]
    for data, ylabel, prefix in variables:
        plot_variable(f.grid, data, f'State at t={t}', 'Grid Points', ylabel, f'{prefix}{t}.png')
    
    #input("Press Enter to continue...")  # Pause for user input before continuing
    return 0

def unsteady_callback(t, f):
    variables = [
        (f.T, 'Temperature (K)', 'temperature_plot_unsteady'),
        (f.velocity, 'Axial Velocity (m/s)', 'axial_velocity_plot_unsteady'),
        (f.L, 'Lambda', 'lambda_plot_unsteady'),
        (f.Uo, 'Oxidizer Velocity Uo (m/s)', 'oxidizer_velocity_plot_unsteady')
    ]
    for data, ylabel, prefix in variables:
        plot_variable(f.grid, data, f'State at t={t}', 'Grid Points', ylabel, f'{prefix}.png')
    
    #input("Press Enter to continue...")  # Pause for user input before continuing
    return 0
    
def plot_process(process_data):
    #Plot process
    
    # 1/a is often used as the parameter to plot flame solution data against in literature
    inverse_strain = [1.0/a for a in flamelet_process_plot['a_mean']]
    
    fig, ax = plt.subplots()
    ax.plot(inverse_strain, flamelet_process_plot['T_max'], 'o-')
    plt.savefig('T_max_vs_strainRate.png')

    fig, ax = plt.subplots()
    ax.plot(inverse_strain, flamelet_process_plot['Q'], 'o-')
    plt.savefig('Q_vs_strainRate.png')

    fig, ax = plt.subplots()
    ax.plot(np.diff(flamelet_process_plot['T_max']), 'o-')
    plt.savefig('Temperature_change.png')
    
    
#-----------------------

def setup_flame(configs):
    p = configs['p']
    fuel_inlet_temp = configs['fuel_inlet_temp']
    oxidizer_inlet_temp = configs['oxidizer_inlet_temp']
    fuel_mdot = configs['fuel_mdot']
    width = configs['width']
    oxidizer_composition = configs['oxidizer_composition']
    fuel_composition = configs['fuel_composition']
    mechanism = configs['mechanism']
    
    gas = ct.Solution(mechanism)

    # Set the composition using mass fractions
    gas.Y = oxidizer_composition
    # Retrieve and print the mole fractions for the oxidizer
    oxidizer_composition = gas.X
    print("Oxidizer composition in mole fractions:", oxidizer_composition)

    # Repeat for the fuel composition
    gas.Y = fuel_composition
    # Retrieve and print the mole fractions for the fuel
    fuel_composition = gas.X
    print("Fuel composition in mole fractions:", fuel_composition)
    
    gas.TPX = fuel_inlet_temp, p, fuel_composition
    density_f = gas.density

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

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

    # 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 * 0.5  
    f.oxidizer_inlet.X = oxidizer_composition
    f.oxidizer_inlet.T = oxidizer_inlet_temp

    #Solver settings
    f.set_refine_criteria(ratio=10, slope=0.6, curve=0.6, prune=0.05)

    # Adjust solver tolerances
    f.flame.set_transient_tolerances(default=(1e-8, 1e-12))
    f.flame.set_steady_tolerances(default=(1e-8, 1e-12))
    
    # Attach callback functions to the solver
    if callable(configs['steady_callback']):
        steady_callback = configs['steady_callback']
        f.set_steady_callback(lambda t: steady_callback(t, f))
    
    #if callable(configs['unsteady_callback']):
    #    unsteady_callback = configs['unsteady_callback']
    #    f.set_time_step_callback(lambda t: unsteady_callback(t, f))

    
    return f

def generate_flamelets(configs):
    print('Solving first solution with no two-point control active')
    loglevel = configs['loglevel']  
    f = setup_flame(configs)
    
    # Solve the initial problem with the standard counterflow solver settings
    f.solve(loglevel)
    print('Current Flamelet State:')
    print('Flame: Tmin: ' + str(np.min(f.T)) + ', Tmax: ' + str(np.max(f.T)))

    # Dict for storing useful data during execution
    flamelet_process = {'T_max': [max(f.T)]}
    flamelet_process['Q'] = [np.trapz(f.heat_release_rate, f.grid)]
    flamelet_process['a_mean'] = [f.strain_rate('mean')]
    
    #Point control
    f.two_point_control_enabled = True
    spacing = configs['spacing']
    maximum_temperature_increment = configs['maximum_temperature_increment']
    Tmax_threshold = 25
    temperature_increment_growth_rate = configs['temperature_increment_growth_rate']
    failed_previous_attempt = False
    temperature_increment = configs['temperature_increment'] #Kelvin
    flamelet_solution_num = 1
    iterations = 0
    create_output_folders()
    while ((np.max(f.T)-np.min(f.T))/np.max(f.T)) > 0.005 and iterations < 1000:        
        control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T)) 
        print('New Control Temperature: ' + str(control_temperature))
        f.set_left_control_point(control_temperature) #Sets a location based on the temperature of the current flame
        f.set_right_control_point(control_temperature)
        print(f'Current Flamelet {flamelet_solution_num} State:')
        print('Flame: Tmin: ' + str(np.min(f.T)) + ', Tmax: ' + str(np.max(f.T)))
        print('tLeft= ' + str(f.tLeft) + ' , tRight= ' + str(f.tRight))
        print('ZLeft= ' + str(f.zLeft) + ' , ZRight= ' + str(f.zRight))
        print('Solving flowfield without any temperature decrement to adjust for shifting temperature control points')
        f.solve(loglevel)
        if iterations == 0:
            iterations +=1
            continue # We want to run this first part onces using two-point algorithm with no temperature decrements before beginning to march
        f.save(f'output/flamelets/flamelet_{flamelet_solution_num}.yaml',overwrite=True) 
        f.save('previous_solution.yaml',overwrite=True) #Backup in case of failure
        # Lower the control point temperatures
        f.tLeft -= temperature_increment
        f.tRight -= temperature_increment
        print(f'Attempting to solve flame with: TLeft = {f.tLeft}, TRight= {f.tRight}, Tincrement= {temperature_increment}')
        iterations +=1
        try:
            f.solve(loglevel)
        except Exception as err:
            print(f'\nerror when using control points: TLeft = {f.tLeft}, TRight= {f.tRight}')
            print(err)

            # Reverse failed temperature control point decrement
            f.tLeft += temperature_increment
            f.tRight += temperature_increment
            print(f'Resetting solution to previous state: TLeft = {f.tLeft}, TRight= {f.tRight}')
            f.restore('previous_solution.yaml') #Backup in case of failure
           
            control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T)) 
            print('Control Temperature: ' + str(control_temperature))
            f.set_left_control_point(control_temperature) 
            f.set_right_control_point(control_temperature)
            print(f'Attempting to solve flame with: TLeft = {f.tLeft}, TRight= {f.tRight}, Tincrement= {temperature_increment}')
            f.solve(loglevel)
            print('Successfully restored state & shifted control points' )
            print('tLeft= ' + str(f.tLeft) + ' , tRight= ' + str(f.tRight))
            print('ZLeft= ' + str(f.zLeft) + ' , ZRight= ' + str(f.zRight))
            
            #Update increment
            print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
            temperature_increment *= 0.5
            if(temperature_increment < 5.0e-2):
                print('Smallest temperature increment reached, exiting.')
                break
            failed_previous_attempt = True
            continue
         
        '''
        if abs(np.max(f.T) - previous_Tmax) > Tmax_threshold: #Do not allow a large maximum temperature jump between solutions
            print('Maximum temperature change violated, ajusting step size')
            print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
            #Update increment
            temperature_increment *= 0.5
            if(temperature_increment < 5.0e-1):
                spacing *= 1.01
                print('No more solutions reached at temperature increment: ' + str(temperature_increment) + ', adjusting spacing to: ' + str(spacing))
                temperature_increment = original_temperature_increment # Reset the increment
            failed_previous_attempt = True
            f.restore('previous_solution.yaml')
            continue
        '''
        if failed_previous_attempt == False:
            if temperature_increment < maximum_temperature_increment:
                temperature_increment *= temperature_increment_growth_rate

        failed_previous_attempt = False
        
        flamelet_process['T_max'].append(max(f.T))
        flamelet_process['Q'].append(np.trapz(f.heat_release_rate, f.grid))
        flamelet_process['a_mean'].append(f.strain_rate('mean'))
        print(f"Flame Information: {flamelet_process['T_max'][-1]:.4g} K (a = {flamelet_process['a_mean'][-1]:.4g}, Q = {flamelet_process['Q'][-1]:.4g})")
        
        plot_solution(f, flamelet_solution_num, control_points=True)
        plot_mixture_fraction_solution(f, flamelet_solution_num, control_points=True)
        flamelet_solution_num += 1


    if iterations == 1000:
        print('Exited due to maximum iterations being exceeded')

if __name__ == '__main__':
    configs = {}
    #'''
    # Input parameters
    configs['p'] = 8.0e5  # pressure [Pascals]
    configs['fuel_inlet_temp'] = 1200.0  # fuel inlet temperature [Kelvin]
    configs['oxidizer_inlet_temp'] = 1200.0  # oxidizer inlet temperature [Kelvin]
    configs['fuel_mdot'] = 0.2 # kg/m^2/s
    configs['width'] = 20e-3 # Distance between inlets [m]
    configs['loglevel'] = 1  # amount of diagnostic output (0 to 9)

    #Mass Fraction composition
    configs['oxidizer_composition'] = 'O2:0.8, N2:0.2'  # oxidizer composition (right)
    configs['fuel_composition'] = 'H2:0.9, N2: 0.1'  # fuel composition (left)
    
    # Reaction mechanism
    configs['mechanism'] = 'mechanism/burke/burke.yaml'
    
    # Two-Point Controls
    configs['spacing'] = 0.95
    configs['temperature_increment'] = 10
    configs['maximum_temperature_increment'] = 25
    configs['temperature_increment_growth_rate'] = 1.03 #Give this as a factor where 1.02 would be a 2% increase per iteration
    
    #Callbacks if they need to be used
    configs['steady_callback'] = steady_callback
    configs['unsteady_callback'] = unsteady_callback
    #'''

    '''
    # Ray Speths' Test
    # Input parameters
    configs['p'] = ct.one_atm  # pressure [Pascals]
    configs['fuel_inlet_temp'] = 300.0  # fuel inlet temperature [Kelvin]
    configs['oxidizer_inlet_temp'] = 300.0  # oxidizer inlet temperature [Kelvin]
    configs['fuel_mdot'] = 0.2 # kg/m^2/s
    configs['width'] = 50e-3 # Distance between inlets [m]
    configs['loglevel'] = 1  # amount of diagnostic output (0 to 9)

    #Mass Fraction composition
    configs['oxidizer_composition'] = 'O2:0.21, AR:0.78'  # oxidizer composition (right)
    configs['fuel_composition'] = 'H2:1, AR:1'  # fuel composition (left)
    
    # Reaction mechanism
    configs['mechanism'] = 'h2o2.yaml'
    
    # Two-Point Controls
    configs['spacing'] = 0.95
    configs['temperature_increment'] = 10
    configs['maximum_temperature_increment'] = 25
    
    #Callbacks if they need to be used
    configs['steady_callback'] = steady_callback
    configs['unsteady_callback'] = unsteady_callback
    '''

    '''
    # Input parameters
    configs['p'] = 5.4e6  # pressure [Pascals] 
    configs['fuel_inlet_temp'] = 800.0  # fuel inlet temperature [Kelvin]
    configs['oxidizer_inlet_temp'] = 711.0  # oxidizer inlet temperature [Kelvin]
    configs['fuel_mdot'] = 0.13 # kg/m^2/s
    configs['width'] = 60e-3 # Distance between inlets [m]

    configs['oxidizer_composition'] = 'O2:0.9449, H2O:0.0551'  # oxidizer composition (right)
    configs['fuel_composition'] = 'H2:0.4018, H2O: 0.5982'  # fuel composition (left)

    configs['loglevel'] = 1  # amount of diagnostic output (0 to 5)

    configs['mechanism'] = 'mechanism/H2.yaml'
    configs['steady_callback'] = steady_callback
    '''

    generate_flamelets(configs)

[wandadars]

I'm getting back to this feature now to get these last few issue resolved.

[wandadars]

@ischoegl It is kind of confusing how internal residual functions set some default BCs and then boundary objects have to come by again and change/adjust the boundaries?

[ischoegl]

@ischoegl It is kind of confusing how internal residual functions set some default BCs and then boundary objects have to come by again and change/adjust the boundaries?

This has been Cantera’s approach for a long time. I don’t really see an obvious alternative, either?

[wandadars]

@ischoegl It is kind of confusing how internal residual functions set some default BCs and then boundary objects have to come by again and change/adjust the boundaries?

This has been Cantera’s approach for a long time. I don’t really see an obvious alternative, either?

Yeah. I've been trying to debug an issue that has popped up with something not being set properly when transitioning from a non-two-point control to the two-point control, which causes the steady solver to fail. I found myself bouncing back and forth between the StFlow and the Boundaries1D to double check that something that I thought had been set didn't get switched by the other class. I don't have any more constructive thoughts on it other than that it seems like the boundary values would be the responsibility of the boundary objects.

Do you/Ray/Bryan have any go-to approaches for debugging a failing steady-state solve? I've set the log level to 8 to get outputs for the solution and residuals to the YAML debug file. There seemed to be some sort of bug where if log level was 8, the failing newton solve section would throw an error because of the part in the Sim1D where it tries to write debug information before and after taking a timestep, but the prescence of the first "debug" section in the YAML file cases the writeHeader to complain that the group already exists when the "after timestepping" part tries to write a "debug" entry to the file. I ended up passing the overwrite flag as true for those writes. I'm not sure if that isn't the correct thing to do there.

I also have a steady-state callback in the python script that I'm using to plot all the variables of interest during each steady-state solve attempt. I'm just running a bland H2/O2 burke mechanism with 1200K boundaries at 8e5 Pa, ideal gas. The initial steady-state solve works fine, but I'm working on trying to figure out when I turn on the two point control and attempt a single steady-state solve, the solver fails.

[wandadars]

I've updated my debugging python script if anyone wants to reproduce the crash scenario.

[wandadars]

Just so I don't forget, I ran across a strange execution where two calls were made to evalLambda, one with jmin=0 and jmax=41 (for a 42 point grid), and then a call with jmin=0 and jmax=1. This second call went into the interior loop at the bottom of the eval method for exactly one iteration because of the coding that sets j0 = max(jmin, 1) and that the for-loop is inclusive of the upper bound. It may not be an actual issue, but I was surprised by it when I saw some output that I wasn't expecting with some print statements I had placed.

Not sure if an update that comes to evalLambda for only once cell, which ends up updating the state of more than the requested cell is a non-desirable effect.

void StFlow::evalLambda(double* x, double* rsd, int* diag,
                        double rdt, size_t jmin, size_t jmax)
{
    if (!m_usesLambda) { // disable this equation
        for (size_t j = jmin; j <= jmax; j++) {
            rsd[index(c_offset_L, j)] = lambda(x, j);
            diag[index(c_offset_L, j)] = 0;
        }
        return;
    }

    if (jmin == 0) { // left boundary
        if (m_twoPointControl) {
            rsd[index(c_offset_L, jmin)] = lambda(x,jmin+1) - lambda(x,jmin);
        } else {
            rsd[index(c_offset_L, jmin)] = -rho_u(x, jmin);
        }
    }

    if (jmax == m_points - 1) { // right boundary
        rsd[index(c_offset_L, jmax)] = lambda(x, jmax) - lambda(x, jmax-1);
        diag[index(c_offset_L, jmax)] = 0;
    }

    // j0 and j1 are constrained to only interior points
    size_t j0 = std::max<size_t>(jmin, 1);
    size_t j1 = std::min(jmax, m_points - 2);
    double epsilon = 1e-8; // Precision threshold for being 'equal' to a coordinate
    bool conditionMet = false;
    writelog("Coming into evalLambda: j0: {:d}, j1: {:d}\n",jmin, jmax);
    for (size_t j = j0; j <= j1; j++) { // interior points
        if (m_twoPointControl) {
            //Print out the first if statement conditions for debugging and include the difference between grid(j) and m_zLeft
            writelog("grid(j) = {:10.4e}, m_zLeft = {:10.4e}, diff = {:10.4e}\n", grid(j), m_zLeft, std::abs(grid(j) - m_zLeft));
            if (std::abs(grid(j) - m_zLeft) < epsilon ) {
                writelog("We made it boys!\n");
                rsd[index(c_offset_L, j)] = T(x,j) - m_tLeft;
                conditionMet = true;
            } else if (grid(j) > m_zLeft) {
                rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
            } else if (grid(j) < m_zLeft) {
                rsd[index(c_offset_L, j)] = lambda(x,j+1) - lambda(x,j);
            }
        } else {
            rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
        }
    }
    //print the booleans using writelog for debugging, using the booleans precision format

    writelog("m_twoPointControl = {:d}, conditionMet = {:d}\n", m_twoPointControl, conditionMet);
    if (m_twoPointControl && !conditionMet) {
        throw CanteraError("StFlow::evalLambda",
            "The left boundary condition for the lambda equation was not met.");
    }
}

[ischoegl]

Yeah. I've been trying to debug an issue that has popped up with something not being set properly when transitioning from a non-two-point control to the two-point control, which causes the steady solver to fail. I found myself bouncing back and forth between the StFlow and the Boundaries1D to double check that something that I thought had been set didn't get switched by the other class. I don't have any more constructive thoughts on it other than that it seems like the boundary values would be the responsibility of the boundary objects.

I can see that, but don't necessarily object to the side effect that StFlow in effect handles most of the cases correctly, which avoids redundant implementations in Boundaries1D. It's definitely a historical choice, but I personally can live with that (and still believe that it is the cleaner implementation as there can be StFlow specializations that should not be propagated into Boundaries1D).

Do you/Ray/Bryan have any go-to approaches for debugging a failing steady-state solve?

I don't, but @speth may. There is some write-up on debugging here, but I have not really done much with that. Instead, I typically create a breakout from the Python API, which was largely sufficient for what I have done thus far.

I've set the log level to 8 to get outputs for the solution and residuals to the YAML debug file. There seemed to be some sort of bug where if log level was 8, the failing newton solve section would throw an error because of the part in the Sim1D where it tries to write debug information before and after taking a timestep, but the prescence of the first "debug" section in the YAML file cases the writeHeader to complain that the group already exists when the "after timestepping" part tries to write a "debug" entry to the file. I ended up passing the overwrite flag as true for those writes. I'm not sure if that isn't the correct thing to do there.

Interesting. I was not aware that there are write steps before and after taking a timestep (and am also wondering whether this saves redundant info in subsequent time steps?). In any case, your workaround looks reasonable; as this presumably is also an issue for the main branch and not only this PR, please feel free to file an issue report.

[wandadars]

@speth I was wondering if you could help me check my understanding of something in the Sim1D.cpp file. What is the logic for what should be stored in the debug_sim1d.yaml file during an unsuccessful steady-state solve?

I was thinking that the save() and saveResidual() calls in the solve() method that happen after a failed steady-state solve would write the state of the system to the yaml file, and then there is the section immediately after where the code tries to take some timesteps, at which time it also prints information to the debug yaml file. I figure that this data would be stored in the yaml file as something like a "state of the system after we take N timesteps" entry in the yaml file. Is the timestep coding trying to take N timesteps from the failed solution state? Or does it take N timesteps from the initial state before the newton solve was attempted?

Is it like this:
(Initial State) --->[NewtonSolve]--->(Failed State)--->[Timestepper]--->(New Initial State)---> |
^---------------------------------------------------------------------------------------------------------------<<|

Or this?
(Initial State)--->[NewtonSolve]--->(Failed State) : (Initial State)--->[Timestepper]--->(New Initial State)---> |
^------------------------------------------------------------------------------------------------------------------------------------<<|

I'm just wondering about whether it is worth it to write an entry in the debug yaml file for the state before the attempted Newton solve so that we can capture the before/after state around the state-state newton solve attempt. I've written a simple Bokeh server visualizer to visualize the full state of the system for each entry in the yaml file. I added new high-level groups for "debug_afterTimestepping" and "residual_afterTimestepping" to see changes. I recently added a "debug_beforeNewton" and "residual_beforeNewton" group by calling the save/saveResidual methods before the call to the newtonSolve() method. I've noticed that when I plot all three, the timestepping result and the beforeNewton seem to be the same, which I suppose means that the first option above is probably what is happening.

Not sure if I'm going about this all wrong, but I was trying to get some sort of more detailed state information that was easy to parse for when I'm adjusting BCs and equations where I can easily see the changes happening in the domain.

[speth]

Is the timestep coding trying to take N timesteps from the failed solution state? Or does it take N timesteps from the initial state before the newton solve was attempted?

I believe that the timestepping starts from the state before the attempted Newton solve. You can see this in Sim1D::solve because the call to do the timestepping is:

dt = timeStep(nsteps, dt, m_state->data(), m_xnew.data(), loglevel-1);

where m_state is the initial condition for the time stepping, and in Sim1D::newtonSolve, m_state is updated only if the Newton solve succeeds.

I'm just wondering about whether it is worth it to write an entry in the debug yaml file for the state before the attempted Newton solve so that we can capture the before/after state around the state-state newton solve attempt.

I think the state written to the debug file and labelled as "After unsuccessful Newton solve" is the state before the attempted Newton solve, since m_state is not updated in this case. You might try swapping m_state and m_xnew, calling save, and then swapping back if you want to get whatever state the Newton solver exited with. I'm not sure what condition that state vector is in when it fails, though. This is a bit like what saveResidual doe to work around the fact that save only uses the contents of m_state.

[wandadars]

I was looking at my "setRightControlPoint() and setLeftControlPoint() methods for the two-point control method, and there's already a much better implementation in the setFixedTemperature() method used for freely propagating flames. I was trying to see if I could find a way to extricate that conditional check in that method that looks at the return value of isFree() and find a way to generalize for any domains that want to inert a specified temperature point in a domain.

The coding in that method has me a bit confused because it is looping over several domains to insert the new point. Is that because the domains use a global numbering and shoving a single point into a domain would mess up the other (n-1) domains around the free-flame domain? Also the extra work the method does of saving the grid and solutions for all domains that are not the freely propagating domain seems like it shouldn't necessarily be required, but I can see that the final section that loops over the domains uses some sort of daisy-chaining of indices to update grids for all domains.

I was thinking that being able to call that setFixedTemperature() twice with the left and right control point temperatures would be nice, but it seemed a bit ham-fisted to just add a third "|| d_free->twoPointControlEnabled()" to the if-statement.

[speth]

Yes, I think it would make sense to adapt the setFixedTemperature method to handle all cases. However, just adding || d_free->twoPointControlEnabled() to the condition under which the fixed point is set is not sufficient. The setFixedTemperature method specifically sets the StFlow::m_zfixed and m_tfixed member variables, which are then used in evaluating the StFlow governing equations. What you need to be able to do here is change where the information about the z and T of the fixed point are stored. I guess you also need to control the order in which the domains / points are scanned, since to set the right control point, you want the rightmost point where the temperature crosses the specified value.

The coding in that method has me a bit confused because it is looping over several domains to insert the new point. Is that because the domains use a global numbering and shoving a single point into a domain would mess up the other (n-1) domains around the free-flame domain?

The idea is that the simulation could consist of multiple flow domains (set up as boundary-flow-boundary-flow-boundary, for example). If that were a freely propagating flame, you would still only set one temperature fixed point, and the location where the temperature reaches the specified value could be in either domain, so you have to loop over all of them.

Also the extra work the method does of saving the grid and solutions for all domains that are not the freely propagating domain seems like it shouldn't necessarily be required, but I can see that the final section that loops over the domains uses some sort of daisy-chaining of indices to update grids for all domains.

The grid and the state vector are both global, covering all domains. So to insert a point, you need to effectively push back the data for all the remaining points in both the grid and state vector.

[wandadars]

An update, I have been able to run a script that loops over all the burning flame configurations (by symmetrically increasing the left and right mdot values) until flame extinction, then reloading the last solution and enabling two-point control to march down the unstable branch. It does seem capable of marching down the unstable branch. One thing that must be accounted for is that the drop in the maximum temperature of the flame is very sensitive to the temperature increment at the control points during the unstable branch walk-down. So in the script that I'm using, I often will be saving previous solutions and restoring them if failures occur due to the 1D solver being unable to solve for a case where the flame changes significantly from the initial condition.

I haven't worked out a particularly satisfying way to re-use that older fixed temperature coding yet. I'm going to look at it again today, but if I can't make progress on that, I don't think it should hold up the rest of this PR.

[wandadars]

I think we can move forward with pushing this PR across the finish line. I've been running some cases to compare the results for the two-point control method versus FlameMaster, and it looks like the feature is able to perform similarly to the arclength continuation method that is implemented in FlameMaster.

FlameMaster S-curve for a Allstar H2/O2 mechanism.

Cantera two-point control S-curve for the same Allstar mechanism.

I couldn't figure out a way to extricate the coding from that setFixedTemperature() method that is used for free-flames. I think for the time being, having the methods that are particular to the two-point control is what we'll do.

Some things that I have learned so far is that for mechanisms, the initial placement of the control points is important. From the 1996 Nishioka paper, they say that the points should be placed in regions of steep temperature gradients, which I had interpreted as being far from the peak location. But from using the feature, it looks like for some cases, having control points close to the peak yields a more stable procession of incremental solutions.

I heavily utilize the save() and restore() methods during a process of walking down the unstable branch. If a particular solution fails to converge, I restore a previous solution and decrease the temperature increment and try again until either success is achieved or a minimum temperature increment threshold is passed at which time the script bails out with a failure message.

When the two-point control is activated, the strain rates are quite high and the flame is incredibly thin compared to the domain width, so the refinement is a critical aspect in order to keep the flame resolved.

The only caveat for people using the feature would be that if you have a solution (call is solutionA) and then attempt to adjust the temperature control points and decrement the temperature at the control points, if that process fails, then you'll have to restore solutionA. When solutionA is restored, the control points will have to re-located on this previous grid. I personally perform the restore, then set the control points for the restored solution, then call the solve() method just to make sure that the flame is in a state where the control points and other variables are all in agreement with each other (Uo and Lambda being consistent with the control points). The exact interplay between the location of the control points, the residual vectors, and the refinement feature is still not perfectly clear to me.

[wandadars]

I'm working on adding a unit test that runs an H2/O2 counterflow diffusion flame, gets the temperature profile, then activates the two-point control, sets control points based on the maximum flame temperature, and the re-solves the governing equations with the two-point control activated (but no temperature decrement at the control points). This should give a solution that is similar to the solution without the two-point control method.

I'm having some issues though with the unit tests hanging when I try to run them. I haven't figured it out yet, but it might have something to do with the output files having two-point control entries added to them and that causing an issue when tests are comparing with files that don't have those entries.

This is the test for reference:

def test_two_point_control(self):

        p = ct.one_atm
        fuel='H2:1.0, AR:1.0'
        T_fuel=300
        mdot_fuel=0.05
        oxidizer='O2:0.21, AR:0.78'
        T_ox=300
        mdot_ox=0.1
        width=0.05

        # Solution object used to compute mixture properties
        gas = ct.Solution("h2o2.yaml")
        gas.TP = T_fuel, p

        # Flame object
        sim = ct.CounterflowDiffusionFlame(gas, width=width)

        # Set properties of the fuel and oxidizer mixtures
        sim.fuel_inlet.mdot = mdot_fuel
        sim.fuel_inlet.X = fuel
        sim.fuel_inlet.T = T_fuel

        sim.oxidizer_inlet.mdot = mdot_ox
        sim.oxidizer_inlet.X = oxidizer
        sim.oxidizer_inlet.T = T_ox

        sim.solve(loglevel=2)
        temperature_1 = sim.T

        sim.two_point_control_enabled = True
        spacing = 0.95
        control_temperature = np.min(sim.T) + spacing * (np.max(sim.T) - np.min(sim.T))
        sim.set_left_control_point(control_temperature)
        sim.set_right_control_point(control_temperature)
        print('HELLO')
        sim.solve(loglevel=2)
        print('HELLO2')
        temperature_2 = sim.T

        # Check difference between the two solutions
        self.assertArrayNear(temperature_1, temperature_2, rtol=1e-4, atol=1e-6)

[wandadars]

I got the unit test for the two-point control to work. Still seeing some other tests fail or hang for some reason.

@speth I have a small description of the two-point method in the include file for the StFlow.h. Do I need to expand on this method in more detail somewhere else in the Cantera documentation? (Like an .rst file or something?).

[ischoegl]

@wandadars ... thank you for your continued work on this. While I won't be able to add much in terms of detailed review at the moment, I wanted to answer your question about documentation: yes, it would be great if you could add some information on the theory in the doc/sphinx/reference/onedim folder. Also, adding an example in samples/python/onedim would be terrific. At this point, you can add BibTeX-style citations in both doxygen and Sphinx documentations.

On a separate note, there appears to be a critical/inadvertent change in computational behavior, as at this point a lot of the unit test suite appears to be timing out due to the 1 hour time limit (rather than running within a few minutes).

[speth]

I pushed a fix for the test that was causing the CI to hang. The test in question was getting stuck regridding a flame over and over, where it would add and then remove the same grid points on successive iterations. I think we were just getting lucky that this didn't happen before, and was only triggered by the addition of the equation for Uo in this PR. While this is essentially a no-op when the two point control is disabled, it does affect the solver a little bit.

[wandadars]

@speth I get a failing test when trying to round-trip the parameters that get saved in the yaml file.

filename = self.test_work_path / "two_point_control.yaml"
    
        original_settings = sim.flame.settings
        sim.save(filename)
        sim.restore(filename)
        filename.unlink()
    
        restored_settings = sim.flame.settings
        assert 'continuation-method' in restored_settings
        assert restored_settings['continuation-method']['type'] == original_settings['continuation-method']['type']
>       assert restored_settings['continuation-method']['left-location'] == original_settings['continuation-method']['left-location']
**E       assert 0.0225 == 0.022500000000000003**

This is a very small difference and I was going to switch the test to an approx, but I wanted to double check with you: Should we expect perfect equality here in the case of reading data from the YAML file?

For reference, this is what gets written:

continuation-method:
      type: two-point
      left-location: 0.0225
      right-location: 0.0275
      left-temperature: 2379.341155103330
      right-temperature: 2516.019408460186

[wandadars]

I'm getting core dump errors for unit tests that deal with the two_point_control_enabled property. Should all the python properties for the two-point control method be in the FlameBase class in onedim.py, or should those be down in the CounterflowDiffusionFlame class?

I haven't quite grasped it yet, but the following coding gives the core dumps when the commented section is uncommented.

//! Sets the status of the two-point control
    void enableTwoPointControl(bool twoPointControl) {
        m_twoPointControl = twoPointControl;
        /*
        if (m_usesLambda){
            m_twoPointControl = twoPointControl;
        } else {
            throw CanteraError("StFlow::enableTwoPointControl",
                "Invalid operation: two-point control can only be used with axisymmetric flames.");
        }
        */
    }

The m_usesLambda variable is protected, but this method is part of the same class, so I didn't think that would be an issue, but it appears to be.

Edit:

While the m_usesLambda is a protected variable, the method in question is also a part of the same class, which I thought meant it could access protected variables. I chained the direct use the variable to the isStrained() method, which returns the value of the m_usesLambda variable, and that stopped the weird core-dumps that I was seeing when calling the method from the Python interface.

Edit 2: I spoke too soon. Now I'm still getting those core-dumps. Before it was in the coding that get/set the two_point_control_enabled property, and changing that direct access of m_usesLambda to a function-based access worked. But now, I get something strange where, earlier in the tests I call sim.right_control_point_temperature and there aren't any issues, but if I set the two_point_control_enabled to false and then try to call that right_control_point_temperature property, the pytest reports a core dump.

[wandadars]

I may not understand how the c++ exceptions make it up to the Python level. That might be what is causing some of the issues that I'm having.

[speth]

I may not understand how the c++ exceptions make it up to the Python level. That might be what is causing some of the issues that I'm having.

For any C++ method that might raise an exception, you need to add except +translate_exception to the declaration in the .pxd file so Cython knows how to convert it to a Python exception.

Pytest's output capturing capabilities unfortunately suppress the useful error message you would get otherwise. If you run the test outside of SCons (may require some fiddling with paths; see) with a command something like:

pytest -raP -s --verbose test/python/test_onedim.py -k two_point

where the -s flag disables output capturing, you'll see something like the following in the output:

terminating due to uncaught exception of type Cantera::CanteraError: 
*******************************************************************************
CanteraError thrown by StFlow::setRightControlPointTemperature:
Invalid operation: two-point control is not enabled.
*******************************************************************************

[speth]

Should we expect perfect equality here in the case of reading data from the YAML file?

I think roundoff-level errors have to be tolerated when comparing values that have been converted to string.

[wandadars]

I may not understand how the c++ exceptions make it up to the Python level. That might be what is causing some of the issues that I'm having.

For any C++ method that might raise an exception, you need to add except +translate_exception to the declaration in the .pxd file so Cython knows how to convert it to a Python exception.

Pytest's output capturing capabilities unfortunately suppress the useful error message you would get otherwise. If you run the test outside of SCons (may require some fiddling with paths; see) with a command something like:

pytest -raP -s --verbose test/python/test_onedim.py -k two_point

where the -s flag disables output capturing, you'll see something like the following in the output:

terminating due to uncaught exception of type Cantera::CanteraError: 
*******************************************************************************
CanteraError thrown by StFlow::setRightControlPointTemperature:
Invalid operation: two-point control is not enabled.
*******************************************************************************

Ah, that's it. Thanks for clearing that up. Is there a way to view what the documentation looks like with the edits to the markdown files?

[speth]

Ah, that's it. Thanks for clearing that up. Is there a way to view what the documentation looks like with the edits to the markdown files?

Yes, though it's a little tricky to find. At the top of this page, Go to the "Checks" tab, the select "CI" from the list on the left. Scroll all the way down on the right until you get to the "Artifacts" section. The "docs" artifact contains the documentation build using your PR.

[wandadars]

@speth I believe I have addressed all the review comments/suggestions that you provided.

[speth]

Thanks, @wandadars, I think this is looking pretty good. I have a few more comments, mainly regarding the documentation, the example, and some minor formatting quibbles. This should be good to go after these are addressed.

[speth]

This looks good to me. I just added a couple of minor formatting changes. Once this passes CI, I'll merge it. Thanks for bringing the long saga of adding this feature to Cantera to a successful end, @wandadars.

[ischoegl]

Likewise thanks!