From de32f09b6b2bf950ea6315ebc8ee5f72fd16df52 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 5 May 2021 10:54:00 -0400 Subject: [PATCH] Updates to EnzoRiemann related to the shift to using primitives to compute the fluxes (and the new "integration quantity" terminology) Other changes include: - Adjusted the max asymmetry residual in the Dual Energy Cloud Test for the HLLC Riemann Solver. - I transitioned from registering the integration quantities with the Riemann Solver (in the Riemann Solver factory method) to having the Riemann Solver directly specify its expected primitives and integrations quantities. This seems to make much more sense than independently coming up with lists of integration and primitive quanities in EnzoMethodMHDVlct and making sure that they are consistent with the Riemann Solver's requirements. This change includes tweaks to EnzoRiemann::construct_riemann and the introduction of the EnzoRiemann::integration_quantities() and EnzoRiemann::primitive_quantities() interface methods. - This change made the EnzoMethodMHDVlct::determine_quantities_ helper method unnecessary, so I deleted it. - I also deleted the now unnecessary method from EnzoEquationOfState, eint_from_primitive. Finally, I updated the website documentation for EnzoRiemann to reflect these changes --- doc/source/devel/hydro_infrastructure.rst | 251 +++++++++++----------- input/vlct/run_dual_energy_cloud_test.py | 2 +- src/Enzo/enzo_EnzoEOSIdeal.cpp | 28 --- src/Enzo/enzo_EnzoEOSIdeal.hpp | 3 - src/Enzo/enzo_EnzoEquationOfState.hpp | 15 -- src/Enzo/enzo_EnzoMethodMHDVlct.cpp | 68 ++---- src/Enzo/enzo_EnzoMethodMHDVlct.hpp | 17 -- src/Enzo/enzo_EnzoRiemann.cpp | 38 ++-- src/Enzo/enzo_EnzoRiemann.hpp | 37 ++-- src/Enzo/enzo_EnzoRiemannHLL.hpp | 6 +- src/Enzo/enzo_EnzoRiemannImpl.hpp | 192 ++++++++--------- src/Enzo/enzo_EnzoRiemannLUT.hpp | 118 +++++----- src/Enzo/enzo_EnzoRiemannUtils.hpp | 47 ++-- 13 files changed, 388 insertions(+), 434 deletions(-) diff --git a/doc/source/devel/hydro_infrastructure.rst b/doc/source/devel/hydro_infrastructure.rst index 055fa506a0..eee43e060d 100644 --- a/doc/source/devel/hydro_infrastructure.rst +++ b/doc/source/devel/hydro_infrastructure.rst @@ -203,27 +203,24 @@ classes include: Each of these operation classes are fairly modular (to allow for selective usage of the frame work components). However, all of the classes require that an instance of ``EnzoEquationOfState`` get's -passed. +passed. The operation classes are also provided with ``PUP`` methods +to allow for easy serialization alongside the ``Method`` class that +makes use of them. Each of the operation classes are designed to be configured upon initialization. The instances can then be used multiple times per time-step (along multiple dimensions if the operation is directional) -and in other time-steps. The operation classes are also provided with -``PUP`` methods to allow for easy migration alongside the ``Method`` -class that makes use them. - -For each operation class (other than ``EnzoEquationOfState``), the -expected integration quantities or primitives (other than passively -advected scalars) are *registered* at construction. - - * The names of all integration quantites that get registered - in the construction of ``EnzoRiemann`` must share a name - with the registered quantities in ``FIELD_TABLE``. - - * All registered integration quantity names in the construction of - ``EnzoRiemann`` or ``EnzoIntegrationQuanUpdate`` must be specified - in ``FIELD_TABLE`` as quantities that are actively advected in - some contexts. +and in other time-steps. Lists (excluding passive scalars) of the +expected primitives and integration quantities are respectively +*registered* during the construction of ``EnzoReconstructor`` and +``EnzoIntegrationQuanUpdate``. These quantities must each share a name +with the registered quantities in ``FIELD_TABLE``. In contrast, +configuration of ``EnzoRiemann``, is less flexible and instances +actually specify the non-passive integration quantities and +non-passive primitives that they require. This difference exists +because the operations encapsulated by ``EnzoReconstructor`` and +``EnzoIntegrationQuanUpdate`` can be applied to individual quantities +in far a more independent manner. Because all fields storing passively advected scalars are not necessarily known when initializing a hydro/MHD integrator (i.e. @@ -589,38 +586,47 @@ static factory method: .. code-block:: c++ - EnzoRiemann* EnzoRiemann::construct_riemann - (std::vector integrable_quantities, - std::string solver); + EnzoRiemann* EnzoRiemann::construct_riemann(std::string solver, bool mhd, + bool internal_energy); + +The factory method requires that we specify the name of the solver (via +``solver``), whether magnetic fields are present (via ``mhd``), and whether +the internal energy flux must be computed (via ``internal_energy``). + +An instance of ``EnzoRiemann`` specifies the required non-passive +integration quantities with: + +.. code-block:: c++ + + const std::vector integration_quantities() const; + +and the required non-passive primitives with: + +.. code-block:: c++ + + const std::vector primitive_quantites() const; -The factory method requires that we both register the names of the -integrable quantities (excluding passively advected scalars), with -``integrable_quantities``, and specify the name of the solver -``solver``. Note that the names of the integrable quantites should -match quantities specified in ``FIELD_TABLE`` that are identified as -being actively advected. For more details about ``FIELD_TABLE``, see -:ref:`Centered-Field-Registry` The main interface function of ``EnzoRiemann`` is: .. code-block:: c++ - void solve - (EnzoEFltArrayMap &prim_map_l, EnzoEFltArrayMap &prim_map_r, - const EFlt3DArray &pressure_array_l, const EFlt3DArray &pressure_array_r, - EnzoEFltArrayMap &flux_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const std::vector &passive_lists, - EFlt3DArray *interface_velocity) + void solve(EnzoEFltArrayMap &prim_map_l, EnzoEFltArrayMap &prim_map_r, + EnzoEFltArrayMap &flux_map, int dim, EnzoEquationOfState *eos, + int stale_depth, const str_vec_t &passive_list, + EFlt3DArray *interface_velocity) const; In this function, the ``prim_map_l`` and ``prim_map_r`` arguments are references to the ``EnzoEFltArrayMap`` objects holding the arrays of -reconstructed left/right integrable quantities and passively advected -scalars. The ``pressure_array_l``/ ``pressure_array_r`` arguments -specify arrays holding the left/right reconstructed pressure. The -``flux_map`` argument holds the face-centered arrays where the -computed fluxes for each integrable quantity and passively advected -scalar will be stored. ``dim`` indicates the dimension along which the -flux should be computed (0,1,2 corresponds to x,y,z). +reconstructed left/right primitive quantities. These maps should have +entries for each primitive listed by +``EnzoRiemann::primitive_quantities()`` and passive scalar in +``passive_list``. The ``flux_map`` argument holds the face-centered +arrays where the computed fluxes for each integration quantity are +written. This should have entries for each quantity listed in +``EnzoRiemann::integration_quantities()`` and passive in +``passive_list``. ``dim`` indicates the dimension along which the flux +should be computed (0,1,2 corresponds to x,y,z). ``interface_velocity`` is an optional argument used to specify a pointer to an array that can be used to store interface velocity values computed by the Riemann Solver (this is primarily used for @@ -646,11 +652,14 @@ solver-specific calculations and is called at every cell-interface. Additionally, the functor also specifies a specialization of the template class ``EnzoRiemannLUT`` that primarily - * Specifies the exact set of actively advected integrable quantities - that a given solver expects + * Specifies the exact set of actively advected integration quantities + and primitive quantities that a given solver expects. Technically, + the primitives and any optional active integration quantities, like + ``"internal_energy"``, are not directly specified by the lookup table, + but ``EnzoRiemannImpl`` accounts for this. * Serves as a compile-time lookup table. It statically maps the names of the all of the components of the relevant actively advected - quantities to unique array indices. + integration quantities to unique array indices. See :ref:`EnzoRiemannLUT-section` for a more detailed description of ``EnzoRiemannLUT`` and @@ -670,26 +679,25 @@ computed. At each location, the following sequence of operations are performed: 1. Retrieve the left and right primitives at the given location from - the input arrays and stores them in stack-allocated arrays of - ``enzo_float`` elements called ``wl`` and ``wr``. As mentioned - above, the values are organized according to the specialization - of ``EnzoRiemannLUT`` provided by the ``ImplFunctor`` - (hereafter, ``ImplFunctor::LUT``) - 2. The left and right pressure values are retrieved from the - temporary fields holding the values that were precomputed from - the reconstructed quantities (presumably using a concrete - subclass of ``EnzoEquationOfState``). The values are stored in + the input arrays and stores them in stack-allocated ``enzo_float`` + arrays called ``wl`` and ``wr``. As mentioned above, the values are + organized according to the specialization of + ``EnzoRiemannLUT`` provided by the ``ImplFunctor`` + (hereafter, ``ImplFunctor::LUT``). *Note: for non-barotropic + equations of state* ``pressure`` *is stored at* + ``ImplFunctor::LUT::total_energy``. + 2. The left and right pressure values are determined (they may have + been precomputed using a concrete subclass of + ``EnzoEquationOfState``). The values are stored in local variables ``pressure_l`` and ``pressure_r``. 3. The conserved forms of the left and right reconstructed - primitives and stored in the arrays called ``Ul`` and - ``Ur``. Primitives that are always in conserved form (e.g. - density or magnetic field). The elements of ``Ul`` / ``Ur`` - are also ordered by ``ImplFunctor::LUT`` (e.g. the index for a - given component of the velocity in ``wl`` / ``wr`` matches the - index for the same component of the momentum in ``Ul`` / ``Ur``). - 4. The standard left and right hydro/MHD fluxes are computed using - the above quantities and stored in ``Fl`` and ``Fr`` (organized by - ``ImplFunctor::LUT``) + integration quantities are computed and stored in the arrays + called ``Ul`` and ``Ur`` (organized by ``ImplFunctor::LUT``) + *Note: There may be some duplication of values between* + ``Ul`` *&* ``Ur`` *and* ``wl`` *&* ``Ur``. + 4. The standard left and right integration quantity fluxes fluxes are + computed using the above quantities and stored in ``Fl`` and ``Fr`` + (organized by ``ImplFunctor::LUT``) 5. These quantities are all passed to the static public ``operator()`` method provided by ``ImplFunctor`` that returns the array of interface fluxes in the array, ``fluxes``. (It also @@ -697,11 +705,9 @@ performed: 6. The interface fluxes and interface velocity are then copied into the output fields. -A separate method is provided to compute the fluxes for the passively -advected quantities. This method will also be compute the fluxes of any -specified quantities that are nominally actively advected, but can fall -back to using passive advection when the solver doesn't explictly support -it (the main example is ``"internal_energy"``) +After computing the fluxes for all of the actively advected integration +quantities at all locations, a helper method is invoked to compute the +fluxes for the passively advected quantities. *Note: Currently EnzoRiemannImpl has only been tested and known to work for 3D problems. Additionally, no solvers (or more specifically, @@ -728,9 +734,9 @@ The class is expected to: * publically define the ``LUT`` type, which should be a specialization of the ``EnzoRiemannLUT`` template class. ``ImplFunctor::LUT`` should indicate which actively advected - quantities are expected by ``ImplFunctor`` and how they organized. - For more details about how how ``EnzoRiemannLUT`` is used, - see :ref:`EnzoRiemannLUT-section` + integration quantities are expected by ``ImplFunctor`` and how they + are organized. For more details about how ``EnzoRiemannLUT`` + is used, see :ref:`EnzoRiemannLUT-section` * provide the const-qualified function call method, ``operator()``. @@ -753,14 +759,17 @@ holding the Riemann Flux at a given cell-interface. Note that ``lutarray`` is actually an alias for ``std::array``. Each of these arrays hold values associated with the components of each relevant -actively advected quantity and are organized according to -``ImplFunctor::LUT`` (again, see :ref:`EnzoRiemannLUT-section` for -more details about the ``LUT`` type). - -``flux_l``/ ``flux_r``, ``prim_l``/ ``prim_r``, and ``cons_l``/ -``cons_r`` store the left/right interface fluxes values, primitive -quantities, and conserved quantities (they are passed ``Fl``/ ``Fr``, -``wl``/ ``wr``, and ``Ul``/ ``Ur``, respectively). +actively advected integration/primitive quantity and are organized +according to ``ImplFunctor::LUT`` (again, see +:ref:`EnzoRiemannLUT-section` for more details about the ``LUT`` type). + +``flux_l``/ ``flux_r`` and ``cons_l``/ ``cons_r`` store the left/right +interface fluxes values and conserved quantities (they are passed +respectively passed ``Fl``/ ``Fr`` and ``Ul``/ ``Ur``, respectively). +``prim_l``/ ``prim_r`` store the left/right interface primitive +values, and are passed ``wl``/ ``wr``. As mentioned before, for +non-barotropic equations of state, ``prim_l``/ ``prim_r`` store +pressure at ``ImplFunctor::LUT::total_energy`` The left and right reconstructed pressure values are passed as ``pressure_l`` and ``pressure_r``. ``barotropic_eos`` indicates @@ -792,45 +801,45 @@ This is a template class that provides the following features at compile time: * a lookup table (LUT) that maps the names of components of a subset - of the actively advected quantities defined in ``FIELD_TABLE`` to - unique, contiguous indices. + of the actively advected integration quantities defined in + ``FIELD_TABLE`` to unique, contiguous indices. - * the number of quantity components included in the table + * the number of integration quantity components included in the table - * a way to iterate over just the conserved quantities or specific + * a way to iterate over just the conserved or specific integration quantities values that are stored in an array using these mapping - * a way to query which of the actively advected quantities in - FIELD_TABLE are not included in the LUT + * a way to query which of the actively advected integration quantities + in FIELD_TABLE are not included in the LUT These feature are provided via the definition of publicly accessible integer constants in every specialization of the template class. All specializations have: - * a constant called ``NEQ`` equal to the number of quantity components - included in the lookup table + * a constant called ``NEQ`` equal to the number of integration quantity + components included in the lookup table * a constant called ``specific_start`` equal to the number of components - of conserved quantities included in the lookup table + of conserved integration quantities included in the lookup table * ``qkey`` constants, which include constants named for the components - of ALL actively advected quantities in FIELD_TABLE. A constant - associated with a SCALAR quantity, ``{qname}``, is simply called - ``{qname}`` while constants associated with a vector quantity + of ALL actively advected integration quantities in FIELD_TABLE. A + constant associated with a SCALAR quantity, ``{qname}``, is simply + called ``{qname}`` while constants associated with a vector quantity ``{qname}`` are called ``{qname}_i``, ``{qname}_j``, and ``{qname}_k``. -The `qkey` constants serve as both the keys of the lookup table and a +The ``qkey`` constants serve as both the keys of the lookup table and a way to check whether a component of an actively advected quantity is included in the table. Their values are satisfy the following conditions: - * All constants named for values corresponding to quantities included - in the table have values of ``-1`` + * All constants named for values corresponding to quantities NOT + included in the lookup table have values of ``-1`` - * All constants named for conserved quantities have unique integer - values in the internal ``[0,specific_start)`` + * All constants named for conserved integration quantities have unique + integer values in the internal ``[0,specific_start)`` - * All constants named for specific quantities have unique integer - values in the interval ``[specific_start, NEQ)`` + * All constants named for specific integration quantities have unique + integer values in the interval ``[specific_start, NEQ)`` The lookup table is always expected to include density and the 3 velocity components. Although it may not be strictly enforced (yet), the lookup @@ -842,17 +851,17 @@ programmatically probe the table's contents at runtime and validate that the above requirements are specified. For the sake of providing some concrete examples about how the code works, -let's assume that we have a class ``MyInputLUT`` that is defined as: +let's assume that we have a class ``MyIntegLUT`` that is defined as: .. code-block:: c++ - struct MyIntLUT { + struct MyIntegLUT { enum vals { density=0, velocity_i, velocity_j, velocity_k, total_energy, NEQ, specific_start = 1}; }; -The template specialization ``EnzoRiemannLUT`` assumes that -all undefined `qkey` constants omitted from ``MyIntLUT`` are not included +The template specialization ``EnzoRiemannLUT`` assumes that +all undefined ``qkey`` constants omitted from ``MyIntegLUT`` are not included in the lookup table and will define them within the template specialization to have values of ``-1``. @@ -861,9 +870,11 @@ velocity, one would evaluate: .. code-block:: c++ - int density_index = EnzoRiemannLUT::density; //=0 - int vj_index = EnzoRiemannLUT::velocity_j; //=2 + int density_index = EnzoRiemannLUT::density; //=0 + int vj_index = EnzoRiemannLUT::velocity_j; //=2 +Additionally, the value of ``EnzoRiemannLUT::bfield_k`` would be +``-1``. It makes more sense to talk about the use of this template class when we have a companion array. For convenience, the alias template @@ -872,16 +883,16 @@ have a companion array. For convenience, the alias template ``std::array::NEQ>;``. As an example, imagine that the total kinetic energy density needs to be -computed at a single location from an values stored in an array, ``prim``, -of type ``lutarray>``: +computed at a single location from an values stored in an array, ``integ``, +of type ``lutarray>``: .. code-block:: c++ - using LUT = EnzoRiemannLUT; - enzo_float v2 = (prim[LUT::velocity_i] * prim[LUT::velocity_i] + - prim[LUT::velocity_j] * prim[LUT::velocity_j] + - prim[LUT::velocity_k] * prim[LUT::velocity_k]); - enzo_float kinetic = 0.5 * prim[LUT::density] * v2; + using LUT = EnzoRiemannLUT; + enzo_float v2 = (integ[LUT::velocity_i] * integ[LUT::velocity_i] + + integ[LUT::velocity_j] * integ[LUT::velocity_j] + + integ[LUT::velocity_k] * integ[LUT::velocity_k]); + enzo_float kinetic = 0.5 * integ[LUT::density] * v2; ``EnzoRiemannLUT``, makes it very easy to @@ -912,24 +923,14 @@ density at a single location for an arbitrary lookup table: Adding new quantites -------------------- -To add support for new actively advected integrable cell-centered +To add support for new actively advected integration cell-centered quantities (e.g. cosmic ray energy/flux), the table of cell-centered quantities (``FIELD_TABLE``) must be updated. See -:ref:`Centered-Field-Registry` -for more details. - -To add support for computing fluxes for such quantities, modifications -must be made to ``EnzoRiemannImpl``. Currently, an abstract base class -called for ``EnzoFluxFunctor`` is provided for this purpose. The idea -is define a subclass to be defined for each additional set of flux -calculations and then in then have the factory method, -``EnzoRiemann::construct_riemann``, pass an array of the relevant -functors to ``EnzoRiemannImpl``. - -*However, because the functors are called as pointers will probably -incur overhead. In reality, the better solution might be to hardcode -in the additonal flux calculation functions in some kind of helper -method of* ``EnzoRiemannImpl``. +:ref:`Centered-Field-Registry` for more details. To add support for +computing fluxes for such quantities, modifications must be made to +either ``EnzoRiemannImpl`` or the ``ImplFunctor`` of an existing +solver. Alternatively, for certain quantities, a brand new solver +may need to be introduced. Adding new solvers ------------------ diff --git a/input/vlct/run_dual_energy_cloud_test.py b/input/vlct/run_dual_energy_cloud_test.py index 395e6dbe0d..8ce4ba51a0 100644 --- a/input/vlct/run_dual_energy_cloud_test.py +++ b/input/vlct/run_dual_energy_cloud_test.py @@ -83,7 +83,7 @@ def analyze_tests(): r += check_cloud_asym('hlld_cloud_0.0625/hlld_cloud_0.0625.block_list', 'hlld_cloud', 5.5e-13) r += check_cloud_asym('hllc_cloud_0.0625/hllc_cloud_0.0625.block_list', - 'hllc_cloud', 3.6e-13) + 'hllc_cloud', 3.8e-13) r += check_cloud_asym('hlle_cloud_0.0625/hlle_cloud_0.0625.block_list', 'hlle_cloud', 3.e-13) n_passed = np.sum(r) diff --git a/src/Enzo/enzo_EnzoEOSIdeal.cpp b/src/Enzo/enzo_EnzoEOSIdeal.cpp index c5d08eb045..5b47247c44 100644 --- a/src/Enzo/enzo_EnzoEOSIdeal.cpp +++ b/src/Enzo/enzo_EnzoEOSIdeal.cpp @@ -186,34 +186,6 @@ void EnzoEOSIdeal::pressure_from_integration //---------------------------------------------------------------------- -void EnzoEOSIdeal::eint_from_primitive(EnzoEFltArrayMap &primitive, - EFlt3DArray &eint, - int stale_depth) const -{ - if (grackle_variable_gamma_()){ - // we would need to model a field with the - ERROR("EnzoEOSIdeal::eint_from_primitive", - "Doesn't currently support spatial variations in gamma"); - } - - EFlt3DArray density = primitive.get("density", stale_depth); - EFlt3DArray pressure = primitive.get("pressure", stale_depth); - - CSlice unstaled(stale_depth,-stale_depth); - EFlt3DArray e = eint.subarray(unstaled, unstaled, unstaled); - enzo_float inv_gm1 = 1./(get_gamma()-1.); - - for (int iz=0; izmethod_grackle_use_grackle){ + // we can remove the following once EnzoMethodGrackle no longer requires + // the internal_energy to be a permanent field + ASSERT("EnzoMethodMHDVlct::determine_quantities_", + ("Grackle cannot currently be used alongside this integrator " + "unless the dual-energy formalism is in use"), + eos_->uses_dual_energy_formalism()); + } +#endif /* CONFIG_USE_GRACKLE */ + // Determine whether magnetic fields are to be used mhd_choice_ = parse_bfield_choice_(mhd_choice); + riemann_solver_ = EnzoRiemann::construct_riemann + (rsolver, mhd_choice_ != bfield_choice::no_bfield, + eos_->uses_dual_energy_formalism()); + // determine integration and primitive quantities std::vector integration_quantities, primitive_quantities; - EnzoMethodMHDVlct::determine_quantities_(eos_, mhd_choice_, - integration_quantities, - primitive_quantities); + integration_quantities = riemann_solver_->integration_quantities(); + primitive_quantities = riemann_solver_->primitive_quantities(); // Initialize the remaining component objects half_dt_recon_ = EnzoReconstructor::construct_reconstructor (primitive_quantities, half_recon_name, (enzo_float)theta_limiter); full_dt_recon_ = EnzoReconstructor::construct_reconstructor (primitive_quantities, full_recon_name, (enzo_float)theta_limiter); - riemann_solver_ = EnzoRiemann::construct_riemann(integration_quantities, - rsolver); + integration_quan_updater_ = new EnzoIntegrationQuanUpdate(integration_quantities, true); @@ -131,51 +144,6 @@ EnzoMethodMHDVlct::bfield_choice EnzoMethodMHDVlct::parse_bfield_choice_ //---------------------------------------------------------------------- -void EnzoMethodMHDVlct::determine_quantities_ -(const EnzoEquationOfState *eos, EnzoMethodMHDVlct::bfield_choice mhd_choice, - std::vector &integration_quantities, - std::vector &primitive_quantities) noexcept -{ -#ifdef CONFIG_USE_GRACKLE - if (enzo::config()->method_grackle_use_grackle){ - // we can remove the following once EnzoMethodGrackle no longer requires - // the internal_energy to be a permanent field - ASSERT("EnzoMethodMHDVlct::determine_quantities_", - ("Grackle cannot currently be used alongside this integrator " - "unless the dual-energy formalism is in use"), - eos->uses_dual_energy_formalism()); - } -#endif /* CONFIG_USE_GRACKLE */ - - std::string common[] {"density", "velocity"}; - for (std::string quantity : common){ - integration_quantities.push_back(quantity); - primitive_quantities.push_back(quantity); - } - - if (mhd_choice != bfield_choice::no_bfield){ - integration_quantities.push_back("bfield"); - primitive_quantities.push_back("bfield"); - } - - if (eos->is_barotropic()){ - ERROR("EnzoMethodMHDVlct::determine_quantities_", - "Not presently equipped to handle barotropic equations of state."); - } else { - // add specific total energy to integration quantities - integration_quantities.push_back("total_energy"); - // add pressure to primitive quantities - primitive_quantities.push_back("pressure"); - // add specific internal energy to integration quantities (if using the - // dual energy formalism) - if (eos->uses_dual_energy_formalism()){ - integration_quantities.push_back("internal_energy"); - } - } -} - -//---------------------------------------------------------------------- - EnzoMethodMHDVlct::~EnzoMethodMHDVlct() { delete eos_; diff --git a/src/Enzo/enzo_EnzoMethodMHDVlct.hpp b/src/Enzo/enzo_EnzoMethodMHDVlct.hpp index a0e89e8a02..c918eba66a 100644 --- a/src/Enzo/enzo_EnzoMethodMHDVlct.hpp +++ b/src/Enzo/enzo_EnzoMethodMHDVlct.hpp @@ -145,23 +145,6 @@ class EnzoMethodMHDVlct : public Method { /// returns the bfield_choice enum that matches the input string bfield_choice parse_bfield_choice_(std::string choice) const noexcept; - /// Determines the integration quantities and primitives from (FIELD_TABLE) - /// that are to be integrated and reconstructed - /// - /// @param[in] eos Pointer to the fluid's EquationOfState object. - /// @param[in] mhd_choice Encodes how the integrator will handle B-fields - /// @param[out] integration_quantities Reference to a vector that get's filled - /// by this function with the integration quantities (matching names in - /// FIELD_TABLE) used by the integrator - /// @param[out] primitive_quantities Reference to a vector that get's - /// filled by this function with the names of primitive quantities - /// (matching names in FIELD_TABLE) that are used by the integrator for - /// reconstruction - static void determine_quantities_ - (const EnzoEquationOfState *eos, bfield_choice mhd_choice, - std::vector &integration_quantities, - std::vector &primitive_quantities) noexcept; - /// Checks that the mesh size is sufficiently large to handle the given ghost /// depth and confirms that the ghost depth is consistent with the /// requirements of the reconstructors diff --git a/src/Enzo/enzo_EnzoRiemann.cpp b/src/Enzo/enzo_EnzoRiemann.cpp index 1263f0b089..1eab721ef7 100644 --- a/src/Enzo/enzo_EnzoRiemann.cpp +++ b/src/Enzo/enzo_EnzoRiemann.cpp @@ -12,30 +12,44 @@ //---------------------------------------------------------------------- -EnzoRiemann* EnzoRiemann::construct_riemann -(std::vector integrable_quantities, std::string solver) +EnzoRiemann* EnzoRiemann::construct_riemann(std::string solver, bool mhd, + bool internal_energy) { // determine the type of solver to construct: // convert string to lower case (https://stackoverflow.com/a/313990) std::string formatted(solver.size(), ' '); std::transform(solver.begin(), solver.end(), formatted.begin(), ::tolower); - EnzoRiemann* out; + EnzoRiemann* out = nullptr; // set to NULL to suppress compiler warnings + + if (formatted == "hll"){ + ASSERT("EnzoRiemann::construct_riemann", + ("An \"HLL\" Riemann solver without magnetic fields isn't " + "currently supported."), mhd); + out = new EnzoRiemannHLLMHD(internal_energy); + + } else if (formatted == "hlle"){ + if (mhd){ + out = new EnzoRiemannHLLEMHD(internal_energy); + } else { + ERROR("EnzoRiemann::construct_riemann", + "The \"HLLE\" Riemann solver without magnetic fields is untested"); + out = new EnzoRiemannHLLE(internal_energy); + } - // Eventually we may want to check for non-MHD Riemann solvers - if (formatted == std::string("hll")){ - out = new EnzoRiemannHLLMHD(integrable_quantities); - } else if (formatted == std::string("hlle")){ - out = new EnzoRiemannHLLEMHD(integrable_quantities); } else if (formatted == std::string("hllc")){ - out = new EnzoRiemannHLLC(integrable_quantities); + ASSERT("EnzoRiemann::construct_riemann", + "The \"HLLC\" Riemann Solver can't support mhd", !mhd); + out = new EnzoRiemannHLLC(internal_energy); + } else if (formatted == std::string("hlld")){ - // could possibly check that MHD fields are included - out = new EnzoRiemannHLLD(integrable_quantities); + ASSERT("EnzoRiemann::construct_riemann", + "The \"HLLD\" Riemann Solver requires magnetic fields", mhd); + out = new EnzoRiemannHLLD(internal_energy); + } else { ERROR("EnzoRiemann::construct_riemann", "The only known solvers are HLL, HLLE, HLLC, & HLLD"); - out = NULL; // Deals with compiler warning } return out; diff --git a/src/Enzo/enzo_EnzoRiemann.hpp b/src/Enzo/enzo_EnzoRiemann.hpp index ffd3e3d14a..58351c6495 100644 --- a/src/Enzo/enzo_EnzoRiemann.hpp +++ b/src/Enzo/enzo_EnzoRiemann.hpp @@ -22,15 +22,15 @@ class EnzoRiemann : public PUP::able /// Factory method for constructing the EnzoRiemann object. (The signature /// may need to be modified as additional physics get added) /// - /// @param integrable_quantities A vector of integration quantities (which - /// must be listed as advected quantities in FIELD_TABLE). This is used - /// to register the quantities that the Riemann Solver operates on. /// @param solver The name of the Riemann solver to use. Valid names include /// "hll", "hlle", and "hlld" - static EnzoRiemann* construct_riemann - (std::vector integrable_quantities, std::string solver); + /// @param mhd Indicates whether magnetic fields are present + /// @param internal_energy Indicates whether the internal energy is an + /// integration quantity + static EnzoRiemann* construct_riemann(std::string solver, bool mhd, + bool internal_energy); - EnzoRiemann() throw() + EnzoRiemann() noexcept {} /// Virtual destructor @@ -54,12 +54,16 @@ class EnzoRiemann : public PUP::able /// Computes the Riemann Fluxes for each conserved field along a given /// dimension, dim /// @param[in] priml_map,primr_map Maps of arrays holding the left/right - /// reconstructed face-centered primitives. This should be - /// face-centered along `dim` (without having values on the exterior - /// faces of the block) and cell-centered along the other dimensions. - /// @param[out] flux_map Holds arrays where the calculated fluxes - /// will be stored. The arrays should be face-centered along `dim` - /// (without having values on the exterior faces of the block) + /// reconstructed face-centered primitives. A list of the expected + /// primitive quantities is provided by the `primitive_quantities` method. + /// These should be face-centered along `dim` (without having values on + /// the exterior faces of the block) and cell-centered along the other + /// dimensions. + /// @param[out] flux_map Holds arrays where the calculated fluxes for the + /// integration quantities will be stored. Fluxes are computed for each + /// integration quantities in the list provided by the + /// ``integration_quantities`` method. The arrays should be face-centered + /// along `dim` (without having values on the exterior faces of the block) /// @param[in] dim Dimension along which to compute Riemann fluxes. /// Values of 0, 1, and 2 correspond to the x, y, and z directions. /// @param[in] eos Instance of the fluid's EnzoEquationOfState object @@ -85,6 +89,15 @@ class EnzoRiemann : public PUP::able int stale_depth, const str_vec_t &passive_list, EFlt3DArray *interface_velocity) const = 0; + /// Return the names of the actively advected integration quantities for + /// which fluxes will be computed. + virtual const std::vector integration_quantities() + const noexcept = 0; + + /// Return the names of the (non passive scalar) primitives that are required + /// to compute the flux. + virtual const std::vector primitive_quantities() + const noexcept = 0; }; #endif /* ENZO_ENZO_RIEMANN_HPP */ diff --git a/src/Enzo/enzo_EnzoRiemannHLL.hpp b/src/Enzo/enzo_EnzoRiemannHLL.hpp index 89a9a2f100..e334a3cd09 100644 --- a/src/Enzo/enzo_EnzoRiemannHLL.hpp +++ b/src/Enzo/enzo_EnzoRiemannHLL.hpp @@ -168,7 +168,7 @@ struct EinfeldtWavespeed template struct DavisWavespeed { - /// @class DavisHLLEWavespeed + /// @class DavisWavespeed /// @ingroup Enzo /// @brief [\ref Enzo] Encapsulates one of the wavespeed estimators /// described by Davis 1998. It was used in Enzo's Runge-Kutta MHD @@ -214,9 +214,9 @@ struct DavisWavespeed template struct HLLImpl { - /// @class EnzoRiemannHLLE + /// @class HLLImpl /// @ingroup Enzo - /// @brief [\ref Enzo] Encapsulates operations of HLLE approximate Riemann + /// @brief [\ref Enzo] Encapsulates operations of HLL approximate Riemann /// Solver. public: diff --git a/src/Enzo/enzo_EnzoRiemannImpl.hpp b/src/Enzo/enzo_EnzoRiemannImpl.hpp index 85b21cdc90..1a6a59b884 100644 --- a/src/Enzo/enzo_EnzoRiemannImpl.hpp +++ b/src/Enzo/enzo_EnzoRiemannImpl.hpp @@ -12,12 +12,6 @@ #include #include // used to check that static methods are defined -// TO DO BEFORE PR: -// - finish sorting out the primitives naming -// - clean up all of the EnzoRiemannUtils docstrings -// - I've very made very inconsistent changes to the rest of the docstrings -// - I need to entirely excise - //---------------------------------------------------------------------- struct HydroLUT { @@ -80,16 +74,15 @@ struct MHDLUT{ /// publicly defined within the scope of `ImplFunctor` (that can be acessed with /// `ImplFunctor::LUT`). /// -/// The arguments, `flux_l`, `flux_r`, `prim_l`, `prim_r`, `cons_l`, and -/// `cons_r` are constant arrays that store values associated with each -/// actively advected quantity, at the indices designated by -/// `ImplFunctor::LUT`. For each quantity, the arrays hold the its associated -/// fluxes, its inherent values, and its values when converted to conserved -/// form. Note that because the quantities are inherently specific or -/// conserved. For quantities that are inherently conserved (e.g. density) the -/// associated values in `prim_l`/`prim_r` are identical those in -/// `cons_l`/`cons_r`. The only exception is pressure. If the conserved -/// quantities hold the total_energy, then the primitives hold the pressure. +/// The arguments, `flux_l`, `flux_r`, `cons_l`, and `cons_r` are constant +/// arrays that store values associated with each actively advected integration +/// quantity, at the indices designated by `ImplFunctor::LUT`. For each +/// quantity, `flux_l`/`flux_r` hold the associated fluxes and `cons_l`/`cons_r` +/// hold the associated values (in conserved form). The arguments `prim_l` and +/// `prim_r` are similar, but hold the associated primitive values instead. The +/// indices are also deignated by `ImplFunctor::LUT`. Note that if +/// `ImplFunctor::LUT::total_energy` points to a valid index, then +/// `prim_l`/`prim_r` store pressure at that location. /// /// The left and right reconstructed pressure values are passed as `pressure_l` /// and `pressure_r`. `barotropic_eos` indicates whether the fluid equation of @@ -163,11 +156,9 @@ class EnzoRiemannImpl : public EnzoRiemann /// Constructor /// - /// @param integrable_groups A vector of integrable quantities (listed as - /// advected quantities in FIELD_TABLE). These are used as group names in - /// the Grouping objects that store field names. In effect this is used - /// to register the quantities operated on by the Riemann Solver - EnzoRiemannImpl(std::vector integrable_groups); + /// @param internal_energy Indicates whether internal_energy is an + /// integration quantity. + EnzoRiemannImpl(bool internal_energy); /// Virtual destructor virtual ~EnzoRiemannImpl(){ }; @@ -188,6 +179,12 @@ class EnzoRiemannImpl : public EnzoRiemann int stale_depth, const str_vec_t &passive_list, EFlt3DArray *interface_velocity) const; + const std::vector integration_quantities() const noexcept + { return integration_quantities_; } + + const std::vector primitive_quantities() const noexcept + { return primitive_quantities_; } + protected : //methods /// Computes the fluxes for the passively advected quantites. @@ -200,8 +197,11 @@ protected : //methods protected: //attributes - /// Names of the quantities to advect - std::vector integrable_groups_; + /// Names of the integration quantities for which the flux is computed + std::vector integration_quantities_; + + /// Names of the primitive quantities that are used in the flux calcualtion + std::vector primitive_quantities_; /// Tracks whether the internal energy needs to be computed bool calculate_internal_energy_flux_; @@ -210,52 +210,26 @@ protected : //methods //---------------------------------------------------------------------- template -EnzoRiemannImpl::EnzoRiemannImpl -(std::vector integrable_groups) +EnzoRiemannImpl::EnzoRiemannImpl(bool internal_energy) : EnzoRiemann() -{ - // Quick sanity check - integrable_groups must have density and velocity - ASSERT("EnzoRiemannImpl","integrable_groups must contain \"density\"", - std::find(integrable_groups.begin(), integrable_groups.end(), - "density") != integrable_groups.end()); - ASSERT("EnzoRiemannImpl","integrable_groups must contain \"velocity\"", - std::find(integrable_groups.begin(), integrable_groups.end(), - "velocity") != integrable_groups.end()); - - LUT::validate(); // sanity check - validate LUT's expected properties - - // Determine all quantities used by LUT - std::set lut_groups = LUT::quantity_names(); - - // TODO: add special treatment of total energy for cases with barotropic EOS - - // Check that all quantities in lut_groups are in integrable_groups - for (std::string group : lut_groups){ - ASSERT1("EnzoRiemannImpl", "\"%s\" must be a registered integrable group", - group.c_str(), - std::find(integrable_groups.begin(), integrable_groups.end(), - group) != integrable_groups.end()); - } - - calculate_internal_energy_flux_ = false; +{ + integration_quantities_ = LUT::integration_quantity_names(); - // Check that all quantities in integrable_groups or if it's "internal energy" - for (std::string group : integrable_groups) { - if (std::find(lut_groups.begin(), lut_groups.end(), group) - != lut_groups.end()){ - integrable_groups_.push_back(group); - } else if (group == "internal_energy") { - calculate_internal_energy_flux_ = true; - } else if (group == "bfield") { + for (std::string quantity : integration_quantities_){ + if (quantity == "total_energy"){ + primitive_quantities_.push_back("pressure"); + } else if (quantity == "internal_energy"){ ERROR("EnzoRiemannImpl", - "This solver isn't equipped to handle bfields. Their fluxes " - "will all be set to 0."); + "No support for a LUT directly containing \"internal_energy\""); } else { - ERROR1("EnzoRiemannImpl", "can't handle the \"%s\" integrable group", - group.c_str()); + primitive_quantities_.push_back(quantity); } } + calculate_internal_energy_flux_ = internal_energy; + if (calculate_internal_energy_flux_){ + integration_quantities_.push_back("internal_energy"); + } } //---------------------------------------------------------------------- @@ -265,11 +239,44 @@ void EnzoRiemannImpl::pup (PUP::er &p) { EnzoRiemann::pup(p); - p|integrable_groups_; + p|integration_quantities_; + p|primitive_quantities_; p|calculate_internal_energy_flux_; } //---------------------------------------------------------------------- + +/// computes the flux of a passive scalar at a given cell interface. +/// +/// @param left,right The passive scalars (as mass fractions) on the left and +/// right sides of the cell interface +/// @param density_flux The density flux at the local interface +static inline enzo_float calc_passive_scalar_flux_(enzo_float left, + enzo_float right, + enzo_float density_flux){ + // next line is equivalent to: upwind = (density_flux > 0) ? left : right; + // but is branchless + enzo_float upwind = (density_flux > 0) * left + (density_flux <= 0) * right; + return upwind * density_flux; +} + +//---------------------------------------------------------------------- + +static inline enzo_float passive_eint_flux_(enzo_float density_l, + enzo_float pressure_l, + enzo_float density_r, + enzo_float pressure_r, + enzo_float gamma, + enzo_float density_flux){ + enzo_float eint_l = EOSStructIdeal::specific_eint(density_l, pressure_l, + gamma); + enzo_float eint_r = EOSStructIdeal::specific_eint(density_r, pressure_r, + gamma); + return calc_passive_scalar_flux_(eint_l, eint_r, density_flux); +} + +//---------------------------------------------------------------------- + template void EnzoRiemannImpl::solve (EnzoEFltArrayMap &prim_map_l, EnzoEFltArrayMap &prim_map_r, @@ -294,10 +301,16 @@ void EnzoRiemannImpl::solve const EFlt3DArray pressure_array_l = prim_map_l.at("pressure"); const EFlt3DArray pressure_array_r = prim_map_r.at("pressure"); - // Check if the interface velocity must be stored - const bool store_interface_vel = (interface_velocity != nullptr); - EFlt3DArray velocity_i_bar_array; - if (store_interface_vel) { velocity_i_bar_array = *interface_velocity; } + // Check if the internal energy flux must be computed + const bool calculate_internal_energy_flux = calculate_internal_energy_flux_; + EFlt3DArray internal_energy_flux, velocity_i_bar_array; + if (calculate_internal_energy_flux){ + internal_energy_flux = flux_map.at("internal_energy"); + ASSERT("EnzoRiemannImpl::solve", + ("interface_velocity is expected to be non-NULL when computing the " + "internal energy flux"), (interface_velocity != nullptr)); + velocity_i_bar_array = *interface_velocity; + } std::array wl_arrays, wr_arrays, flux_arrays; using enzo_riemann_utils::load_array_of_fields; @@ -305,17 +318,13 @@ void EnzoRiemannImpl::solve wr_arrays = load_array_of_fields(prim_map_r, dim, true); flux_arrays = load_array_of_fields(flux_map, dim, false); - const int mz = flux_arrays[0].shape(0); - const int my = flux_arrays[0].shape(1); - const int mx = flux_arrays[0].shape(2); - ImplFunctor func; // compute the flux at all non-stale cell interfaces const int sd = stale_depth; - for (int iz = sd; iz < mz - sd; iz++) { - for (int iy = sd; iy < my - sd; iy++) { - for (int ix = sd; ix < mx - sd; ix++) { + for (int iz = sd; iz < flux_arrays[0].shape(0) - sd; iz++) { + for (int iy = sd; iy < flux_arrays[0].shape(1) - sd; iy++) { + for (int ix = sd; ix < flux_arrays[0].shape(2) - sd; ix++) { lutarray wl, wr; // get the fluid fields @@ -329,7 +338,8 @@ void EnzoRiemannImpl::solve enzo_float pressure_r = pressure_array_r(iz,iy,ix); // get the conserved quantities - lutarray Ul = enzo_riemann_utils::compute_conserved(wl, gamma); + lutarray Ul = enzo_riemann_utils::compute_conserved(wl, + gamma); lutarray Ur = enzo_riemann_utils::compute_conserved(wr, gamma); @@ -350,42 +360,22 @@ void EnzoRiemannImpl::solve flux_arrays[field_ind](iz,iy,ix) = fluxes[field_ind]; } - if (store_interface_vel){ + if (calculate_internal_energy_flux){ velocity_i_bar_array(iz,iy,ix) = interface_velocity_i; + internal_energy_flux(iz,iy,ix) = passive_eint_flux_ + (wl[LUT::density], pressure_l, wr[LUT::density], pressure_r, + gamma, fluxes[LUT::density]); } } } } - // NEED TO HANDLE INTERNAL ENERGY - - // It would be better to always have the Riemann Solver directly return an - // internal energy flux (just like the interface velocity) since the left and - // right states are basically computed for free - if (calculate_internal_energy_flux_){ - EnzoEFltArrayMap eint_l_map, eint_r_map; - eint_l_map["internal_energy"] = EFlt3DArray(mz, my, mx); - eos->eint_from_primitive(prim_map_l, eint_l_map["internal_energy"], - stale_depth); - eint_r_map["internal_energy"] = EFlt3DArray(mz, my, mx); - eos->eint_from_primitive(prim_map_r, eint_r_map["internal_energy"], - stale_depth); - - solve_passive_advection_(eint_l_map, eint_r_map, flux_map, - flux_arrays[LUT::density], stale_depth, - {"internal_energy"}); - - } - solve_passive_advection_(prim_map_l, prim_map_r, flux_map, flux_arrays[LUT::density], stale_depth, passive_list); if (LUT::has_bfields()){ // If Dedner Fluxes are required, they might get handled here - // - It would probably be better to handle this separately from the - // Riemann Solver since we already precompute L & R conserved - // AND it doesn't require wavespeed information. } } @@ -421,8 +411,8 @@ void EnzoRiemannImpl::solve_passive_advection_ enzo_float dens_flux = density_flux(iz,iy,ix); enzo_float wl = wl_arrays[key_ind](iz,iy,ix); enzo_float wr = wr_arrays[key_ind](iz,iy,ix); - enzo_float w_upwind = (dens_flux>0) ? wl : wr; - flux_arrays[key_ind](iz,iy,ix) = w_upwind * dens_flux; + flux_arrays[key_ind](iz,iy,ix) = + calc_passive_scalar_flux_(wl, wr, dens_flux); } } diff --git a/src/Enzo/enzo_EnzoRiemannLUT.hpp b/src/Enzo/enzo_EnzoRiemannLUT.hpp index 302b18da74..cfd26ee735 100644 --- a/src/Enzo/enzo_EnzoRiemannLUT.hpp +++ b/src/Enzo/enzo_EnzoRiemannLUT.hpp @@ -98,42 +98,42 @@ struct EnzoRiemannLUT{ /// @class EnzoRiemannLUT /// @ingroup Enzo /// @brief [\ref Enzo] Provides a compile-time lookup table that maps the - /// components of a subset of the actively advected quantities from - /// FIELD_TABLE to unique indices + /// components of a subset of the actively advected integration + /// quantities from FIELD_TABLE to unique indices /// /// This is a template class that provides the following features at compile /// time: /// - a lookup table (LUT) that maps the names of components of a subset - /// of the actively advected quantities defined in FIELD_TABLE to - /// unique, contiguous indices. + /// of the actively advected integration quantities defined in + /// FIELD_TABLE to unique, contiguous indices. /// - the number of quantity components included in the table - /// - a way to iterate over just the conserved quantities or specific + /// - a way to iterate over just the conserved or specific integration /// quantities values that are stored in an array using these mapping - /// - a way to query which of the actively advected quantities in - /// FIELD_TABLE are not included in the LUT + /// - a way to query which of the actively advected integration quantities + /// in FIELD_TABLE are not included in the LUT /// /// These feature are provided via the definition of publicly accessible /// integer constants in every specialization of the template class. All /// specializations have: - /// - a constant called `NEQ` equal to the number of quantity + /// - a constant called `NEQ` equal to the number of integration quantity /// components included in the lookup table /// - a constant called `specific_start` equal to the number of components - /// of conserved quantities included in the lookup table + /// of conserved integration quantities included in the lookup table /// - `qkey` constants, which include constants named for the components - /// of ALL actively advected quantities in FIELD_TABLE. A constant - /// associated with a SCALAR quantity, `{qname}`, is simply called - /// `{qname}` while constants associated with a vector quantity - /// `{qname}` are called `{qname}_i`, `{qname}_j`, and `{qname}_k`. + /// of ALL actively advected integration quantities in FIELD_TABLE. A + /// constant associated with a SCALAR quantity, `{qname}`, is simply + /// called `{qname}` while constants associated with a vector quantity, + /// `{qname}`, are called `{qname}_i`, `{qname}_j`, and `{qname}_k`. /// /// The `qkey` constants serve as both the keys of the lookup table and a /// way to check whether a component of an actively advected quantity is - /// included in the table. Their values are satisfy the following conditions: - /// - All constants named for values corresponding to quantities included - /// in the table have values of `-1` - /// - All constants named for conserved quantities have unique integer - /// values in the internal `[0,specific_start)` - /// - All constants named for specific quantities have unique integer - /// values in the interval `[specific_start, NEQ)` + /// included in the table. Their values satisfy the following conditions: + /// - All constants named for values corresponding to quantities NOT + /// included in the lookup table have values of `-1` + /// - All constants named for conserved integration quantities have unique + /// integer values in the internal `[0,specific_start)` + /// - All constants named for specific integration quantities have unique + /// integer values in the interval `[specific_start, NEQ)` /// /// The lookup table is always expected to include density and the 3 velocity /// components. Although it may not be strictly enforced (yet), the lookup @@ -145,16 +145,16 @@ struct EnzoRiemannLUT{ /// the above requirements are specified. /// /// @tparam InputLUT Type that defines the `NEQ` constant, the - /// `specific_start` constants and the `qkey` constants that correspond - /// to the actively advected quantities that are actually included in the - /// lookup table. All of the constant values will be directly reused by - /// the resulting template specialization and therefore must meet the - /// criteria defined above. Any undefined `qkey` constants are assumed to - /// not be included in the lookup table and will be defined within the - /// template specialization to have values of `-1`. This type can be - /// implemented with either an unscoped enum OR a class that includes - /// some combination of publicly accessible unscoped enums and static - /// constant member variables. + /// `specific_start` constant, and the `qkey` constants that correspond + /// to the actively advected integration quantities that are actually + /// included in the lookup table. All of the constant values will be + /// directly reused by the resulting template specialization and + /// therefore must meet the criteria defined above. Any undefined `qkey` + /// constants are assumed to not be included in the lookup table and will + /// be defined within the template specialization to have values of `-1`. + /// This type can be implemented with either an unscoped enum OR a class + /// that includes some combination of publicly accessible unscoped enums + /// and static constant member variables. /// /// @par Examples /// For the sake of the example, let's assume we have a type `MyInputLUT` @@ -181,15 +181,15 @@ struct EnzoRiemannLUT{ /// /// @par /// As an example, imagine that the total kinetic energy density needs to be - /// computed at a single location from an values stored in an array, `prim`, - /// of type `lutarray>`. The resulting code to + /// computed at a single location from the values stored in an array, `integ`, + /// of type `lutarray>`. The resulting code to /// do that might look something like: /// @code - /// using LUT = EnzoRiemannLUT; - /// enzo_float v2 = (prim[LUT::velocity_i] * prim[LUT::velocity_i] + - /// prim[LUT::velocity_j] * prim[LUT::velocity_j] + - /// prim[LUT::velocity_k] * prim[LUT::velocity_k]); - /// enzo_float kinetic = 0.5 * prim[LUT::density] * v2; + /// using LUT = EnzoRiemannLUT; + /// enzo_float v2 = (integ[LUT::velocity_i] * integ[LUT::velocity_i] + + /// integ[LUT::velocity_j] * integ[LUT::velocity_j] + + /// integ[LUT::velocity_k] * integ[LUT::velocity_k]); + /// enzo_float kinetic = 0.5 * integ[LUT::density] * v2; /// @endcode /// /// @par @@ -218,20 +218,20 @@ struct EnzoRiemannLUT{ /// /// @par /// This final example will show the value of grouping the conserved and - /// specific quantities. To implement some Riemann solvers, it is useful to - /// have a generic, reusable function that constructs an array of all - /// quantities that are in the lookup table in their conserved form. The + /// specific integration quantities. To implement some Riemann solvers, it is + /// useful to have a generic, reusable function that constructs an array of + /// all quantities that are in the lookup table in their conserved form. The /// following function was used to accomplish this at one point /// @code /// template - /// lutarray compute_conserved(const lutarray prim) noexcept + /// lutarray compute_conserved(const lutarray integr) noexcept /// { /// lutarray cons; /// for (std::size_t i = 0; i < LUT::specific_start; i++) { - /// cons[i] = prim[i]; + /// cons[i] = integr[i]; /// } /// for (std::size_t i = LUT::specific_start; i < LUT::NEQ; i++) { - /// cons[i] = prim[i] * prim[LUT::density]; + /// cons[i] = integr[i] * integr[LUT::density]; /// } /// return cons; /// } @@ -249,8 +249,9 @@ struct EnzoRiemannLUT{ }; /// The entry with the minimum (non-negative) index corresponding to a - /// specific quantity. All conserved quantity entries must have smaller - /// indices (defaults to -1 if not explicitly specified in InputLUT) + /// specific integration quantity. All conserved integration quantity entries + /// must have smaller indices (defaults to -1 if not explicitly specified in + /// InputLUT) static constexpr std::size_t specific_start = LUTIndexForward_::forward_specific_start::value; @@ -281,16 +282,18 @@ struct EnzoRiemannLUT{ public: //associated static functions - /// returns a set of integrable quantities included in the InputLUT - static std::set quantity_names() noexcept; + /// returns a vector of integration quantity names included in the InputLUT + static std::vector integration_quantity_names() noexcept; /// returns whether the LUT has any bfields static constexpr bool has_bfields(){ return (qkey::bfield_i>=0) || (qkey::bfield_j>=0) || (qkey::bfield_k>=0); } - /// for each quantity in FIELD_TABLE, this passes the associated passes the - /// associated name and index from the wrapped InputLUT to the function `fn` + /// for each actively advected scalar integration quantity and component of + /// an actively advected vector integration quantity in FIELD_TABLE, this + /// passes the associated name and index from the wrapped InputLUT to the + /// function `fn` /// /// The function should expect (std::string name, int index). In cases where /// quantites in FIELD_TABLE do not appear in the wrapped InputLUT, an index @@ -299,7 +302,7 @@ struct EnzoRiemannLUT{ static void for_each_entry(Function fn) noexcept; /// a function that performs a check to make sure that the InputLUT satisfies - /// all assumptions. If it doesn't, an error is raised + /// all assumptions. If it doesn't, the program aborts (with an Error message) static void validate() noexcept; }; @@ -338,18 +341,23 @@ void EnzoRiemannLUT::for_each_entry(Function fn) noexcept{ //---------------------------------------------------------------------- template -std::set EnzoRiemannLUT::quantity_names() noexcept +std::vector EnzoRiemannLUT::integration_quantity_names() + noexcept { - std::set set; + std::vector vec; auto fn = [&](std::string name, int index) { if (index < 0) {return;} - set.insert(EnzoCenteredFieldRegistry::get_actively_advected_quantity_name - (name,true)); + std::string quantity = + EnzoCenteredFieldRegistry::get_actively_advected_quantity_name(name, + true); + if (std::find(vec.begin(), vec.end(), quantity) == vec.end()){ + vec.push_back(quantity); + } }; EnzoRiemannLUT::for_each_entry(fn); - return set; + return vec; } //---------------------------------------------------------------------- diff --git a/src/Enzo/enzo_EnzoRiemannUtils.hpp b/src/Enzo/enzo_EnzoRiemannUtils.hpp index f318f5575c..1d6659d559 100644 --- a/src/Enzo/enzo_EnzoRiemannUtils.hpp +++ b/src/Enzo/enzo_EnzoRiemannUtils.hpp @@ -18,6 +18,31 @@ //---------------------------------------------------------------------- +struct EOSStructIdeal{ + // In the future it might make sense to use this under the hood of the + // EnzoEOSIdeal class. + // + // It might also make sense to make this into POD struct that get's passed to + // an Impl Functor's operator() method. If the structure is simple enough, + // the overhead of constructing/destroying can be optimized out + + static inline enzo_float specific_eint(enzo_float density, + enzo_float pressure, + enzo_float gamma) noexcept + { return pressure / ( (gamma - 1.0) * density); } + + static inline enzo_float eint_dens(enzo_float density, enzo_float pressure, + enzo_float gamma) noexcept + { return pressure / (gamma - 1.0); } + + static inline enzo_float sound_speed(enzo_float density, enzo_float pressure, + enzo_float gamma) noexcept + { return std::sqrt(gamma * pressure / density); } + +}; + +//---------------------------------------------------------------------- + namespace enzo_riemann_utils{ template @@ -60,14 +85,15 @@ namespace enzo_riemann_utils{ } if (LUT::total_energy >= 0) { // overwrite the total energy index + enzo_float density = prim[LUT::density]; enzo_float pressure = prim[LUT::total_energy]; - enzo_float internal_edens = pressure / (gamma - 1.); + enzo_float internal_edens = EOSStructIdeal::eint_dens(density, pressure, + gamma); - enzo_float density = prim[LUT::density]; enzo_float vi = prim[LUT::velocity_i]; enzo_float vj = prim[LUT::velocity_j]; enzo_float vk = prim[LUT::velocity_k]; - enzo_float kinetic_edens = 0.5 * density * (vi*vi + vj*vj + vk *vk); + enzo_float kinetic_edens = 0.5 * density * (vi*vi + vj*vj + vk*vk); enzo_float magnetic_edens = mag_pressure(prim); @@ -79,13 +105,6 @@ namespace enzo_riemann_utils{ //---------------------------------------------------------------------- - template - inline enzo_float sound_speed(const lutarray prim_vals, - enzo_float pressure, enzo_float gamma) noexcept - { return std::sqrt(gamma * pressure / prim_vals[LUT::density]); } - - //---------------------------------------------------------------------- - template inline enzo_float fast_magnetosonic_speed(const lutarray prim_vals, enzo_float pressure, @@ -98,7 +117,9 @@ namespace enzo_riemann_utils{ // TODO: optimize calc of cs2 to omit sqrt and pow // can also skip the calculation of B2 by checking if // LUT::bfield_i, LUT::bfield_j, LUT::bfield_k are all negative - enzo_float cs2 = std::pow(sound_speed(prim_vals, pressure, gamma),2); + enzo_float cs = EOSStructIdeal::sound_speed(prim_vals[LUT::density], + pressure, gamma); + enzo_float cs2 = std::pow(cs,2); enzo_float B2 = (bi*bi + bj*bj + bk *bk); if (B2 == 0){ return std::sqrt(cs2); @@ -127,7 +148,9 @@ namespace enzo_riemann_utils{ // TODO: optimize calc of cs2 to omit sqrt and pow // can also skip the calculation of B2 by checking if // LUT::bfield_i, LUT::bfield_j, LUT::bfield_k are all negative - enzo_float cs2 = std::pow(sound_speed(prim_vals, pressure, gamma),2); + enzo_float cs = EOSStructIdeal::sound_speed(prim_vals[LUT::density], + pressure, gamma); + enzo_float cs2 = std::pow(cs,2); enzo_float B2 = (bi*bi + bj*bj + bk *bk); if (B2 == 0){ return std::sqrt(cs2);