WIP: 1d point control

[wandadars]

This pull request is based on the pull request #766.

Changes proposed in this pull request

• The goal is to add a type of continuation method that allows for 1D counter-flow diffusion flame solutions that exist along the unstable burning branch.
• Minor updates to the original pull request
• An attempt to address issues that were brought up in the original pull request

The original work appears to be based on the method described in the following paper:
https://www.yang.gatech.edu/publications/Journal/C&F%20(2014,%20Huo)%20-%20S%20curve.pdf

Reference information: https://www.sciencedirect.com/science/article/abs/pii/S0010218014001825

The original comments by @speth regarding work that needed to be done with this feature center around numerical stability when using the feature.

I would also like to understand the significance of the Strain Equation option that was present in the original work (#541), and why it was not included in this update of #766. To understand if this has anything to do with the stability issues that are being seen when using the feature in its current form.

Another relevant reference that discusses the governing equations used in the previous reference in more detail is: https://www.sciencedirect.com/science/article/abs/pii/S0010218008001569 so that I can compare with what Cantera uses (https://cantera.org/science/flames.html)

[wandadars]

Is there a reference for just the general coordinate/layout of the stagnation flow problem in the StFlow class? I know it's based on the stagnation flow sections in Kee's book. I was just wondering if there was anything regarding the coordinate system, where the origin is located, etc.

[speth]

I'm not quite sure what you mean. I assume you've seen https://cantera.org/science/flames.html? The location of the origin (z=0) is arbitrary, though as a practical matter the first grid point is set as z=0 and positive z values for subsequent grid points in all of the pre-configured flame types in Cantera.

[wandadars]

Thanks Ray. I reviewing the 1D numerical coding in the StFlow.cpp file with regards to what the 'left' and 'right' boundaries were in the context of the axisymmetric stagnation flow configuration. Things like the unwinding scheme that Kee describes in his book, and how the Lambda equation(c_offset_L) is treated. To make sure I understand, the fuel bc is z=width and the oxidizer bc is z=0 for the counterflow flame configuration?

[wandadars]

There was an issue in the Sim1D.cpp implementation of the method for finding the coordinate of the control point based on temperature, as well as the coding in StFlow.cpp implementation for the residuals at the control points. The coding swapped from @axbriones approach of identifying the grid index of the point in the domain that was to be used for controlling the flame, and so the loop over the grid indices checked for integer equality. The newer work changed to a coordinate-based identification of the flame control point location, using floating point numbers, but the other parts of the code that used to check for grid index equality between integers just had the integers swapped out for the floating point numbers. Floating point comparisons are definitely not the way to go. I added an epsilon for comparing the equality of two floating point numbers. I'm not sure why the original grid index based way of locating the control point was moved away from. Possibly related to grid refinement changing the number of grid points?

[wandadars]

A gif of the temperature field during the two-point adjustment.

[wandadars]

Having some issues with convergence when the temperature gets to around 1500K. Not sure yet what is causing the issue yet.

Below is the script that I'm currently using to try and march through the flame solutions using the two-point control.

import os
import time
from collections import defaultdict

import matplotlib.pyplot as plt
import numpy as np

import cantera as ct


def create_folder(path):
    try:
        if os.path.isdir(path):
            print("Error: The directory you're attempting to create already exists") # or just pass
        else:
            os.makedirs(path)
    except IOError as exception:
        raise IOError('%s: %s' % (path, exception.strerror))
        return None


def plot_solution(f, solution_num):
    #Plot temperature
    fig, ax = plt.subplots()
    ax.plot(f.grid, f.T, 'o-')
    ax.set_ylim([300, 2800])
    plt.savefig('output/temperature/flamelet_' + str(solution_num+1) + '.png')
    plt.close()
    
    #Plot Axial oxidizer velocity
    fig, ax = plt.subplots()
    ax.plot(f.grid, f.Uo, 'o-')
    plt.savefig('output/Uo/flamelet_' + str(solution_num+1) + '.png')
    plt.close()
    
    #Plot Axial oxidizer velocity
    fig, ax = plt.subplots()
    ax.plot(f.grid, f.L, 'o-')
    plt.savefig('output/L/flamelet_' + str(solution_num+1) + '.png')
    plt.close()
    
    #Plot Axial velocity
    fig, ax = plt.subplots()
    ax.plot(f.grid, f.velocity, 'o-')
    plt.savefig('output/u/flamelet_' + str(solution_num+1) + '.png')
    plt.close()
    
    #Plot spread rate
    fig, ax = plt.subplots()
    ax.plot(f.grid, f.spread_rate, 'o-')
    plt.savefig('output/V/flamelet_' + str(solution_num+1) + '.png')
    plt.close()
 
    fig = plt.figure(figsize=(10,10))
    axs = []
    axs.append(plt.subplot2grid((3,2), (0,0))) # top left
    axs.append(plt.subplot2grid((3, 2), (1, 0)))  # middle left
    axs.append(plt.subplot2grid((3, 2), (2, 0)))  # bottom left
    axs.append(plt.subplot2grid((3, 2), (0, 1)))  # top right
    axs.append(plt.subplot2grid((3, 2), (1, 1)))  # middle right
    
    axs[0].plot(f.grid, f.T, 'o-')
    axs[0].set_xlabel('Distance, m')
    axs[0].set_ylabel('Temperature, K')
    axs[0].set_ylim([200, 2900])
    
    axs[1].plot(f.grid, f.velocity, 'o-')
    axs[1].set_xlabel('Distance, m')
    axs[1].set_ylabel('Axial Velocity, m/s')
    
    axs[2].plot(f.grid, f.spread_rate, 'o-')
    axs[2].set_xlabel('Distance, m')
    axs[2].set_ylabel('Spread Rate, 1/s')
    
    axs[3].plot(f.grid, f.Uo, 'o-')
    axs[3].set_xlabel('Distance, m')
    axs[3].set_ylabel('Oxidizer Velocity, m/s')
    
    axs[4].plot(f.grid, f.L, 'o-')
    axs[4].set_xlabel('Distance, m')
    axs[4].set_ylabel('Lambda, m/s')
    fig.tight_layout()
    plt.savefig('output/combined/flamelet_' + str(solution_num+1) + '.png')
    plt.close()


def create_output_folders():
    create_folder('output/temperature')
    create_folder('output/Uo')
    create_folder('output/L')
    create_folder('output/u')
    create_folder('output/V')
    create_folder('output/combined')


# Input parameters
p = ct.one_atm  # pressure
tin_f = 300.0  # fuel inlet temperature
tin_o = 300.0  # oxidizer inlet temperature
mdot_f = 0.2  # kg/m^2/s
width = 0.03 # Distance between inlets

comp_o = 'O2:0.21, AR:0.78'  # oxidizer composition
comp_f = 'H2:1.0'  # fuel composition

loglevel = 1  # amount of diagnostic output (0 to 5)

gas = ct.Solution('h2o2.yaml')
gas.TPX = tin_f, p, comp_f
density_f = gas.density

gas.TPX = tin_o, p, comp_o
density_o = gas.density

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

# Set the state of the two inlets
f.fuel_inlet.mdot = mdot_f
f.fuel_inlet.X = comp_f
f.fuel_inlet.T = tin_f

f.oxidizer_inlet.mdot = (mdot_f / density_f) * density_o * 0.5  
f.oxidizer_inlet.X = comp_o
f.oxidizer_inlet.T = tin_o

#f.set_refine_criteria(ratio=4, slope=0.2, curve=0.3, prune=0.04)
f.set_refine_criteria(ratio=10, slope=0.1, curve=0.2, prune=0.04)

# Solve the problem
f.solve(loglevel, auto=True)
f.save('previous_solution.yaml',overwrite=True)

Tmax = max(f.T)
previous_Tmax = Tmax
print(f"The maximum temperature is {Tmax:.4g} K")

Q = [np.trapz(f.heat_release_rate, f.grid)]
T_max = [Tmax]
a_mean = [f.strain_rate('mean')]

#Point control
f.twoPointControl_enabled = True
spacing= 0.90
Tmax_threshold = 15

f.flame.set_transient_tolerances(default=(1e-6, 1e-12)) 
f.flame.set_steady_tolerances(default=(1e-6, 1e-12)) 
original_temperature_increment = 10 #Kelvin
failed_previous_attempt = False
temperature_increment = original_temperature_increment  #Starting value
flamelet_solution_num = 1
iterations = 1
create_output_folders()
while np.max(f.T) > 300 and iterations < 1000:
    print('Flame: Tmin: ' + str(np.min(f.T)) + ', Tmax: ' + str(np.max(f.T)))
    new_temp = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T)) 
    print('Desired Internal Boundary Temperature: ' + str(new_temp))
    f.set_fuel_side_boundary(new_temp) #Sets a location based on the temperature of the current flame
    f.set_oxid_side_boundary(new_temp)
    print('Cloest found values: Fuel Temp: ' + str(f.tFuel) + '   Oxidizer Temp: ' + str(f.tOxid))
    #Lowers the control points temperature
    f.tFuel -= temperature_increment
    f.tOxid -= temperature_increment
    print(f'Setting temperature points to T = {f.tFuel}')

    # Use callback function to collect data after each timestep
    # This data will be available after a failed solution attempt
    grids = {}
    TT = defaultdict(list)
    YY = defaultdict(list)
    def callback(x):
        N = len(f.grid)
        grids[N] = f.grid
        TT[N].append(f.T)
        YY[N].append(f.Y)
        return 0
    callback(0)
    f.set_time_step_callback(callback)

    iterations +=1
    try:
        f.solve(loglevel)
    except Exception as err:
        print(f'error at T = {f.tFuel}')
        print(err)
        print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
        #Update increment
        temperature_increment *= 0.5
        if(temperature_increment < 0.1):
            spacing = spacing*1.01 # Reduce spacing to closer to flame to see if that helps
            print('No more solutions reached at temperature increment: ' + str(temperature_increment) + ', adjusting spacing to: ' + str(spacing))
            temperature_increment = original_temperature_increment # Reset the increment
            if spacing > 0.999:
                print('No more solutions reached. Stopping at temperature increment: ' + str(temperature_increment) + ', and spacing: ' + str(spacing))
                break
        failed_previous_attempt = True
        #f.restore('previous_solution.yaml')
        continue
        
#    if abs(np.max(f.T) - previous_Tmax) > Tmax_threshold:
#        print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
#        #Undo previous failed increment to internal bcs
#        f.tFuel += temperature_increment
#        f.tOxid = f.tFuel
#        #Update increment
#        temperature_increment *= 0.5
#        if(temperature_increment < 0.1):
#            print('No more solutions reached at temperature increment: ' + str(temperature_increment) + ', adjusting spacing to: ' + str(spacing))
#            temperature_increment = original_temperature_increment # Reset the increment
#            break
#        failed_previous_attempt = True
#        f.restore('previous_solution.yaml')
#        continue


    failed_previous_attempt = False
    
    #f.save('previous_solution.yaml', overwrite=True) #Save a backup in case a solution fails
    Tmax = np.max(f.T)
    previous_Tmax = Tmax
    q = np.trapz(f.heat_release_rate, f.grid)
    a = f.strain_rate('mean')
    print(f"The maximum temperature is {Tmax:.4g} K (a = {a:.4g}, Q = {q:.4g})")
    T_max.append(Tmax)
    a_mean.append(a)
    Q.append(q)
    
    flamelet_solution_num += 1
    plot_solution(f, flamelet_solution_num)


#Plot process
fig, ax = plt.subplots()
ax.plot(a_mean, T_max, 'o-')
plt.savefig('tMax_vs_strainRate.png')

fig, ax = plt.subplots()
ax.plot(a_mean, Q, 'o-')
plt.savefig('Q_vs_strainRate.png')

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

Here's a plot of the history of the main variables. I'm seeing a Jacobian is singular for domain flame, component Uo at point 173 error prior to the failure to converge. That discontinuity in the spread_rate solution doesn't look good.

[wandadars]

I observed that the ability to get solutions down near the 1000K temperature seems to be dependent on the grid not refining. In this approach the interior control points are located using that set_fuel_side_boundary() method. I'm wondering if during a refinement, we would need to re-call those methods after points have been deleted/added to the domain.

Some notes I made about the issue:

[wandadars]

Having the grid refinement disabled results in a more stable march down the unstable burning branch. Ideally, in the solve() method we would need to re-evaluate the location of the control points at the start of each refinement loop. I see that there's a steady callback that can be called, which is where I would put the calls to method that set the internal boundary locations, but I don't quite understand if that would be something I could do. I don't have any experience with callbacks.

[wandadars]

Have I been rebasing correctly on this PR @speth? It's not something I've done very often. I see that @ischoegl has been doing some work in the oneD part of the code, and if I have been rebasing incorrectly, I may open a fresh PR on top of the newest version.

[speth]

It doesn't look like you have pushed anything to this branch recently -- the history here says the last time you pushed was July 6. Have you rebased more recently than that?

[wandadars]

I have not. I had noticed that Mr. Schoegl had done some work recently and I was thinking about trying to rebase my changes again against his. I'm just not sure if my previous rebasing is actually being done correctly, so if I'm messing up the commit history in some way, I wanted to see if you could take a look to double-check.

[ischoegl]

Some of the (presumably extensive) merge conflicts are due to this PR introducing major code reorganizations in StFlow, which involves a series of new methods, i.e.

void StFlow::evalLeftBoundary(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalInterior(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)

void StFlow::evalRadialMomentumEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalContinuityEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalSpeciesEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalEnergyEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalLEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalUoEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

I can definitely see why you're doing this, but relatively modest changes in StFlow that get merged into upstream main will unfortunately spiral out of control. I've done similar - albeit less extensive - things to split the very same (quite unwieldy) function in some of my own external work, and changed the way I approach things since.

That said, I don't believe that hard-coding individual governing equations into separate functions is the way to go any longer.

[wandadars]

Some of the (presumably extensive) merge conflicts are due to this PR introducing a major code reorganizations in StFlow, which involves a series of new methods, i.e.

void StFlow::evalLeftBoundary(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalInterior(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)

void StFlow::evalRadialMomentumEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalContinuityEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalSpeciesEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalEnergyEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalLEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

void StFlow::evalUoEquation(double* x, double* rsd, int* diag, double rdt, size_t j)

I can definitely see why you're doing this, but relatively modest changes in StFlow that get merged into upstream main will unfortunately spiral out of control. I've done similar - albeit less extensive - things to split the very same (quite unwieldy) function in some of my own external work, and changed the way I approach things since.

That said, I don't believe that hard-coding individual governing equations into separate functions is the way to go any longer.

Yes. I realized that I polluted the continuation work with residual equation refactoring instead of having them separate. I do have a better grasp on what the original PR was doing, so I believe I can start again with what's in the current branch and make only changes to add that feature.

Having the single method to evaluate every equation was a bit overwhelming for me. I was also confused as to why there was an evalRightBoundary method but no corresponding evalLeftBoundary method in StFlow.cpp.

[ischoegl]

Having the single method to evaluate every equation was a bit overwhelming for me. I was also confused as to why there was an evalRightBoundary method but no corresponding evalLeftBoundary method in StFlow.cpp.

Agreed. It is overwhelming! My hunch is that this bit of code has survived for this long with relatively small changes due to being overly complex. I hope that this will change eventually ... you may have noticed that I cleaned up the - rather confusing - boundary conditions (apart from locating several hardcoded indices), in a series of incremental steps. On a separate note, we replaced all doublereal instances recently as we're getting close to releasing Cantera 3.0.

[yeanment]

Having the grid refinement disabled results in a more stable march down the unstable burning branch. Ideally, in the solve() method we would need to re-evaluate the location of the control points at the start of each refinement loop. I see that there's a steady callback that can be called, which is where I would put the calls to method that set the internal boundary locations, but I don't quite understand if that would be something I could do. I don't have any experience with callbacks.

I am also interested in the flame control method in S-curve generation, and adapted the code based on your commit. Basically, I did one improvement: change the setFuelInternalBoundary and setOxidizerInternalBoundary function based on setFixedTemperature function, so that the control point has the exact temperature of provided Tfuel and Toxid, and this point will not be removed. Details of my implementation could be found in my scurve branch.

I am more interested in the counterflow premixed flame continuation, and conducted some test cases with my version. The script is based on the version provided above. However, I found the iterations are stucked at a later stage. The temperature does not decrease any more. I may also share a short movie on the temperature change throughout this process.

[wandadars]

Thanks @yeanment . I've been wanting to have a feature like this in Cantera for quite some time now, so any input from others who also share this goal is more than welcome. Getting a feature like this into Cantera in a robust form would be fantastic for expanding the usage of Cantera for flamelet modeling. From my understanding, the method first takes an existing solution, which has a temperature field T(xj), and then a user sets a temperature that they wish to use for locating the control point along the x-axis. We loop through and find the closest points to that desired point and store the location and the temperature. Then the user sets an updated temperature to use at those control points, and the algorithm them solves the governing equations with the two extra equations to drive the temperature at the original control point locations to the new temperature values.

In your change, you are taking the original starting temperature field T(xj), and then finding either a point that is sufficient close that it can be considered to the same, or you find the bounding temperature values around the desired value and interpolate to the exact value and insert an extra point into the domain at that location? Would that be a correct summary of your update to the algorithm? There's definitely some issue for when the flame gets strained strongly and becomes very thin, that the points along the flame profile become sparse when solving without grid refinement. I think the trick would be to find a way to make the grid refinement play nice with the control point locator methods. I believe this is the current issue with failure to converge at lower maximum temperatures.

[yeanment]

In your change, you are taking the original starting temperature field T(xj), and then finding either a point that is sufficient close that it can be considered to the same, or you find the bounding temperature values around the desired value and interpolate to the exact value and insert an extra point into the domain at that location? Would that be a correct summary of your update to the algorithm? There's definitely some issue for when the flame gets strained strongly and becomes very thin, that the points along the flame profile become sparse when solving without grid refinement. I think the trick would be to find a way to make the grid refinement play nice with the control point locator methods. I believe this is the current issue with failure to converge at lower maximum temperatures.

@wandadars I also want this feature for a long time, but does not have enough time to implement it till recently. Your summary is correct. Can you elaborate your point on " I think the trick would be to find a way to make the grid refinement play nice with the control point locator methods.? From my understanding, one fixed temperature will lead to one specific inlet velocity in the case of premixed counterflow twin flame. In my previous experience on s-curve continuation with arc-length continuation method, the grid does play an important role in getting the results convergent fast. If the arc-length continuation fails, I may change the computational domain or minimal grid size to make the computation continues. But in the flame-control method, I do not know how to make it work. Can you give me some suggestions on this?

Besides, I merge the main modification with the recent Cantera master branch. To ensure a minimal change to the original code, I only consider one-point continue, otherwise a lot of tests will fail due to the addition of a new equation. I am not sure if there are some better options in adding additional equations in Cantera without making test failed repeatedly. In my recent implementation, zero gradient is applied to the lambda equation except at the zFuel point (modified in StFlow.cpp). At this position, the fixed temperature boundary condition is imposed and the fixed inlet velocity boundary is thereby released for the left inlet point (modified in Boundary1D.cpp). You may check my recent commit on scurve branch. The code can run continuation. However, when I check the results, I found that the lambda value at the inlet is not the same with that at other points. Do you have any idea on this problem?

Below is some representative results. The first is the temperature during continuation. The second is Tmax versus inlet velocity. The third is one result during continuation. Almost everything is the same with your script except I plot the Qdot in Uo subplot since one-point continuation is considered for now.

[wandadars]

@yeanment From my understanding, if we have a set of boundary conditions, and we have an initial condition and a set of control points, the 1D solver sets up the system of equations for the standard governing equations and the Uo and Lamda equations. These last two equations are zero gradient everywhere except where the control points are located. What I meant by playing nice was that the control points are initially located prior to the solver being active. These points are located via a temperature, and in the algorithm here, those points are taken to be the points closest to that reference temperature. During the solution process, the solver may refine the grid and add several points that are closer to the reference temperature value. I'm not sure if this is an issue, but I was wondering if I could make a call to the functions that place the control points between each refinement, if that would improve the stability. As it is currently, if I do not disable grid refinement, the script fails to generate even a moderate number of solutions down the unstable branch, but with refinement disabled, I see that it can get down much closer to the extinction point.

I see that in your plots, the flame appears to be very close to the right inlet, and this might be causing issues. You may need to adjust your mass flow rates or domain width. Ideally you want that spike in the temperature field to be near the center of your domain.

[yeanment]

@yeanment From my understanding, if we have a set of boundary conditions, and we have an initial condition and a set of control points, the 1D solver sets up the system of equations for the standard governing equations and the Uo and Lamda equations. These last two equations are zero gradient everywhere except where the control points are located. What I meant by playing nice was that the control points are initially located prior to the solver being active. These points are located via a temperature, and in the algorithm here, those points are taken to be the points closest to that reference temperature. During the solution process, the solver may refine the grid and add several points that are closer to the reference temperature value. I'm not sure if this is an issue, but I was wondering if I could make a call to the functions that place the control points between each refinement, if that would improve the stability. As it is currently, if I do not disable grid refinement, the script fails to generate even a moderate number of solutions down the unstable branch, but with refinement disabled, I see that it can get down much closer to the extinction point.

I see that in your plots, the flame appears to be very close to the right inlet, and this might be causing issues. You may need to adjust your mass flow rates or domain width. Ideally you want that spike in the temperature field to be near the center of your domain.

@wandadars That might make a point. But if we do not know how to choose the best one, we cannot make the program choose the best. It's pretty difficult.

[wandadars]

Did you have anything in mind for reorganizing the 1D solver method @ischoegl ? Something like splitting up equations as a set of classes that can be looped over to evaluate? Like an Equation base class and then EnergyEquation, RadialMomentumEquation, ContinuityEquation, and these derived classes would be responsible for providing the residual evaluation at the interior of the domain. It definitely feels like much of the coding has content that is pushed together into large methods. I suppose the separate equation objects approach would cause issues with the boundary objects, because they would need to know which equation they are adjusting.

I'm just curious if any major adjustments are being brainstormed for the oneD coding. If large changes aren't coming any time soon, I will re-work the one and two point continuation coding into a smaller pull request without all those other changes that I made to the oneD classes regarding the equations.

In the current implementation the fuel boundary is assumed to be on the left and the oxidizer on the right, based on the coding in the Sim1D.cpp file for the setFuelInternalBoundary and setOxidizerInternalBoundary methods. Is generally the case for the setup of the counterflow configurations?

[ischoegl]

Did you have anything in mind for reorganizing the 1D solver method @ischoegl ? [...]

I have ideas, but I'm not convinced that I have the bandwidth in the near future. It would be nice to add equations dynamically during the setup process, but I am not convinced that there is much of a benefit short-term. If efforts like yours become more frequent (and we start seeing large blocks of 'reserved' but often unused entries) then it may make sense eventually.

I'm just curious if any major adjustments are being brainstormed for the oneD coding. If large changes aren't coming any time soon, I will re-work the one and two point continuation coding into a smaller pull request without all those other changes that I made to the oneD classes regarding the equations.

There really isn't much of a discussion at the moment. I personally tried to clean up some of the BC's and data entries so they can be 'reused' for extensions I am working on (mostly stuff that works like IonFlow, where I add custom code in a specialization). The work on Cantera itself is done, and I don't anticipate needing anything else in the near future. Please go ahead and implement what you have in mind within the current infrastructure. Methods are bloated as they are already, so it won't matter much if they grow a little more.

In the current implementation the fuel boundary is assumed to be on the left and the oxidizer on the right, based on the coding in the Sim1D.cpp file for the setFuelInternalBoundary and setOxidizerInternalBoundary methods. Is generally the case for the setup of the counterflow configurations?

There aren't any rules afaik, so fuel and oxidizer can be on either side - it would be good to either generalize or implement checks that users do what the code expects ...

[wandadars]

Work on this is shifted to #1622