|
| 1 | +import h5py |
| 2 | +import torch |
| 3 | +import numpy as np |
| 4 | +import cantera as ct |
| 5 | + |
| 6 | +from dfode_kit.data.contracts import MECHANISM_ATTR, require_h5_attr, read_scalar_field_datasets |
| 7 | +from dfode_kit.data.io_hdf5 import get_TPY_from_h5, touch_h5 |
| 8 | +from dfode_kit.utils import BCT, inverse_BCT |
| 9 | + |
| 10 | + |
| 11 | +def advance_reactor(gas, state, reactor, reactor_net, time_step): |
| 12 | + """Advance the reactor simulation for a given state.""" |
| 13 | + state = state.flatten() |
| 14 | + |
| 15 | + expected_shape = (2 + gas.n_species,) |
| 16 | + if state.shape != expected_shape: |
| 17 | + raise ValueError( |
| 18 | + f"Expected state shape {expected_shape}, got {state.shape}" |
| 19 | + ) |
| 20 | + |
| 21 | + gas.TPY = state[0], state[1], state[2:] |
| 22 | + |
| 23 | + reactor.syncState() |
| 24 | + reactor_net.reinitialize() |
| 25 | + reactor_net.advance(time_step) |
| 26 | + reactor_net.set_initial_time(0.0) |
| 27 | + |
| 28 | + return gas |
| 29 | + |
| 30 | + |
| 31 | +@torch.no_grad() |
| 32 | +def load_model(model_path, device, model_class, model_layers): |
| 33 | + state_dict = torch.load(model_path, map_location='cpu') |
| 34 | + |
| 35 | + model = model_class(model_layers) |
| 36 | + model.load_state_dict(state_dict['net']) |
| 37 | + |
| 38 | + model.eval() |
| 39 | + model.to(device=device) |
| 40 | + |
| 41 | + return model |
| 42 | + |
| 43 | + |
| 44 | +@torch.no_grad() |
| 45 | +def predict_Y(model, model_path, d_arr, mech, device): |
| 46 | + gas = ct.Solution(mech) |
| 47 | + n_species = gas.n_species |
| 48 | + expected_dims = 2 + n_species |
| 49 | + if d_arr.shape[1] != expected_dims: |
| 50 | + raise ValueError( |
| 51 | + f"Expected input with {expected_dims} columns, got {d_arr.shape[1]}" |
| 52 | + ) |
| 53 | + |
| 54 | + state_dict = torch.load(model_path, map_location='cpu') |
| 55 | + |
| 56 | + Xmu0 = state_dict['data_in_mean'] |
| 57 | + Xstd0 = state_dict['data_in_std'] |
| 58 | + Ymu0 = state_dict['data_target_mean'] |
| 59 | + Ystd0 = state_dict['data_target_std'] |
| 60 | + |
| 61 | + d_arr = np.clip(d_arr, 0, None) |
| 62 | + d_arr[:, 1] *= 0 |
| 63 | + d_arr[:, 1] += 101325 |
| 64 | + |
| 65 | + orig_Y = d_arr[:, 2:].copy() |
| 66 | + in_bct = d_arr.copy() |
| 67 | + in_bct[:, 2:] = BCT(in_bct[:, 2:]) |
| 68 | + in_bct_norm = (in_bct - Xmu0) / Xstd0 |
| 69 | + |
| 70 | + input = torch.from_numpy(in_bct_norm).float().to(device=device) |
| 71 | + |
| 72 | + output = model(input) |
| 73 | + |
| 74 | + out_bct = output.cpu().numpy() * Ystd0 + Ymu0 + in_bct[:, 2:-1] |
| 75 | + next_Y = orig_Y.copy() |
| 76 | + next_Y[:, :-1] = inverse_BCT(out_bct) |
| 77 | + next_Y[:, :-1] = next_Y[:, :-1] / np.sum(next_Y[:, :-1], axis=1, keepdims=True) * (1 - next_Y[:, -1:]) |
| 78 | + |
| 79 | + return next_Y |
| 80 | + |
| 81 | + |
| 82 | +@torch.no_grad() |
| 83 | +def nn_integrate(orig_arr, model_path, device, model_class, model_layers, time_step, mech, frozen_temperature=510): |
| 84 | + model = load_model(model_path, device, model_class, model_layers) |
| 85 | + |
| 86 | + mask = orig_arr[:, 0] > frozen_temperature |
| 87 | + infer_arr = orig_arr[mask, :] |
| 88 | + |
| 89 | + next_Y = predict_Y(model, model_path, infer_arr, mech, device) |
| 90 | + |
| 91 | + new_states = np.hstack((np.zeros((orig_arr.shape[0], 1)), orig_arr)) |
| 92 | + new_states[:, 0] += time_step |
| 93 | + new_states[:, 2] = orig_arr[:, 1] |
| 94 | + new_states[mask, 3:] = next_Y |
| 95 | + |
| 96 | + setter_gas = ct.Solution(mech) |
| 97 | + getter_gas = ct.Solution(mech) |
| 98 | + new_T = np.zeros_like(next_Y[:, 0]) |
| 99 | + |
| 100 | + for idx, (state, next_y) in enumerate(zip(infer_arr, next_Y)): |
| 101 | + try: |
| 102 | + setter_gas.TPY = state[0], state[1], state[2:] |
| 103 | + h = setter_gas.enthalpy_mass |
| 104 | + |
| 105 | + getter_gas.Y = next_y |
| 106 | + getter_gas.HP = h, state[1] |
| 107 | + |
| 108 | + new_T[idx] = getter_gas.T |
| 109 | + |
| 110 | + except ct.CanteraError: |
| 111 | + continue |
| 112 | + new_states[mask, 1] = new_T |
| 113 | + |
| 114 | + return new_states |
| 115 | + |
| 116 | + |
| 117 | +def integrate_h5( |
| 118 | + file_path, |
| 119 | + save_path1, |
| 120 | + save_path2, |
| 121 | + time_step, |
| 122 | + cvode_integration=True, |
| 123 | + nn_integration=False, |
| 124 | + model_settings=None, |
| 125 | +): |
| 126 | + """Process scalar-field datasets and save CVODE / NN integration outputs.""" |
| 127 | + with h5py.File(file_path, 'r') as f: |
| 128 | + mech = require_h5_attr(f, MECHANISM_ATTR) |
| 129 | + |
| 130 | + data_dict = read_scalar_field_datasets(file_path) |
| 131 | + |
| 132 | + if cvode_integration: |
| 133 | + gas = ct.Solution(mech) |
| 134 | + reactor = ct.Reactor(gas, name='Reactor1', energy='off') |
| 135 | + reactor_net = ct.ReactorNet([reactor]) |
| 136 | + reactor_net.rtol, reactor_net.atol = 1e-6, 1e-10 |
| 137 | + |
| 138 | + processed_data_dict = {} |
| 139 | + |
| 140 | + for name, data in data_dict.items(): |
| 141 | + processed_data = np.empty((data.shape[0], data.shape[1] + 1)) |
| 142 | + for i, state in enumerate(data): |
| 143 | + gas = advance_reactor(gas, state, reactor, reactor_net, time_step) |
| 144 | + |
| 145 | + new_state = np.array([time_step, gas.T, gas.P] + list(gas.Y)) |
| 146 | + |
| 147 | + processed_data[i, :] = new_state |
| 148 | + |
| 149 | + processed_data_dict[name] = processed_data |
| 150 | + |
| 151 | + with h5py.File(save_path1, 'a') as f: |
| 152 | + cvode_group = f.create_group('cvode_integration') |
| 153 | + |
| 154 | + for dataset_name, processed_data in processed_data_dict.items(): |
| 155 | + cvode_group.create_dataset(dataset_name, data=processed_data) |
| 156 | + print(f'Saved processed dataset: {dataset_name} in cvode_integration group') |
| 157 | + |
| 158 | + if nn_integration: |
| 159 | + processed_data_dict = {} |
| 160 | + if model_settings is None: |
| 161 | + raise ValueError("model_settings must be provided for neural network integration.") |
| 162 | + |
| 163 | + for name, data in data_dict.items(): |
| 164 | + try: |
| 165 | + processed_data = nn_integrate(data, **model_settings) |
| 166 | + processed_data_dict[name] = processed_data |
| 167 | + except Exception as e: |
| 168 | + print(f"Error processing dataset '{name}': {e}") |
| 169 | + |
| 170 | + with h5py.File(save_path2, 'a') as f: |
| 171 | + if 'nn_integration' in f: |
| 172 | + del f['nn_integration'] |
| 173 | + nn_group = f.create_group('nn_integration') |
| 174 | + |
| 175 | + for dataset_name, processed_data in processed_data_dict.items(): |
| 176 | + nn_group.create_dataset(dataset_name, data=processed_data) |
| 177 | + print(f'Saved processed dataset: {dataset_name} in nn_integration group') |
| 178 | + |
| 179 | + |
| 180 | +def calculate_error( |
| 181 | + mech_path, |
| 182 | + save_path1, |
| 183 | + save_path2, |
| 184 | + error='RMSE' |
| 185 | +): |
| 186 | + gas = ct.Solution(mech_path) |
| 187 | + |
| 188 | + with h5py.File(save_path1, 'r') as f1, h5py.File(save_path2, 'r') as f2: |
| 189 | + cvode_group = f1['cvode_integration'] |
| 190 | + nn_group = f2['nn_integration'] |
| 191 | + |
| 192 | + common_datasets = set(cvode_group.keys()) & set(nn_group.keys()) |
| 193 | + |
| 194 | + sorted_datasets = sorted(common_datasets, key=lambda x: float(x)) |
| 195 | + results = {} |
| 196 | + |
| 197 | + for ds_name in sorted_datasets: |
| 198 | + cvode_data = cvode_group[ds_name][:, 3:] |
| 199 | + nn_data = nn_group[ds_name][:, 3:] |
| 200 | + |
| 201 | + if error == "RMSE": |
| 202 | + rmse_per_dim = np.sqrt(np.mean((cvode_data - nn_data) ** 2, axis=0)) |
| 203 | + results[ds_name] = rmse_per_dim |
| 204 | + |
| 205 | + print(f"RMSE of ataset: {ds_name}") |
| 206 | + for dim_idx, rmse_val in enumerate(rmse_per_dim, start=1): |
| 207 | + id = gas.species_names[dim_idx - 3] |
| 208 | + print(f" Species {id}: {rmse_val:.6e}") |
| 209 | + print() |
| 210 | + |
| 211 | + return results |
0 commit comments