|
| 1 | +import types |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +import yt |
| 5 | + |
| 6 | +LIBYT_VERSION = (0, 2, 0) |
| 7 | + |
| 8 | + |
| 9 | +def create_libyt_stub( |
| 10 | + simulation: str, |
| 11 | + fig_base_name: str, |
| 12 | + test_data: str, |
| 13 | + get_code_params: dict, |
| 14 | + field_list: dict, |
| 15 | + particle_list: dict, |
| 16 | + simulation_field_to_yt_field, |
| 17 | +) -> types.ModuleType: |
| 18 | + """ |
| 19 | + Returns a stub module that mimics libyt with a specific simulation. |
| 20 | +
|
| 21 | + :note: yt store the data structure in 3d, which is if the simulation is 2d, an additional dimension is added. |
| 22 | +
|
| 23 | + hierarchy is stored in 3d as well in a continuous array with 0 as the first grid, if dimensionality is 2, |
| 24 | + the third dimension is set to 1. |
| 25 | + grid_data stores the grid id starting from the index offset. |
| 26 | +
|
| 27 | + When extracting the data using yt, the data is always in 3d, so if the simulation is 2d, we need to extract it. |
| 28 | +
|
| 29 | + :param simulation: simulation name, e.g., "gamer", "enzo", etc. |
| 30 | + :param fig_base_name: figure base name |
| 31 | + :param test_data: the absolute path to the test data. |
| 32 | + :param get_code_params: the code parameters defined in the simulation frontend, and how to get it. |
| 33 | + { |
| 34 | + "code_params": [(param_name, default_value), ...], |
| 35 | + "method": (function to get the parameter) |
| 36 | + "expected_error": Exception type that is expected to be raised if the parameter is not found |
| 37 | + } |
| 38 | + :param field_list: libyt-v0.2 defined field list |
| 39 | + :param particle_list: libyt-v0.2 defined particle list |
| 40 | + :param simulation_field_to_yt_field: mapping field_list and particle_list name to yt field name |
| 41 | + """ |
| 42 | + # Mock libyt module based on libyt version 0.x.0 |
| 43 | + stub = types.ModuleType("libyt") |
| 44 | + stub.libyt_info = { |
| 45 | + "version": LIBYT_VERSION, |
| 46 | + "SERIAL_MODE": True, |
| 47 | + "INTERACTIVE_MODE": False, |
| 48 | + "JUPYTER_KERNEL": False, |
| 49 | + "SUPPORT_TIMER": False, |
| 50 | + } |
| 51 | + stub.param_yt = { |
| 52 | + "frontend": simulation, |
| 53 | + "fig_basename": fig_base_name, |
| 54 | + "current_time": None, |
| 55 | + "current_redshift": None, |
| 56 | + "omega_lambda": None, |
| 57 | + "omega_matter": None, |
| 58 | + "hubble_constant": None, |
| 59 | + "length_unit": None, |
| 60 | + "mass_unit": None, |
| 61 | + "time_unit": None, |
| 62 | + "magnetic_unit": None, |
| 63 | + "cosmological_simulation": None, |
| 64 | + "dimensionality": None, |
| 65 | + "refine_by": None, |
| 66 | + "velocity_unit": None, |
| 67 | + "domain_left_edge": None, |
| 68 | + "domain_right_edge": None, |
| 69 | + "periodicity": None, |
| 70 | + "domain_dimensions": None, |
| 71 | + "num_grids": None, |
| 72 | + "index_offset": None, |
| 73 | + "field_list": None, |
| 74 | + "particle_list": None, |
| 75 | + } |
| 76 | + stub.param_user = {} |
| 77 | + stub.hierarchy = { |
| 78 | + "grid_left_edge": None, |
| 79 | + "grid_right_edge": None, |
| 80 | + "grid_dimensions": None, |
| 81 | + "grid_parent_id": None, |
| 82 | + "grid_levels": None, |
| 83 | + "proc_num": None, |
| 84 | + # "par_count_list": None, # This is optional, which is a bad design choice. |
| 85 | + } |
| 86 | + stub.grid_data = {} |
| 87 | + stub.particle_data = {} |
| 88 | + |
| 89 | + ds = yt.load(test_data) |
| 90 | + # Fill in param_user |
| 91 | + for param in get_code_params["code_params"]: |
| 92 | + try: |
| 93 | + stub.param_user[param[0]] = get_code_params["method"](ds, param[0]) |
| 94 | + except get_code_params["expected_error"]: |
| 95 | + stub.param_user[param[0]] = param[1] |
| 96 | + |
| 97 | + # Fill in param_yt |
| 98 | + for param in ds.__dict__.keys(): |
| 99 | + if param in stub.param_yt and stub.param_yt[param] is None: |
| 100 | + stub.param_yt[param] = getattr(ds, param) |
| 101 | + if stub.param_yt["velocity_unit"] is None: |
| 102 | + stub.param_yt["velocity_unit"] = stub.param_yt["length_unit"] / stub.param_yt["time_unit"] |
| 103 | + stub.param_yt["domain_left_edge"] = ds.domain_left_edge |
| 104 | + stub.param_yt["domain_right_edge"] = ds.domain_right_edge |
| 105 | + stub.param_yt["periodicity"] = ds.periodicity |
| 106 | + stub.param_yt["domain_dimensions"] = ds.domain_dimensions |
| 107 | + stub.param_yt["num_grids"] = ds.index.num_grids |
| 108 | + stub.param_yt["index_offset"] = ds._index_class.grid._id_offset |
| 109 | + stub.param_yt["field_list"] = field_list |
| 110 | + stub.param_yt["particle_list"] = particle_list |
| 111 | + |
| 112 | + # Fill in hierarchy |
| 113 | + stub.hierarchy["grid_left_edge"] = ds.index.grid_left_edge |
| 114 | + stub.hierarchy["grid_right_edge"] = ds.index.grid_right_edge |
| 115 | + stub.hierarchy["grid_dimensions"] = ds.index.grid_dimensions |
| 116 | + stub.hierarchy["grid_levels"] = ds.index.grid_levels |
| 117 | + stub.hierarchy["proc_num"] = np.zeros((stub.param_yt["num_grids"], 1), dtype=np.int32) |
| 118 | + stub.hierarchy["grid_parent_id"] = np.ones((stub.param_yt["num_grids"], 1), dtype=np.int32) * -1 |
| 119 | + for g in ds.index.grids: |
| 120 | + if g.Parent is not None: |
| 121 | + stub.hierarchy["grid_parent_id"][g.id - stub.param_yt["index_offset"]] = g.Parent.id |
| 122 | + |
| 123 | + # Fill in grid_data |
| 124 | + for gid in range(stub.param_yt["num_grids"]): |
| 125 | + if gid + stub.param_yt["index_offset"] not in stub.grid_data: |
| 126 | + stub.grid_data[gid + stub.param_yt["index_offset"]] = {} |
| 127 | + for field in field_list.keys(): |
| 128 | + # TODO: assume cell-centered fields |
| 129 | + ghost_cells = field_list[field]["ghost_cell"] |
| 130 | + allocate_dim = stub.hierarchy["grid_dimensions"][gid][ |
| 131 | + : stub.param_yt["dimensionality"] |
| 132 | + ].copy() |
| 133 | + if field_list[field]["contiguous_in_x"]: |
| 134 | + allocate_dim = np.flip(allocate_dim) |
| 135 | + for d in range(2 * stub.param_yt["dimensionality"]): |
| 136 | + allocate_dim[int(d / 2)] += ghost_cells[d] |
| 137 | + stub.grid_data[gid + stub.param_yt["index_offset"]][field] = ( |
| 138 | + np.ones(allocate_dim, dtype=np.float64) * np.nan |
| 139 | + ) |
| 140 | + |
| 141 | + if field_list[field]["contiguous_in_x"]: |
| 142 | + if stub.param_yt["dimensionality"] == 3: |
| 143 | + stub.grid_data[gid + stub.param_yt["index_offset"]][field][ |
| 144 | + ghost_cells[0] : allocate_dim[0] - ghost_cells[1], |
| 145 | + ghost_cells[2] : allocate_dim[1] - ghost_cells[3], |
| 146 | + ghost_cells[4] : allocate_dim[2] - ghost_cells[5], |
| 147 | + ] = ( |
| 148 | + ds.index.grids[gid][simulation_field_to_yt_field[field]] |
| 149 | + .swapaxes(0, 2) |
| 150 | + .in_base("code") |
| 151 | + ) |
| 152 | + elif stub.param_yt["dimensionality"] == 2: |
| 153 | + stub.grid_data[gid + stub.param_yt["index_offset"]][field][ |
| 154 | + ghost_cells[0] : allocate_dim[0] - ghost_cells[1], |
| 155 | + ghost_cells[2] : allocate_dim[1] - ghost_cells[3], |
| 156 | + ] = np.squeeze( |
| 157 | + ds.index.grids[gid][simulation_field_to_yt_field[field]].in_base("code") |
| 158 | + ).swapaxes( |
| 159 | + 0, 1 |
| 160 | + ) |
| 161 | + else: |
| 162 | + stub.grid_data[gid + stub.param_yt["index_offset"]][field][ |
| 163 | + ghost_cells[0] : allocate_dim[0] - ghost_cells[1] |
| 164 | + ] = ds.index.grids[gid][simulation_field_to_yt_field[field]].in_base("code") |
| 165 | + |
| 166 | + return stub |
0 commit comments