corecel/grid/TwodSubgridCalculator.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file corecel/grid/TwodSubgridCalculator.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/data/Collection.hh"
0011
0012 #include "FindInterp.hh"
0013 #include "NonuniformGrid.hh"
0014 #include "TwodGridData.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Do bilinear interpolation on a 2D grid with the x value preselected.
0021  *
0022  * This is usually not called directly but rather given as the return result of
0023  * the TwodGridCalculator.
0024  */
0025 class TwodSubgridCalculator
0026 {
0027   public:
0028     //!@{
0029     //! \name Type aliases
0030     using Values
0031         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0032     using InterpT = FindInterp<real_type>;
0033     //!@}
0034
0035   public:
0036     // Construct with grid data, backend values, and lower X data.
0037     inline CELER_FUNCTION TwodSubgridCalculator(TwodGridData const& grid,
0038                                                 Values const& storage,
0039                                                 InterpT x_loc);
0040
0041     // Calculate the value at the given y coordinate
0042     inline CELER_FUNCTION real_type operator()(real_type y) const;
0043
0044     //! Index of the preselected lower x value
0045     CELER_FORCEINLINE_FUNCTION size_type x_index() const
0046     {
0047         return x_loc_.index;
0048     }
0049
0050     //! Fraction between the lower and upper x grid values
0051     CELER_FORCEINLINE_FUNCTION real_type x_fraction() const
0052     {
0053         return x_loc_.fraction;
0054     }
0055
0056   private:
0057     TwodGridData const& grids_;
0058     Values const& storage_;
0059     InterpT const x_loc_;
0060
0061     inline CELER_FUNCTION real_type at(size_type x_idx, size_type y_idx) const;
0062 };
0063
0064 //---------------------------------------------------------------------------//
0065 // INLINE DEFINITIONS
0066 //---------------------------------------------------------------------------//
0067 /*!
0068  * Construct with grid data, backend values, and lower X data.
0069  *
0070  * This is typically constructed from a TwodGridCalculator. The interpolated x
0071  * location could be extended to allow a fractional value of 1 to support
0072  * interpolating on the highest value of the x grid.
0073  */
0074 CELER_FUNCTION
0075 TwodSubgridCalculator::TwodSubgridCalculator(TwodGridData const& grids,
0076                                              Values const& storage,
0077                                              InterpT x_loc)
0078     : grids_{grids}, storage_(storage), x_loc_(x_loc)
0079 {
0080     CELER_EXPECT(grids);
0081     CELER_EXPECT(grids.values.back() < storage.size());
0082     CELER_EXPECT(x_loc.index + 1 < grids.x.size());
0083     CELER_EXPECT(x_loc.fraction >= 0 && x_loc_.fraction < 1);
0084 }
0085
0086 //---------------------------------------------------------------------------//
0087 /*!
0088  * Calculate the value at the given y coordinate for preselected x.
0089  *
0090  * This uses *bilinear* interpolation and and therefore exactly represents
0091  * functions that are a linear combination of 1, x, y, and xy.
0092  */
0093 CELER_FUNCTION real_type TwodSubgridCalculator::operator()(real_type y) const
0094 {
0095     NonuniformGrid<real_type> const y_grid{grids_.y, storage_};
0096     CELER_EXPECT(y >= y_grid.front() && y < y_grid.back());
0097
0098     InterpT const y_loc = find_interp(y_grid, y);
0099     auto at_corner = [this, y_loc](size_type xo, size_type yo) -> real_type {
0100         return this->at(x_loc_.index + xo, y_loc.index + yo);
0101     };
0102
0103     return (1 - x_loc_.fraction)
0104                * ((1 - y_loc.fraction) * at_corner(0, 0)
0105                   + (y_loc.fraction) * at_corner(0, 1))
0106            + (x_loc_.fraction)
0107                  * ((1 - y_loc.fraction) * at_corner(1, 0)
0108                     + (y_loc.fraction) * at_corner(1, 1));
0109 }
0110
0111 //---------------------------------------------------------------------------//
0112 /*!
0113  * Get the value at the specified x/y coordinate.
0114  *
0115  * NOTE: this must match TwodGridData::index.
0116  */
0117 CELER_FUNCTION real_type TwodSubgridCalculator::at(size_type x_idx,
0118                                                    size_type y_idx) const
0119 {
0120     return storage_[grids_.at(x_idx, y_idx)];
0121 }
0122
0123 //---------------------------------------------------------------------------//
0124 }  // namespace celeritas