|
| 1 | +//------------------------------- -*- C++ -*- -------------------------------// |
| 2 | +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details |
| 3 | +// SPDX-License-Identifier: (Apache-2.0 OR MIT) |
| 4 | +//---------------------------------------------------------------------------// |
| 5 | +//! \file celeritas/field/LoadCovfieField.covfie.cc |
| 6 | +//---------------------------------------------------------------------------// |
| 7 | +#include "LoadCovfieField.hh" |
| 8 | + |
| 9 | +#include <cstddef> |
| 10 | +#include <fstream> |
| 11 | +#include <string> |
| 12 | +#include <covfie/core/backend/primitive/array.hpp> |
| 13 | +#include <covfie/core/backend/transformer/affine.hpp> |
| 14 | +#include <covfie/core/backend/transformer/linear.hpp> |
| 15 | +#include <covfie/core/backend/transformer/strided.hpp> |
| 16 | +#include <covfie/core/field.hpp> |
| 17 | +#include <covfie/core/vector.hpp> |
| 18 | + |
| 19 | +#include "corecel/Assert.hh" |
| 20 | +#include "corecel/cont/Range.hh" |
| 21 | +#include "corecel/io/Logger.hh" |
| 22 | +#include "geocel/Types.hh" |
| 23 | +#include "celeritas/Types.hh" |
| 24 | + |
| 25 | +namespace celeritas |
| 26 | +{ |
| 27 | +namespace |
| 28 | +{ |
| 29 | +//---------------------------------------------------------------------------// |
| 30 | +// Stateless covfie interpolation backend used solely for deserialization. |
| 31 | +// Any stateless backend (linear, nearest_neighbour, etc.) produces an |
| 32 | +// identical binary format because stateless backends write no IO headers. |
| 33 | +// We use linear as an arbitrary choice; the file contents are the same |
| 34 | +// regardless. |
| 35 | +template<class B> |
| 36 | +using deserialization_interp_t = covfie::backend::linear<B>; |
| 37 | + |
| 38 | +// Covfie 3D field type for Cartesian maps. |
| 39 | +using storage3_t = covfie::backend::array<covfie::vector::float3>; |
| 40 | +using strided3_t = covfie::backend::strided<covfie::vector::size3, storage3_t>; |
| 41 | +using file_cart_t |
| 42 | + = covfie::field<covfie::backend::affine<deserialization_interp_t<strided3_t>>>; |
| 43 | + |
| 44 | +// Covfie 2D field type for BrBz cylindrical maps. |
| 45 | +using storage2_t = covfie::backend::array<covfie::vector::float2>; |
| 46 | +using strided2_t = covfie::backend::strided<covfie::vector::size2, storage2_t>; |
| 47 | +using file_rz_t |
| 48 | + = covfie::field<covfie::backend::affine<deserialization_interp_t<strided2_t>>>; |
| 49 | + |
| 50 | +//---------------------------------------------------------------------------// |
| 51 | +/*! |
| 52 | + * Recover world-coordinate grid from a covfie affine axis. |
| 53 | + * |
| 54 | + * The affine transform maps world coordinates to grid indices: |
| 55 | + * grid_index = scale * world_pos + translation |
| 56 | + * so world_pos_min = -t/s and world_pos_max = (n - 1 - t) / s. |
| 57 | + */ |
| 58 | +AxisGrid<real_type> from_affine(float scale, float translation, std::size_t n) |
| 59 | +{ |
| 60 | + CELER_VALIDATE(scale > 0, |
| 61 | + << "covfie affine transform has non-positive scale factor " |
| 62 | + << scale); |
| 63 | + CELER_VALIDATE(n > 1, |
| 64 | + << "covfie field grid axis is degenerate with " << n |
| 65 | + << " points"); |
| 66 | + |
| 67 | + AxisGrid<real_type> result; |
| 68 | + result.min = -translation / scale; |
| 69 | + result.max = (static_cast<float>(n) - 1.f - translation) / scale; |
| 70 | + result.num = static_cast<size_type>(n); |
| 71 | + return result; |
| 72 | +} |
| 73 | + |
| 74 | +//---------------------------------------------------------------------------// |
| 75 | +} // namespace |
| 76 | + |
| 77 | +//---------------------------------------------------------------------------// |
| 78 | +/*! |
| 79 | + * Load a Cartesian magnetic field map from a binary covfie file. |
| 80 | + */ |
| 81 | +CartMapFieldInput load_covfie_cart_field(std::string const& filename) |
| 82 | +{ |
| 83 | + std::ifstream ifs(filename, std::ifstream::binary); |
| 84 | + CELER_VALIDATE(ifs.good(), |
| 85 | + << "failed to open covfie field file '" << filename << "'"); |
| 86 | + |
| 87 | + file_cart_t file_field(ifs); |
| 88 | + ifs.close(); |
| 89 | + |
| 90 | + // Access the affine backend for the transform matrix |
| 91 | + auto const& affine_data = file_field.backend(); |
| 92 | + auto const& mat = affine_data.get_configuration(); |
| 93 | + |
| 94 | + // Access the strided backend for grid dimensions [nx, ny, nz] |
| 95 | + auto const& strided_data = affine_data.get_backend().get_backend(); |
| 96 | + auto const sizes = strided_data.get_configuration(); |
| 97 | + |
| 98 | + // Recover world-coordinate grid from affine transform |
| 99 | + CartMapFieldInput inp; |
| 100 | + inp.x = from_affine(mat(0, 0), mat(0, 3), sizes[0]); |
| 101 | + inp.y = from_affine(mat(1, 1), mat(1, 3), sizes[1]); |
| 102 | + inp.z = from_affine(mat(2, 2), mat(2, 3), sizes[2]); |
| 103 | + |
| 104 | + // Allocate field data: layout is [X][Y][Z][3] |
| 105 | + std::size_t const nx = inp.x.num; |
| 106 | + std::size_t const ny = inp.y.num; |
| 107 | + std::size_t const nz = inp.z.num; |
| 108 | + std::size_t const num_components = static_cast<std::size_t>(Axis::size_); |
| 109 | + inp.field.resize(num_components * nx * ny * nz); |
| 110 | + |
| 111 | + // Read grid node values directly from the strided backend |
| 112 | + file_cart_t::view_t field_view{file_field}; |
| 113 | + auto const& strided_view = field_view.backend().get_backend().get_backend(); |
| 114 | + |
| 115 | + for (auto ix : range(nx)) |
| 116 | + { |
| 117 | + for (auto iy : range(ny)) |
| 118 | + { |
| 119 | + for (auto iz : range(nz)) |
| 120 | + { |
| 121 | + auto const bvec = strided_view.at({ix, iy, iz}); |
| 122 | + |
| 123 | + auto const base = ((ix * ny + iy) * nz + iz) * num_components; |
| 124 | + |
| 125 | + inp.field[base + 0] = bvec[0]; |
| 126 | + inp.field[base + 1] = bvec[1]; |
| 127 | + inp.field[base + 2] = bvec[2]; |
| 128 | + } |
| 129 | + } |
| 130 | + } |
| 131 | + |
| 132 | + CELER_LOG(debug) << "Loaded covfie Cartesian field: " << nx << "x" << ny |
| 133 | + << "x" << nz << " grid, x=[" << inp.x.min << ", " |
| 134 | + << inp.x.max << "], y=[" << inp.y.min << ", " << inp.y.max |
| 135 | + << "], z=[" << inp.z.min << ", " << inp.z.max << "]"; |
| 136 | + CELER_ENSURE(inp); |
| 137 | + return inp; |
| 138 | +} |
| 139 | + |
| 140 | +//---------------------------------------------------------------------------// |
| 141 | +/*! |
| 142 | + * Load a 2D BrBz cylindrical magnetic field map from a binary covfie file. |
| 143 | + */ |
| 144 | +RZMapFieldInput load_covfie_rz_field(std::string const& filename) |
| 145 | +{ |
| 146 | + std::ifstream ifs(filename, std::ifstream::binary); |
| 147 | + CELER_VALIDATE(ifs.good(), |
| 148 | + << "failed to open covfie field file '" << filename << "'"); |
| 149 | + |
| 150 | + file_rz_t file_field(ifs); |
| 151 | + ifs.close(); |
| 152 | + |
| 153 | + // Access the affine backend for the [2x3] transform matrix |
| 154 | + auto const& affine_data = file_field.backend(); |
| 155 | + auto const& mat = affine_data.get_configuration(); |
| 156 | + |
| 157 | + // Access the strided backend for grid dimensions [nr, nz] |
| 158 | + auto const& strided_data = affine_data.get_backend().get_backend(); |
| 159 | + auto const sizes = strided_data.get_configuration(); |
| 160 | + |
| 161 | + // Recover world-coordinate grids from affine transform |
| 162 | + // Note: 2D affine has translation in column 2, not 3 |
| 163 | + auto grid_r = from_affine(mat(0, 0), mat(0, 2), sizes[0]); |
| 164 | + auto grid_z = from_affine(mat(1, 1), mat(1, 2), sizes[1]); |
| 165 | + |
| 166 | + RZMapFieldInput inp; |
| 167 | + inp.num_grid_r = grid_r.num; |
| 168 | + inp.num_grid_z = grid_z.num; |
| 169 | + inp.min_r = grid_r.min; |
| 170 | + inp.max_r = grid_r.max; |
| 171 | + inp.min_z = grid_z.min; |
| 172 | + inp.max_z = grid_z.max; |
| 173 | + |
| 174 | + // Allocate field data: layout is [Z][R] (R has stride 1) |
| 175 | + std::size_t const nr = grid_r.num; |
| 176 | + std::size_t const nz = grid_z.num; |
| 177 | + std::size_t const grid_size = nr * nz; |
| 178 | + inp.field_r.resize(grid_size); |
| 179 | + inp.field_z.resize(grid_size); |
| 180 | + |
| 181 | + // Read grid node values directly from the strided backend |
| 182 | + file_rz_t::view_t field_view{file_field}; |
| 183 | + auto const& strided_view = field_view.backend().get_backend().get_backend(); |
| 184 | + |
| 185 | + for (auto iz : range(nz)) |
| 186 | + { |
| 187 | + for (auto ir : range(nr)) |
| 188 | + { |
| 189 | + auto const bvec = strided_view.at({ir, iz}); |
| 190 | + |
| 191 | + // Index: [Z][R] with R stride 1 |
| 192 | + std::size_t const idx = iz * nr + ir; |
| 193 | + |
| 194 | + inp.field_r[idx] = bvec[0]; |
| 195 | + inp.field_z[idx] = bvec[1]; |
| 196 | + } |
| 197 | + } |
| 198 | + |
| 199 | + CELER_LOG(debug) << "Loaded covfie BrBz field: " << nr << "x" << nz |
| 200 | + << " grid, r=[" << grid_r.min << ", " << grid_r.max |
| 201 | + << "], z=[" << grid_z.min << ", " << grid_z.max << "]"; |
| 202 | + CELER_ENSURE(inp); |
| 203 | + return inp; |
| 204 | +} |
| 205 | + |
| 206 | +//---------------------------------------------------------------------------// |
| 207 | +} // namespace celeritas |
0 commit comments