From 59bac0e73aa16e2f1cd734124e1c14628251c939 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Wed, 23 Apr 2025 10:23:02 +0100 Subject: [PATCH 1/6] first working version of RT fields in MEGATRON --- yt/frontends/ramses/data_structures.py | 4 +++ yt/frontends/ramses/field_handlers.py | 7 +++++ yt/frontends/ramses/io_utils.pyx | 39 ++++++++++++++++++-------- 3 files changed, 39 insertions(+), 11 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 4918d31ea2a..38cf6992cfc 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -429,6 +429,9 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): # Initializing data container for field in fields: data[field] = np.zeros(cell_count, "float64") + + # Get the size of the data, single precision if RT field, double otherwise + field_size = "single" if any(field.startswith("Photon") for field in fields) else "double" # Do an early exit if the cell count is null if cell_count == 0: @@ -452,6 +455,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): fields, data, oct_handler, + field_size # Hardcoded for now ) return data diff --git a/yt/frontends/ramses/field_handlers.py b/yt/frontends/ramses/field_handlers.py index fafd18b5967..087e1cb95a8 100644 --- a/yt/frontends/ramses/field_handlers.py +++ b/yt/frontends/ramses/field_handlers.py @@ -251,6 +251,11 @@ def offset(self): *structure* of your fluid file is non-canonical, change this. """ nvars = len(self._detected_field_list[self.ds.unique_identifier]) + + first_field = self.field_list[0][1] + # Get the size of the data, single precision if RT field, double otherwise + field_size = "single" if first_field.startswith("Photon") else "double" + with FortranFile(self.fname) as fd: # Skip headers nskip = len(self.attrs) @@ -269,6 +274,7 @@ def offset(self): # # So there are 8 * nvars records each with length (nocts, ) # at each (level, cpus) + offset, level_count = read_offset( fd, min_level, @@ -276,6 +282,7 @@ def offset(self): self.parameters[self.ds.unique_identifier]["nvar"], self.domain.amr_header, Nskip=nvars * 8, + size=field_size, ) self._level_count = level_count diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index 643099a1704..b55ea023c06 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -19,9 +19,13 @@ ctypedef np.float64_t DOUBLE_t cdef int INT32_SIZE = sizeof(np.int32_t) cdef int INT64_SIZE = sizeof(np.int64_t) cdef int DOUBLE_SIZE = sizeof(np.float64_t) +cdef int SINGLE_SIZE = sizeof(np.float32_t) -cdef inline int skip_len(int Nskip, int record_len) noexcept nogil: +cdef inline int skip_len(int Nskip, int record_len, str size="double") noexcept nogil: + # If the data is single precision, we need to skip 4 bytes + if size == "single": + return Nskip * (record_len * SINGLE_SIZE + INT64_SIZE) return Nskip * (record_len * DOUBLE_SIZE + INT64_SIZE) @@ -111,7 +115,7 @@ def read_amr(FortranFile f, dict headers, @cython.wraparound(False) @cython.cdivision(True) @cython.nonecheck(False) -cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip): +cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip, str size="double"): cdef np.ndarray[np.int64_t, ndim=2] offset, level_count cdef INT64_t ndim, twotondim, nlevelmax, n_levels, nboundary, ncpu, ncpu_and_bound @@ -153,7 +157,7 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n if ilevel >= min_level: offset_view[icpu, ilevel - min_level] = f.tell() level_count_view[icpu, ilevel - min_level] = file_ncache - f.seek(skip_len(Nskip, file_ncache), SEEK_CUR) + f.seek(skip_len(Nskip, file_ncache, size), SEEK_CUR) return offset, level_count @@ -172,7 +176,9 @@ def fill_hydro(FortranFile f, INT64_t ndim, list all_fields, list fields, dict tr, RAMSESOctreeContainer oct_handler, - np.ndarray[np.int32_t, ndim=1] domain_inds=np.array([], dtype='int32')): + str size, + np.ndarray[np.int32_t, ndim=1] domain_inds=np.array([], dtype='int32'), + ): cdef INT64_t offset cdef dict tmp cdef str field @@ -197,7 +203,8 @@ def fill_hydro(FortranFile f, # The ordering is very important here, as we'll write directly into the memory # address the content of the files. - cdef np.float64_t[::1, :, :] buffer + cdef np.float32_t[::1, :, :] buffer_single + cdef np.float64_t[::1, :, :] buffer_double jump_len = 0 j = 0 @@ -211,7 +218,11 @@ def fill_hydro(FortranFile f, jumps[j] = jump_len cdef int first_field_index = jumps[0] - buffer = np.empty((level_count.max(), twotondim, nfields_selected), dtype="float64", order='F') + # The buffer is different depending on the size of the data + if size == "single": + buffer_single = np.empty((level_count.max(), twotondim, nfields_selected), dtype="float32", order='F') + else: + buffer_double = np.empty((level_count.max(), twotondim, nfields_selected), dtype="float64", order='F') # Precompute which levels we need to read Ncells = len(level_inds) @@ -231,7 +242,7 @@ def fill_hydro(FortranFile f, offset = offsets[icpu, ilevel] if offset == -1: continue - f.seek(offset + skip_len(first_field_index, nc), SEEK_SET) + f.seek(offset + skip_len(first_field_index, nc, size), SEEK_SET) # We have already skipped the first fields (if any) # so we "rewind" (this will cancel the first seek) @@ -241,9 +252,12 @@ def fill_hydro(FortranFile f, for j in range(nfields_selected): jump_len += jumps[j] if jump_len > 0: - f.seek(skip_len(jump_len, nc), SEEK_CUR) + f.seek(skip_len(jump_len, nc, size), SEEK_CUR) jump_len = 0 - f.read_vector_inplace('d', &buffer[0, i, j]) + if size == "single": + f.read_vector_inplace('f', &buffer_single[0, i, j]) + else: + f.read_vector_inplace('d', &buffer_double[0, i, j]) jump_len += jumps[nfields_selected] @@ -251,12 +265,15 @@ def fill_hydro(FortranFile f, # but since we're doing an absolute seek at the beginning of # the loop on CPUs, we can spare one seek here ## if jump_len > 0: - ## f.seek(skip_len(jump_len, nc), SEEK_CUR) + ## f.seek(skip_len(jump_len, nc, size), SEEK_CUR) # Alias buffer into dictionary tmp = {} for i, field in enumerate(fields): - tmp[field] = buffer[:, :, i] + if size == "single": + tmp[field] = np.asarray(buffer_single[:, :, i], dtype="float64") + else: + tmp[field] = buffer_double[:, :, i] if ncpu_selected > 1: oct_handler.fill_level_with_domain( From 88a95b97aff3dfd048f8adef68eb3f465475aa46 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Wed, 23 Apr 2025 10:53:28 +0100 Subject: [PATCH 2/6] changed field check for single precision to a boolean instead of string --- yt/frontends/ramses/data_structures.py | 4 ++-- yt/frontends/ramses/field_handlers.py | 4 ++-- yt/frontends/ramses/io_utils.pyx | 22 +++++++++++----------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index 38cf6992cfc..b8926931baa 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -431,7 +431,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): data[field] = np.zeros(cell_count, "float64") # Get the size of the data, single precision if RT field, double otherwise - field_size = "single" if any(field.startswith("Photon") for field in fields) else "double" + is_single = True if any(field.startswith("Photon") for field in fields) else False # Do an early exit if the cell count is null if cell_count == 0: @@ -455,7 +455,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): fields, data, oct_handler, - field_size # Hardcoded for now + is_single # Hardcoded for now ) return data diff --git a/yt/frontends/ramses/field_handlers.py b/yt/frontends/ramses/field_handlers.py index 087e1cb95a8..587a45a0495 100644 --- a/yt/frontends/ramses/field_handlers.py +++ b/yt/frontends/ramses/field_handlers.py @@ -254,7 +254,7 @@ def offset(self): first_field = self.field_list[0][1] # Get the size of the data, single precision if RT field, double otherwise - field_size = "single" if first_field.startswith("Photon") else "double" + is_single = True if first_field.startswith("Photon") else False with FortranFile(self.fname) as fd: # Skip headers @@ -282,7 +282,7 @@ def offset(self): self.parameters[self.ds.unique_identifier]["nvar"], self.domain.amr_header, Nskip=nvars * 8, - size=field_size, + single=is_single, ) self._level_count = level_count diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index b55ea023c06..a740f62529c 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -22,9 +22,9 @@ cdef int DOUBLE_SIZE = sizeof(np.float64_t) cdef int SINGLE_SIZE = sizeof(np.float32_t) -cdef inline int skip_len(int Nskip, int record_len, str size="double") noexcept nogil: +cdef inline int skip_len(int Nskip, int record_len, bint single) noexcept nogil: # If the data is single precision, we need to skip 4 bytes - if size == "single": + if single: return Nskip * (record_len * SINGLE_SIZE + INT64_SIZE) return Nskip * (record_len * DOUBLE_SIZE + INT64_SIZE) @@ -115,7 +115,7 @@ def read_amr(FortranFile f, dict headers, @cython.wraparound(False) @cython.cdivision(True) @cython.nonecheck(False) -cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip, str size="double"): +cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip, bint single): cdef np.ndarray[np.int64_t, ndim=2] offset, level_count cdef INT64_t ndim, twotondim, nlevelmax, n_levels, nboundary, ncpu, ncpu_and_bound @@ -157,7 +157,7 @@ cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t n if ilevel >= min_level: offset_view[icpu, ilevel - min_level] = f.tell() level_count_view[icpu, ilevel - min_level] = file_ncache - f.seek(skip_len(Nskip, file_ncache, size), SEEK_CUR) + f.seek(skip_len(Nskip, file_ncache, single), SEEK_CUR) return offset, level_count @@ -176,7 +176,7 @@ def fill_hydro(FortranFile f, INT64_t ndim, list all_fields, list fields, dict tr, RAMSESOctreeContainer oct_handler, - str size, + bint single, np.ndarray[np.int32_t, ndim=1] domain_inds=np.array([], dtype='int32'), ): cdef INT64_t offset @@ -219,7 +219,7 @@ def fill_hydro(FortranFile f, cdef int first_field_index = jumps[0] # The buffer is different depending on the size of the data - if size == "single": + if single: buffer_single = np.empty((level_count.max(), twotondim, nfields_selected), dtype="float32", order='F') else: buffer_double = np.empty((level_count.max(), twotondim, nfields_selected), dtype="float64", order='F') @@ -242,7 +242,7 @@ def fill_hydro(FortranFile f, offset = offsets[icpu, ilevel] if offset == -1: continue - f.seek(offset + skip_len(first_field_index, nc, size), SEEK_SET) + f.seek(offset + skip_len(first_field_index, nc, single), SEEK_SET) # We have already skipped the first fields (if any) # so we "rewind" (this will cancel the first seek) @@ -252,9 +252,9 @@ def fill_hydro(FortranFile f, for j in range(nfields_selected): jump_len += jumps[j] if jump_len > 0: - f.seek(skip_len(jump_len, nc, size), SEEK_CUR) + f.seek(skip_len(jump_len, nc, single), SEEK_CUR) jump_len = 0 - if size == "single": + if single: f.read_vector_inplace('f', &buffer_single[0, i, j]) else: f.read_vector_inplace('d', &buffer_double[0, i, j]) @@ -265,12 +265,12 @@ def fill_hydro(FortranFile f, # but since we're doing an absolute seek at the beginning of # the loop on CPUs, we can spare one seek here ## if jump_len > 0: - ## f.seek(skip_len(jump_len, nc, size), SEEK_CUR) + ## f.seek(skip_len(jump_len, nc, single), SEEK_CUR) # Alias buffer into dictionary tmp = {} for i, field in enumerate(fields): - if size == "single": + if single: tmp[field] = np.asarray(buffer_single[:, :, i], dtype="float64") else: tmp[field] = buffer_double[:, :, i] From fa416e14a90ca5e62fbc791eddb0a49764f87d3f Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 16 May 2025 12:16:40 +0100 Subject: [PATCH 3/6] read energies in rt groups of the info_rt_file_xxxxx.txt --- yt/frontends/ramses/field_handlers.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/yt/frontends/ramses/field_handlers.py b/yt/frontends/ramses/field_handlers.py index 587a45a0495..6b6c3398bcf 100644 --- a/yt/frontends/ramses/field_handlers.py +++ b/yt/frontends/ramses/field_handlers.py @@ -549,10 +549,13 @@ def detect_fields(cls, ds): rheader = {} - def read_rhs(cast): + def read_rhs(cast, group=None): line = f.readline() p, v = line.split("=") - rheader[p.strip()] = cast(v) + if group is not None: + rheader[f"Group {group} {p.strip()}"] = cast(v) + else: + rheader[p.strip()] = cast(v) with open(fname) as f: # Read nRTvar, nions, ngroups, iions @@ -574,7 +577,7 @@ def read_rhs(cast): read_rhs(float) f.readline() - # Reat unit_np, unit_pfd + # Read unit_np, unit_pf for _ in range(2): read_rhs(float) @@ -587,9 +590,20 @@ def read_rhs(cast): # Read n star, t2star, g_star for _ in range(3): read_rhs(float) - - # Touchy part, we have to read the photon group properties - mylog.debug("Not reading photon group properties") + + f.readline() + f.readline() + + # Get global group properties (groupL0, groupL1, spec2group) + for _ in range(3): + read_rhs(lambda line: [float(e) for e in line.split()]) + + # get egy for each group (to get proper energy densities) + for _ in range(4): + group = int(f.readline().split()[1]) + read_rhs(lambda line: [float(e) for e in line.split()], group) + f.readline() + f.readline() cls.rt_parameters[ds.unique_identifier] = rheader From b0c1817c46e889c00804887794eebc045c784b46 Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 16 May 2025 12:20:23 +0100 Subject: [PATCH 4/6] specifit that rt fields are single precision (temp megatron fix) --- yt/frontends/ramses/data_structures.py | 6 +++++- yt/frontends/ramses/io_utils.pyx | 6 +++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py index b8926931baa..a0923c3a163 100644 --- a/yt/frontends/ramses/data_structures.py +++ b/yt/frontends/ramses/data_structures.py @@ -455,7 +455,7 @@ def _fill_no_ghostzones(self, fd, fields, selector, file_handler): fields, data, oct_handler, - is_single # Hardcoded for now + is_single, ) return data @@ -478,6 +478,9 @@ def _fill_with_ghostzones( # Initializing data container for field in fields: tr[field] = np.zeros(cell_count, "float64") + + # Get the size of the data, single precision if RT field, double otherwise + is_single = True if any(field.startswith("Photon") for field in fields) else False # Do an early exit if the cell count is null if cell_count == 0: @@ -511,6 +514,7 @@ def _fill_with_ghostzones( fields, tr, oct_handler, + is_single, domain_inds=domain_inds, ) return tr diff --git a/yt/frontends/ramses/io_utils.pyx b/yt/frontends/ramses/io_utils.pyx index a740f62529c..600b18c84cb 100644 --- a/yt/frontends/ramses/io_utils.pyx +++ b/yt/frontends/ramses/io_utils.pyx @@ -22,7 +22,7 @@ cdef int DOUBLE_SIZE = sizeof(np.float64_t) cdef int SINGLE_SIZE = sizeof(np.float32_t) -cdef inline int skip_len(int Nskip, int record_len, bint single) noexcept nogil: +cdef inline int skip_len(int Nskip, int record_len, np.npy_bool single) noexcept nogil: # If the data is single precision, we need to skip 4 bytes if single: return Nskip * (record_len * SINGLE_SIZE + INT64_SIZE) @@ -115,7 +115,7 @@ def read_amr(FortranFile f, dict headers, @cython.wraparound(False) @cython.cdivision(True) @cython.nonecheck(False) -cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip, bint single): +cpdef read_offset(FortranFile f, INT64_t min_level, INT64_t domain_id, INT64_t nvar, dict headers, int Nskip, np.npy_bool single): cdef np.ndarray[np.int64_t, ndim=2] offset, level_count cdef INT64_t ndim, twotondim, nlevelmax, n_levels, nboundary, ncpu, ncpu_and_bound @@ -176,7 +176,7 @@ def fill_hydro(FortranFile f, INT64_t ndim, list all_fields, list fields, dict tr, RAMSESOctreeContainer oct_handler, - bint single, + np.npy_bool single, np.ndarray[np.int32_t, ndim=1] domain_inds=np.array([], dtype='int32'), ): cdef INT64_t offset From 1208d4963f875695faee2cadf0c4caa7bbeec5cf Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Fri, 16 May 2025 16:57:10 +0100 Subject: [PATCH 5/6] there are 8 energy groups in MEGATRON, not 4. --- yt/frontends/ramses/field_handlers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/frontends/ramses/field_handlers.py b/yt/frontends/ramses/field_handlers.py index 6b6c3398bcf..e54fff7a0e7 100644 --- a/yt/frontends/ramses/field_handlers.py +++ b/yt/frontends/ramses/field_handlers.py @@ -599,7 +599,7 @@ def read_rhs(cast, group=None): read_rhs(lambda line: [float(e) for e in line.split()]) # get egy for each group (to get proper energy densities) - for _ in range(4): + for _ in range(8): group = int(f.readline().split()[1]) read_rhs(lambda line: [float(e) for e in line.split()], group) f.readline() From 9599d96e4a08b3671af6b4b5e19b8a997390cb2f Mon Sep 17 00:00:00 2001 From: Anatole Storck Date: Tue, 27 May 2025 17:09:46 +0100 Subject: [PATCH 6/6] add note: rt field paramater readin only works for up to 8 energy bands --- yt/frontends/ramses/field_handlers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yt/frontends/ramses/field_handlers.py b/yt/frontends/ramses/field_handlers.py index e54fff7a0e7..f53fa0ded3c 100644 --- a/yt/frontends/ramses/field_handlers.py +++ b/yt/frontends/ramses/field_handlers.py @@ -590,7 +590,6 @@ def read_rhs(cast, group=None): # Read n star, t2star, g_star for _ in range(3): read_rhs(float) - f.readline() f.readline() @@ -598,7 +597,8 @@ def read_rhs(cast, group=None): for _ in range(3): read_rhs(lambda line: [float(e) for e in line.split()]) - # get egy for each group (to get proper energy densities) + # get egy for each group (used to get proper energy densities) + # NOTE: this assumes there are 8 energy groups for _ in range(8): group = int(f.readline().split()[1]) read_rhs(lambda line: [float(e) for e in line.split()], group)