-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathintegrator_selector.py
More file actions
360 lines (311 loc) · 12.3 KB
/
integrator_selector.py
File metadata and controls
360 lines (311 loc) · 12.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
import numpy as np
import cantera as ct
import time
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import Dict, List, Tuple, Optional
from ember import Config, Paths, InitialCondition, StrainParameters, General, \
Times, TerminationCondition, ConcreteConfig, Debug, RK23Tolerances, QssTolerances, CvodeTolerances, _ember
import os
@dataclass
class ToleranceSettings:
"""Settings for different integrator tolerances"""
# CVODE settings
cvode_rtol: float = 1e-6
cvode_momentum_atol: float = 1e-7
cvode_energy_atol: float = 1e-8
cvode_species_atol: float = 1e-13
cvode_min_timestep: float = 1e-18
cvode_max_timestep: float = 1e-5
# RK23 settings
rk_rtol: float = 1e-6
rk_atol: float = 1e-8
rk_min_timestep: float = 1e-10
rk_max_timestep: float = 1e-4
rk_max_steps: int = 100000
# QSS settings
qss_eps_min: float = 1e-2
qss_eps_max: float = 1e1
qss_min_timestep: float = 1e-16
qss_max_timestep: float = 1e-6
qss_atol: float = 1e-11
qss_min_val: float = 1e-60
@dataclass
class SimulationSettings:
"""General simulation settings"""
output_dir: str = 'run/diffusion_benchmark'
n_threads: int = 2
global_timestep: float = 1e-6
profile_interval: int = 20
t_end: float = 0.08
n_points: int = 100
x_left: float = -0.02
x_right: float = 0.02
T_fuel: float = 600
T_oxidizer: float = 1200
pressure: float = 101325
strain_rate: float = 100
@dataclass
class IntegratorConfig:
name: str
config: ConcreteConfig
constant_integrator: Optional[str] # For constant integration mode
color: str # For plotting
@dataclass
class SimulationResult:
x: np.ndarray # Spatial coordinates
times: np.ndarray
temperatures: np.ndarray
species: Dict[int, np.ndarray]
total_time: float
integrator_distribution: Dict[str, int] # Count of each integrator type used
def create_base_config(sim_settings: SimulationSettings, tol_settings: ToleranceSettings) -> Config:
"""Create base configuration with given settings"""
return Config(
Paths(outputDir=sim_settings.output_dir),
General(nThreads=sim_settings.n_threads, chemistryIntegrator='cvode'),
InitialCondition(
Tfuel=sim_settings.T_fuel,
Toxidizer=sim_settings.T_oxidizer,
centerWidth=0.0,
equilibrateCounterflow=False,
flameType='diffusion',
slopeWidth=0.0,
xLeft=sim_settings.x_left,
xRight=sim_settings.x_right,
pressure=sim_settings.pressure,
nPoints=sim_settings.n_points
),
StrainParameters(final=sim_settings.strain_rate, initial=sim_settings.strain_rate),
Times(globalTimestep=sim_settings.global_timestep,
profileStepInterval=sim_settings.profile_interval),
TerminationCondition(
abstol=0.0,
dTdtTol=0,
steadyPeriod=1.0,
tEnd=sim_settings.t_end,
tolerance=0.0
),
RK23Tolerances(
absoluteTolerance=tol_settings.rk_atol,
relativeTolerance=tol_settings.rk_rtol,
minimumTimestep=tol_settings.rk_min_timestep,
maximumTimestep=tol_settings.rk_max_timestep,
maxStepsNumber=tol_settings.rk_max_steps
),
QssTolerances(
epsmin=tol_settings.qss_eps_min,
epsmax=tol_settings.qss_eps_max,
dtmin=tol_settings.qss_min_timestep,
dtmax=tol_settings.qss_max_timestep,
iterationCount=1,
abstol=tol_settings.qss_atol,
minval=tol_settings.qss_min_val,
stabilityCheck=False
),
CvodeTolerances(
relativeTolerance=tol_settings.cvode_rtol,
momentumAbsTol=tol_settings.cvode_momentum_atol,
energyAbsTol=tol_settings.cvode_energy_atol,
speciesAbsTol=tol_settings.cvode_species_atol,
minimumTimestep=tol_settings.cvode_min_timestep,
maximumTimestep=tol_settings.cvode_max_timestep
),
Debug(veryVerbose=False)
)
def create_configs(sim_settings: SimulationSettings) -> Tuple[IntegratorConfig, List[IntegratorConfig]]:
"""Create configurations for baseline and test integrators"""
# Create stricter tolerance settings for baseline
baseline_tol = ToleranceSettings(
cvode_rtol=1e-8,
cvode_momentum_atol=1e-8,
cvode_energy_atol=1e-8,
cvode_species_atol=1e-13
)
# Regular tolerance settings for test cases
test_tol = ToleranceSettings()
# Baseline configuration (CVODE with tight tolerances)
baseline_config = IntegratorConfig(
name="CVODE-baseline",
config=ConcreteConfig(create_base_config(sim_settings, baseline_tol)),
constant_integrator='cvode',
color='k'
)
# Test configurations
base_config = create_base_config(sim_settings, test_tol)
test_configs = [
IntegratorConfig(
name="CVODE",
config=ConcreteConfig(base_config),
constant_integrator='cvode',
color='b'
),
IntegratorConfig(
name="BoostRK",
config=ConcreteConfig(base_config),
constant_integrator='boostRK',
color='g'
),
IntegratorConfig(
name="QSS",
config=ConcreteConfig(base_config),
constant_integrator='qss',
color='r'
),
IntegratorConfig(
name="Heuristic",
config=ConcreteConfig(base_config),
constant_integrator=None, # None means use heuristic selection
color='m'
)
]
return baseline_config, test_configs
def setup_species_indices() -> Dict[int, str]:
"""Setup the species indices to track"""
gas = ct.Solution('gri30.yaml')
names = gas.species_names
return {
names.index('CH4'): 'CH4',
names.index('O2'): 'O2',
names.index('CO2'): 'CO2',
names.index('H2O'): 'H2O',
names.index('CO'): 'CO',
names.index('OH'): 'OH'
}
def set_integrator_heuristic(solver) -> List[str]:
"""Set integrator type based on temperature and equivalence ratio"""
nPoints = len(solver.T)
# Start with temperature-based decision
integ = np.where(solver.T <= 600.0, 'boostRK', 'cvode')
try:
# Get equivalence ratio
phi = solver.phi
# Use boostRK for extreme conditions
integ = np.where(phi == -1, 'boostRK', integ) # invalid phi
integ = np.where(phi <= 1e-8, 'boostRK', integ) # oxidizer-dominated
integ = np.where(phi >= 1e4, 'boostRK', integ) # fuel-dominated
# Create boolean mask for CVODE points
cvode_mask = (integ == 'cvode')
# Include neighboring points
cvode_mask_left = np.roll(cvode_mask, 1)
cvode_mask_right = np.roll(cvode_mask, -1)
cvode_mask_left[0] = False
cvode_mask_right[-1] = False
use_cvode = cvode_mask | cvode_mask_left | cvode_mask_right
integ = np.where(use_cvode, 'cvode', 'boostRK')
except Exception as e:
print(f"Warning: Could not calculate phi for integrator selection: {e}")
return integ.tolist()
def run_simulation(config: IntegratorConfig) -> SimulationResult:
"""Run simulation with given configuration and return results"""
solver = _ember.FlameSolver(config.config)
solver.initialize()
start_time = time.time()
integrator_counts = {'cvode': 0, 'boostRK': 0, 'qss': 0}
done = False
times = []
while not done:
if config.constant_integrator:
integrator_types = [config.constant_integrator] * len(solver.T)
else:
integrator_types = set_integrator_heuristic(solver)
# Update integrator counts
for int_type in integrator_types:
integrator_counts[int_type] += 1
solver.set_integrator_types(integrator_types)
done = solver.step()
times.append(solver.t)
total_time = time.time() - start_time
return SimulationResult(
x=solver.x,
times=np.array([solver.t]), # Single time point for steady-state
temperatures=solver.T,
species={idx: solver.Y[idx] for idx in setup_species_indices()},
total_time=total_time,
integrator_distribution=integrator_counts
)
def calculate_errors(baseline: SimulationResult, result: SimulationResult) -> Dict[str, float]:
"""Calculate various error metrics compared to baseline"""
# Interpolate result to baseline spatial points for comparison
temp_interp = np.interp(baseline.x, result.x, result.temperatures)
errors = {
'temp_max_error': np.max(np.abs((temp_interp - baseline.temperatures) / baseline.temperatures)),
'temp_mean_error': np.mean(np.abs((temp_interp - baseline.temperatures) / baseline.temperatures)),
'temp_rms_error': np.sqrt(np.mean(((temp_interp - baseline.temperatures) / baseline.temperatures) ** 2))
}
return errors
def plot_results(baseline: SimulationResult, results: Dict[str, SimulationResult],
species_indices: Dict[int, str], configs: List[IntegratorConfig],
output_dir: str):
"""Create comparison plots"""
# Temperature plot
plt.figure(figsize=(10, 6))
plt.plot(baseline.x, baseline.temperatures, 'k--', label='Baseline')
for name, result in results.items():
config = next(c for c in configs if c.name == name)
plt.plot(result.x, result.temperatures, color=config.color, label=name)
plt.xlabel('Position (m)')
plt.ylabel('Temperature (K)')
plt.legend()
plt.savefig(os.path.join(output_dir, 'temperature_comparison_1d.png'))
plt.close()
# Species plots
plt.figure(figsize=(12, 8))
for idx, species_name in species_indices.items():
plt.subplot(2, 3, list(species_indices.keys()).index(idx) + 1)
plt.plot(baseline.x, baseline.species[idx], 'k--', label='Baseline')
for name, result in results.items():
config = next(c for c in configs if c.name == name)
plt.plot(result.x, result.species[idx], color=config.color, label=name)
plt.xlabel('Position (m)')
plt.ylabel(f'{species_name} Mass Fraction')
if idx == list(species_indices.keys())[0]: # Only show legend once
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'species_comparison_1d.png'))
plt.close()
def main():
# Create simulation settings
sim_settings = SimulationSettings(
output_dir='run/diffusion_benchmark_test',
n_threads=2,
global_timestep=1e-6,
profile_interval=20,
t_end=0.08,
n_points=100
)
# Setup
baseline_config, test_configs = create_configs(sim_settings)
species_indices = setup_species_indices()
# Create output directory if it doesn't exist
os.makedirs(sim_settings.output_dir, exist_ok=True)
# Run baseline
print("Running baseline simulation...")
baseline_result = run_simulation(baseline_config)
# Run test cases
results = {}
errors = {}
print("\nRunning test simulations...")
for config in test_configs:
print(f"\nRunning {config.name}...")
result = run_simulation(config)
results[config.name] = result
errors[config.name] = calculate_errors(baseline_result, result)
# Print performance comparison
print("\nPerformance Comparison:")
print(f"{'Integrator':<15} {'Total Time (s)':<15} {'Max Error (%)':<15} {'Mean Error (%)':<15} {'RMS Error (%)':<15}")
print("-" * 75)
for name, result in results.items():
error = errors[name]
print(f"{name:<15} {result.total_time:<15.3f} "
f"{error['temp_max_error']*100:<15.3f} "
f"{error['temp_mean_error']*100:<15.3f} "
f"{error['temp_rms_error']*100:<15.3f}")
# Print integrator distribution for heuristic case
if 'Heuristic' in results:
print("\nHeuristic Integrator Distribution:")
dist = results
print(dist)
plot_results(baseline_result, results, species_indices, test_configs, sim_settings.output_dir)
if __name__ == "__main__":
main()