Missing files
Kofi Abusah
David Ketcheson
Kofi Abusah
import numpy as np
from clawpack.pyclaw import controller
from clawpack.pyclaw import ClawSolver2D
from clawpack.pyclaw import solution
from clawpack.pyclaw import state # Import the state module to access State
from clawpack.pyclaw.geometry import Dimension, Grid, Domain
from clawpack.pyclaw.solver import BC # Import BC directly from the solver module
# Corrected: Assuming 'riemann_reactive_euler_2D' is the name of the module created by f2py
# and that it contains a function with the same name.
from riemann_reactive_euler_2D import riemann_reactive_euler_2D
# Define the number of conserved quantities
NUM_EQ = 5 # Density, momentum (x,y), energy, reactant mass fraction
# Physical constants (adjust as needed for MMH/RFNA)
GAMMA = 1.4 # Ratio of specific heats for ideal gas
R = 287.05 # Gas constant
# Define the source term function for the chemical reaction
def reactive_source(solver, state_obj, dt):
"""
Source term for chemical reactions. This is a placeholder and
needs to be replaced with a realistic reaction rate model for MMH/RFNA.
This function modifies state_obj.q in place.
"""
rho = state_obj.q[0, :, :]
e = state_obj.q[NUM_EQ - 1, :, :] # Energy (assuming last equation is energy)
Y = state_obj.q[NUM_EQ - 1, :, :] # Reactant mass fraction (assuming last equation is reactant)
# Placeholder reaction rate (e.g., Arrhenius kinetics)
k = 1e9 * np.exp(-10000 / (e / (rho * R + 1e-9))) # Avoid division by zero
dYdt = -k * Y # Example rate of consumption of reactant
# Update energy due to heat release (Q is heat of reaction)
Q = 5e6 # Placeholder heat of reaction (J/kg)
dedt = Q * dYdt # Example energy change
# Apply the source term update to state_obj.q
state_obj.q[3, :, :] += dedt * dt # Update energy
state_obj.q[4, :, :] += dYdt * dt # Update reactant mass fraction
# This Python Riemann solver is explicitly not used when kernel_language is 'fortran'.
# Its existence here is for context, but the Fortran version will be called.
# def reactive_euler_riemann(q_l, q_r, aux_l, aux_r, problem_data):
# raise NotImplementedError("This Python Riemann solver is not used when kernel_language is 'fortran' in multi-D.")
# Create the PyClaw solution object
def create_solution(domain_obj, initial_conditions):
"""
Creates the PyClaw Solution object with the initial conditions.
"""
state_obj = state.State(domain_obj, NUM_EQ)
state_obj.problem_data['gamma'] = GAMMA
state_obj.q = initial_conditions
return solution.Solution(state_obj, domain_obj)
# Define the solver
def create_solver():
"""
Creates the PyClaw solver object configured to use Fortran kernels.
"""
# Pass the actual compiled Fortran Riemann solver function to ClawSolver2D.
solver = ClawSolver2D(riemann_reactive_euler_2D)
solver.kernel_language = 'fortran' # Crucial for multi-D problems
solver.num_waves = NUM_EQ
# Source term handling: assign the function to step_source
solver.step_source = reactive_source
# Set boundary conditions (e.g., wall boundary, outflow).
solver.bc_lower = {0: BC.wall, 1: BC.wall}
solver.bc_upper = {0: BC.extrap, 1: BC.extrap}
return solver
# Main simulation function
def simulate_rde():
"""
Runs the RDE simulation using PyClaw.
"""
mx = 100
my = 100
x_lower = 0.0
x_upper = 1.0
y_lower = 0.0
y_upper = 0.1
x = Dimension(x_lower, x_upper, mx, name='x')
y = Dimension(y_lower, y_upper, my, name='y')
domain_obj = Domain((x, y))
initial_conditions = np.zeros((NUM_EQ, mx, my))
initial_conditions[0, :, :] = 1.0
initial_conditions[1, :, :] = 0.0
initial_conditions[2, :, :] = 0.0
initial_conditions[3, :, :] = 2.5
initial_conditions[4, :, :] = 1.0
sol = create_solution(domain_obj, initial_conditions)
solver = create_solver()
claw = controller.Controller()
claw.keep_copy = True
claw.num_output_times = 10
claw.solution = sol
claw.solver = solver
claw.tfinal = 0.001
status = claw.run()
if __name__ == '__main__':