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

# 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

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")