diff --git a/doc/source/devel/hydro_infrastructure.rst b/doc/source/devel/hydro_infrastructure.rst index c8cfa85032..1abdad8845 100644 --- a/doc/source/devel/hydro_infrastructure.rst +++ b/doc/source/devel/hydro_infrastructure.rst @@ -28,36 +28,49 @@ Shorthand Terms We briefly define a few terms that are used throughout the documentation and codebase -Integrable/Reconstructable Quantities +Integration/Primitive Quantities ------------------------------------- Throughout this guide and the relevant sections of the codebase, we -refer to quantities as reconstructable (i.e. they are used for -reconstructing left/right interface states) and integrable (i.e. the -primary quantities that are evolved by the integrator). There is a -high degree of overlap between these categories and the precise -categorization depends on the equation of state. As of now all -integrable quantities are "conserved" or "specific" (a quantity like -velocity that when multiplied by density becomes conserved). +categorize quantities as integration quantities and primitives. -As an example, the categorization of the quantities for an ideal, -adiabatic gas are: +The integration quantities are the cell-centered quantities that +Enzo-E, as a whole, uses to describe the state of the fluid. In other +words, these are the hydro/MHD quantities that Enzo-E expects to be +integrated from one time-step to the next. The integration quantities +include all passive scalars in their "conserved" form (i.e. as +densities). All integration quantities can basically be +subcategorized as either "conserved" or "specific" (a specific +quantity like velocity that becomes conserved after multiplication by +density). In cases using the dual energy formalism, we treat the +specific internal energy as an integration quantity even though the +internal energy density is not technically conserved (it requires +source terms). - * density - reconstructable and integrable +The primitive quantities follow the normal textbook definition. We use +the primitives internally within the hydro/MHD solver for +reconstruction. The primitives include all passive scalars in their +"specific" form (i.e. as mass fractions). Note that some quantities +(like density or velocity) are categorized as both an integration +quantity and a primitive. - * velocity - reconstructable and integrable +To provide a more concrete example, we categorize the quantities related +to an ideal, adiabatic gas: - * pressure - only reconstructable + * density - integration and primitive - * (specific) total energy - only integrable + * velocity - integration and primitive - * magnetic field - reconstructable and integrable + * pressure - only primitive -*Note: Both the reconstructable and integrable quantities are -frequently referred to as primitives throughout the codebase as -primitives. This is mostly historical and mainly refers to the fact -that the collection of quantities are not all of the quantities are -"conserved" (at least some subset of them are "specific").* + * (specific) total energy - only integration + + * magnetic field - integration and primitive + +If using the dual energy formalism, we would categorize the (specific) +internal energy as an integration quantity. As of now, there wouldn't +be a primitive counterpart to the internal energy since it can be +computed from the reconstructed density and pressure. stale depth ----------- @@ -147,13 +160,11 @@ the row for the "velocity" quantity registers the ``"velocity_x"``, At present, the registry currently provides operations: - * for building ``Grouping`` objects that contain registered quantity - fields. - * to access quantity properties registerred in ``FIELD_TABLE`` at + * to access quantity properties registered in ``FIELD_TABLE`` at runtime - * provide a list of known groups that can be used in the input file + * to provide a list of known groups that can be used in the input file to identify fields as passively advected scalars (as of now, the - only such group is ``"colour"``). + only such group is ``"color"``). ============== General Design @@ -172,17 +183,15 @@ classes include: * ``EnzoEquationOfState`` - encapsulates many of the operations related to the fluid's equation of state (e.g. computing pressure, - converting reconstructable quantities to integrable quantities and - vice-versa) + converting the integration quantities or primitives) * ``EnzoReconstructor`` - encapsulates interpolation algorithms to - reconstruct left/right interface states of from cell-centered - values + reconstruct left/right interface states of cell-centered values * ``EnzoRiemann`` - encapsulates various Rimann Solver algorithms - * ``EnzoIntegrableUpdate`` - encapsulates the operation of updating - integrable quantities after a (partial) time-step. + * ``EnzoIntegrationQuanUpdate`` - encapsulates the operation of + updating integration quantities after a (partial) time-step. * ``EnzoBfieldMethod`` - encapsulates operations related to integrating magnetic fields that are not performed by the other operation classes. @@ -191,27 +200,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 integrable or reconstructable quantities (other than passively -advected scalars) are *registered* at construction. - - * The names of all reconstructable quantites that get registered - in the construction of ``EnzoRiemann`` must share a name - with the registered quantities in ``FIELD_TABLE``. - - * All registered integrable quantity names in the construction of - ``EnzoRiemann`` or ``EnzoIntegrableUpdate`` 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 a far more independent manner. Because all fields storing passively advected scalars are not necessarily known when initializing a hydro/MHD integrator (i.e. @@ -247,54 +253,62 @@ Specific Usage In the context of this toolkit, the keys of an ``EnzoEFltArrayMap`` -are usually the names of a scalar quantity (like ``"density"``) -or component of a vector quantity (like ``"velocity_x"``). Each key -is paired with an instance of ``EFlt3DArray`` that stores data -associated data. Below, we provide a description of the main -uses of ``EnzoEFltArrayMap`` by the provided operation classes: +are usually the names of a scalar quantity (like ``"density"``) or +component of a vector quantity (like ``"velocity_x"``). Each key is +paired with an instance of ``EFlt3DArray`` that stores associated +data. To simplify logic, arrays are not aliased between separate maps. +Below, we provide a description of the main uses of +``EnzoEFltArrayMap`` by the provided operation classes: + + 1. Map of cell-centered integration quantities. + + * This has keys named for all integration scalar quantities and + components of integration vector quantities. The associated + arrays hold the values of the cell-centered quantities at a + given time. + + * This also contains key-value pairs for passively advected + scalars. In this context, the passive scalars are stored in + "conserved" form. - 1. Map of cell-centered quantities. + * In a predictor-corrector scheme (like VL+CT), we might have + multiple maps used to store values at different partial + timesteps. - * This has quantity/quantity-component keys named for all - integrable and reconstructable quantities used by the - integrator. The associated arrays hold the values of the - cell-centered quantities at a given time. We currently store - integrable and reconstructed quantities together due to the - high degree of overlap between each category. + 2. Map of cell-centered primitive quantities. + + * This map is used to temporarily store the cell-centered + primitive quantities for use in reconstruction. * This also contains key-value pairs for passively advected - scalars. Note that in this context, the passive scalars are - usually represented in "specific" form. For reference, the - general convention throughout Enzo-E is to represent primarily - passively advected scalars in "conserved" form (as mass - densities) outside of hydrodynamic integrator methods and to - convert them to "specific" form (mass fractions) within - hydrodynamic integrator methods. - - 2. Map of temporary cell-centered for tracking the total change - in a quantity over a timestep. - - * This map holds key-array pairs named for all integrable - quantities and groups of passively advected scalars. For each - (partial) timestep, these arrays are used to accumulate the - total change in the conserved form of each quantity. This - includes the flux divergence and the contributions from source - terms. At the end of the (partial) timestep, these are used to - actually update the values of the integrable quantities - - 3. Map of reconstructed left/right quantites - - * 2 instances of ``EnzoEFltArrayMap`` are used to respectively hold - the reconstructed left and right interface quantities. This should - share have the same keys that are described for the first category - of maps. - * These maps are frequently passed to instances of ``EnzoReconstructor`` - to store the reconstructed passively advected scalars and - reconstructable quantities. They are then usually passed to - ``EnzoEquationOfState`` to compute and store the reconstructed - integrable quantities and reconstructed pressure. Then, these are - frequently passed to ``EnzoRiemann`` to compute fluxes from the - integrable quantities and the passively advected scalars. + scalars. In this context, the passive scalars are stored in + "specific" form. + + * Quantities in both the primitive map and integration map should + NOT be aliases of each other. They should be deepcopies instead. + + 3. Map of temporary cell-centered values for tracking the total + change in a quantity over a timestep. + + * This map holds key-array pairs named for all integration + quantities. For each (partial) timestep, these arrays are used + to accumulate the total change in the conserved form of each + quantity. This includes the flux divergence and the + contributions from source terms. At the end of the (partial) + timestep, these are used to actually update the values of the + integration quantities + + 4. Map of reconstructed left/right primitive quantites + + * 2 instances of ``EnzoEFltArrayMap`` are used to respectively + hold the reconstructed left and right interface primitive + quantities. This should share have the same keys that are + described for the second category of maps. + * These maps are frequently passed to instances of + ``EnzoReconstructor`` to store the reconstructed passively + advected scalars and primitive quantities. Then, these are + frequently passed to ``EnzoRiemann`` to compute fluxes for + the integration quantities and passively advected scalars. * Although this inherently represents data centered on the faces of the mesh, the contained arrays should formally have the shape required to hold cell-centered data. This is done to facillitate the reuse of @@ -302,7 +316,7 @@ uses of ``EnzoEFltArrayMap`` by the provided operation classes: means that there is always some unused allocated memory at the end of one of the dimensions. - 4. Maps of Riemann Flux fields + 5. Maps of Riemann Flux fields * An instance of this kind of map is required for each dimension and is used to hold the face-centered fluxes along @@ -312,22 +326,20 @@ uses of ``EnzoEFltArrayMap`` by the provided operation classes: normally holds ``n`` elements (including ghost zones) along axis ``i``, then an array used to store fluxes along axis ``i`` should hold ``n-1`` elements along axis ``i``. + * This should have all of the same keys that are in the the first + category of maps. * This kind of map should contain keys named for all passively advected - scalars and registered integrable quantities. The set of keys in these + scalars and registered integration quantities. The set of keys in these maps should be identical to the set of keys in the first category of maps, regardless of whether a quantity is "specific" or "conserved" (e.g. the map will hold a "velocity_x" key even though the associated array stores the x-component of the momentum density flux). -Note that the ``EnzoEquationOfState`` and ``EnzoIntegrableUpdate`` -classes additionally require a ``EnzoEFltArrayMap`` object that hold the -passively advected scalars in conserved form. - -In general, the use of ``EnzoEFltArrayMap`` objects with common sets of keys -helps simplify the implementation of various methods (e.g. the -cell-centered array associated with "density" is used to compute the -reconstruct values that are stored in the fields of the "density" -group in the reconstructed grouping). +In general, the use of ``EnzoEFltArrayMap`` objects with common sets +of keys helps simplify the implementation of various methods (e.g. the +cell-centered array associated with "density" is used to reconstruct +values that are stored in the fields of the "density" +array in the primitive map). ================= @@ -393,28 +405,23 @@ Returns the thermal pressure floor. .. code-block:: c++ - apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integrable_map, + apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integration_map, int stale_depth); This method applies the applies the pressure floor to the total_energy -array specified in ``integrable_map``. If using the dual-energy formalism +array specified in ``integration_map``. If using the dual-energy formalism the floor is also applied to the internal energy (also specified in -``integrable_map``) and synchronizes the internal energy with the total +``integration_map``) and synchronizes the internal energy with the total energy. If the equation of state is barotropic, this should do nothing. .. code-block:: c++ - void pressure_from_integrable(EnzoEFltArrayMap &integrable_map, - const EFlt3DArray &pressure, - EnzoEFltArrayMap &conserved_passive_map, - int stale_depth); + void pressure_from_integration(EnzoEFltArrayMap &integration_map, + const EFlt3DArray &pressure, + int stale_depth); -This method computes the pressure from the integrable quantities -(stored in ``integrable_map``) and stores the result in ``pressure``. -``conserved_passive_map`` should include the passive scalars in -conserved form. The last argument currently doesn't do anything -and will only be important if Grackle is in use (with -``primordial_chemistry>1``). +This method computes the pressure from the integration quantities +(stored in ``integration_map``) and stores the result in ``pressure``. *In principle this should wrap* ``EnzoComputePressure``, *but currently that is not the case. Some minor refactoring is needed to @@ -426,59 +433,17 @@ internal energy.* .. code-block:: c++ - void pressure_from_reconstructable(EnzoEFltArrayMap &reconstructable, - EFlt3DArray &pressure, - int stale_depth); - -This method computes the pressure from the reconstructable quantities -(stored in ``reconstructable``) and stores the result in ``pressure``. - -Note: for a non-barotropic equation of state, pressure is considered a -reconstructable quantity. In that case, if the pressure array in -``reconstructable`` is an alias of ``pressure``, nothing happens. -However, if the arrays aren't aliases, then the values are simply -copied between arrays. - -.. code-block:: c++ - - void reconstructable_from_integrable - (EnzoEFltArrayMap &integrable, EnzoEFltArrayMap &reconstructable, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth, - const std::vector &passive_list); - -This method is responsible for computing the reconstructable -quantities (to be held in ``reconstructable``) from the -integrable quantities (stored in ``integrable``). Note that -because of the high degree of overlap between the quantities in each -category, the overlapping quantities are assumed to be represented by -aliased arrays. Entries in ``passive_list`` specifies the passively -advected scalars that should be present in both maps. -The conserved form of the passively advected scalars -must be provided (stored in ``conserved_passive_map``) in case the -equation of state is barotropic and Grackle is in use. - -For a barotropic equation of state, this nominally does nothing while -for a non-barotropic equation of state, this just computes pressure by -calling ``EnzoEquationOfState::pressure_from_integrable``. - -.. code-block:: c++ - - void integrable_from_reconstructable - (EnzoEFltArrayMap &reconstructable, EnzoEFltArrayMap &integrable, + void primitive_from_integration + (EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &primitive_map, int stale_depth, const std::vector &passive_list); -This method computes the integrable quantities (to be held in -``integrable``) from the reconstructable quantities (stored in -``reconstructable``). Again, because of the high degree of -overlap between the quantities in each category, the overlapping -quantities are assumed to be represented by aliased quantities. -Entries in ``passive_list`` specifies the passively -advected scalars that should be present in both maps. - -For a barotropic equation of state, this nominally does nothing, while -for a non-barotropic equation of state, this nominally just computes -specific total energy. If the dual-energy formalism is in use this -also computes the internal energy. +This method is responsible for computing the primitive quantities (to +be held in ``primitive_map``) from the integration quantities (stored +in ``integration_map``). Non-passive scalar quantities appearing in +both ``integration_map`` and ``primitive_map`` are simply deepcopied +and passive scalar quantities are converted from conserved-form to +specific form. For a non-barotropic EOS, this also computes pressure +(by calling ``EnzoEquationOfState::pressure_from_integration``) *In the future, it might be worth considering making this into a subclass of Cello's ``Physics`` class. If that is done, it may be advisable to @@ -514,11 +479,12 @@ To get a pointer to an instance of a concrete implementation of std::string name, enzo_float theta_limiter); The factory method requires that we register the names of the -reconstructable quantities via ``active_reconstructable_quantities`` -and specify the name of the reconstruction algorithm, ``name``. Note -that the names of the reconstructable quantites should match -quantities specified in ``FIELD_TABLE`` ; for more details about -``FIELD_TABLE``, see :ref:`Centered-Field-Registry` +non-passive scalar primitive quantities that are are to be +reconstructed via ``active_reconstructed_quantities``. We specify +the name of the reconstruction algorithm, ``name``. Note that the +names of the primitive quantites should match quantities specified in +``FIELD_TABLE`` ; for more details about ``FIELD_TABLE``, see +:ref:`Centered-Field-Registry` Public Interface ---------------- @@ -531,26 +497,26 @@ The main interface function provided by this class is: EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, int stale_depth, const std::vector& passive_list); -This function takes the cell-centered reconstructable primtive -quantities (specified by the contents of ``prim_map``) and computes -the left and right reconstructed states (the results are stored in -``priml_map`` and ``primr_map``) along the dimension specifed by -``dim``. If dim has a value of ``0``/ ``1``/ ``2`` then the values are -reconstructed along the x-/y-/z-axis. ``stale_depth`` indicates the -current stale_depth for the supplied cell-centered quantities (prior -to reconstruction). Note that the arrays in ``priml_map`` and -``primr_map`` should have arrays that are large enough to store -cell-centered quantitites so that they can be reused to hold the -face-centered fields along each dimension. ``passive_list`` is used to -specify the names (keys) of the passively advected quantities that are -to be reconstructed. +This function takes the cell-centered primtive quantities (specified +by the contents of ``prim_map``) and computes the left and right +reconstructed states (the results are stored in ``priml_map`` and +``primr_map``) along the dimension specifed by ``dim``. If dim has a +value of ``0``/ ``1``/ ``2`` then the values are reconstructed along +the x-/y-/z-axis. ``stale_depth`` indicates the current stale_depth +for the supplied cell-centered quantities (prior to +reconstruction). Note that the arrays in ``priml_map`` and +``primr_map`` should be large enough to store cell-centered +quantitites so that they can be reused to hold the face-centered +fields along each dimension. ``passive_list`` is used to specify the +names (keys) of the passively advected quantities that are to be +reconstructed. The ``int EnzoReconstructor::immediate_staling_rate()`` method is provided to determine the amount by which the stale depth increases immediately after reconstruction, for a given algorithm. The ``int EnzoReconstructor::delayed_staling_rate()`` method returns how much the stale depth increases after adding flux divergence, computed from -the reconstructed values, to the integrable quantities (this is +the reconstructed values, to the integration quantities (this is normally 1). Finally ``int EnzoReconstructor::total_staling_rate()`` gives the sum of the results yielded by the prior 2 methods. @@ -617,38 +583,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 @@ -674,11 +649,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 @@ -698,26 +676,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 @@ -725,11 +702,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, @@ -756,9 +731,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()``. @@ -781,14 +756,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 @@ -820,45 +798,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 @@ -870,17 +848,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``. @@ -889,9 +867,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 @@ -900,16 +880,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 @@ -940,24 +920,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 ------------------ @@ -982,41 +952,41 @@ The additional steps for implementing a new Riemann solver by speciallizing ``EnzoRiemannImpl`` that uses the new ``ImplFunctor`` class (e.g. ``using EnzoRiemannHLLD = EnzoRiemannImpl;``). -============================== -Updating integrable quantities -============================== +=============================== +Updating integration quantities +=============================== -The ``EnzoIntegrableUpdate`` class has been provided to encapsulate -the operation of updating integrable quantities after a (partial) +The ``EnzoIntegrationQuanUpdate`` class has been provided to encapsulate +the operation of updating integration quantities after a (partial) time-step. The operation was factored out of the ``EnzoMethodMHDVlct`` class since it appear in all Godunov solvers. -The constructor for ``EnzoIntegrableUpdate`` has the following +The constructor for ``EnzoIntegrationQuanUpdate`` has the following signature: .. code-block:: c++ - EnzoIntegrableUpdate(std::vector integrable_groups, - bool skip_B_update) + EnzoIntegrationQuanUpdate(std::vector integration_groups, + bool skip_B_update) The function requires that we: - * register the names of the integrable quantities (with - ``integrable_groups``) + * register the names of the integration quantities (with + ``integration_groups``) * indicate whether the update to the magnetic field should be skipped. -The names of the integrable quantites should match the names specified +The names of the integration quantites should match the names specified in ``FIELD_TABLE``; see :ref:`Centered-Field-Registry` for more details. The update to the magnetic field should be skipped when Constrained Transport is in use (since the magnetic field update is handled separately). If the magnetic field is not specified as an -integrable quantity, then the value specified for ``skip_B_update`` is +integration quantity, then the value specified for ``skip_B_update`` is unimportant The following method is used to compute the change in (the conserved -form of) the integrable and passively advected quantites due to the -flux divergence along dimension ``dim`` over the (partial) imestep +form of) the integration and passively advected quantites due to the +flux divergence along dimension ``dim`` over the (partial) timestep ``dt``. The arrays in ``dUcons_map`` are used to accumulate the total changes in these quantities. ``passive_list`` lists the names (keys) of the passively advected scalars. @@ -1038,27 +1008,27 @@ argument is used in the same way as the previous function. void clear_dUcons_group(EnzoEFltArrayMap &dUcons_map, enzo_float value, const std::vector &passive_list) const; -The method used to actually add the accumulated change in the integrable +The method used to actually add the accumulated change in the integration (specified in ``dUcons_map``) to the values of the -integrable quantities from the start of the timestep (specificed by -``initial_integrable_map``) has the following signature: +integration quantities from the start of the timestep (specificed by +``initial_integration_map``) has the following signature: .. code-block:: c++ void update_quantities - (EnzoEFltArrayMap &initial_integrable_map, EnzoEFltArrayMap &dUcons_map, - EnzoEFltArrayMap &out_integrable_map, + (EnzoEFltArrayMap &initial_integration_map, EnzoEFltArrayMap &dUcons_map, + EnzoEFltArrayMap &out_integration_map, EnzoEFltArrayMap &out_conserved_passive_scalar, EnzoEquationOfState *eos, int stale_depth, const std::vector &passive_list) const; The fields included in ``dUcons_map`` should include contributions from both the flux divergence AND source terms. The results for the -actively advected quanties are stored in ``out_integrable_map`` and +actively advected quanties are stored in ``out_integration_map`` and the results for the passively advected scalars are stored in conserved form in the arrays held by ``out_conserved_passive_scalar`` (note that the initial values of the passive scalars specified in -``initial_integrable_map`` are in specific form). +``initial_integration_map`` are in specific form). ========================== Magnetic Field Integration @@ -1191,9 +1161,9 @@ magnetic field values. In ``EnzoBfieldMethodCT`` this will also update the face-centered magnetic field values (it assumes that ``identify_upwind`` was called once for each dimension and uses the stored data). When using this -alongside ``EnzoIntegrableUpdate``, care needs to be taken about the +alongside ``EnzoIntegrationQuanUpdate``, care needs to be taken about the order in which this method is called relative to -``EnzoIntegrableUpdate::update_quantities`` that accounts for the time +``EnzoIntegrationQuanUpdate::update_quantities`` that accounts for the time when floors are applied to the total energy. Descriptor Methods diff --git a/input/vlct/run_dual_energy_cloud_test.py b/input/vlct/run_dual_energy_cloud_test.py index 5018666cb6..a3069ec1f2 100644 --- a/input/vlct/run_dual_energy_cloud_test.py +++ b/input/vlct/run_dual_energy_cloud_test.py @@ -83,9 +83,9 @@ 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) + 'hlle_cloud', 3.2e-13) n_passed = np.sum(r) n_tests = len(r) print("{:d} Tests passed out of {:d} Tests.".format(n_passed,n_tests)) diff --git a/src/Enzo/_enzo.hpp b/src/Enzo/_enzo.hpp index 3acbd607e7..6d3a18bdea 100644 --- a/src/Enzo/_enzo.hpp +++ b/src/Enzo/_enzo.hpp @@ -210,7 +210,7 @@ extern "C" { #include "enzo_EnzoLazyPassiveScalarFieldList.hpp" #include "enzo_EnzoEquationOfState.hpp" #include "enzo_EnzoEOSIdeal.hpp" -#include "enzo_EnzoIntegrableUpdate.hpp" +#include "enzo_EnzoIntegrationQuanUpdate.hpp" #include "enzo_EnzoReconstructor.hpp" #include "enzo_EnzoReconstructorNN.hpp" #include "enzo_EnzoReconstructorPLM.hpp" diff --git a/src/Enzo/enzo.CI b/src/Enzo/enzo.CI index 057bffa04e..9e4c171be5 100644 --- a/src/Enzo/enzo.CI +++ b/src/Enzo/enzo.CI @@ -99,7 +99,7 @@ module enzo { PUPable EnzoMethodMHDVlct; PUPable EnzoEOSIdeal; - PUPable EnzoIntegrableUpdate; + PUPable EnzoIntegrationQuanUpdate; PUPable EnzoReconstructorNN; // EnzoReconstructorPLMEnzoRKLim PUPable EnzoReconstructorPLM; diff --git a/src/Enzo/enzo_EnzoBfieldMethod.hpp b/src/Enzo/enzo_EnzoBfieldMethod.hpp index 99d987cd6d..8087bc157c 100644 --- a/src/Enzo/enzo_EnzoBfieldMethod.hpp +++ b/src/Enzo/enzo_EnzoBfieldMethod.hpp @@ -147,10 +147,10 @@ class EnzoBfieldMethod : public PUP::able /// For constrained transport, this includes both the face-centered and the /// cell-centered bfields (the face-centered bfield updates are handled /// completely internally). The velocity and bfield components stored in this - /// `cur_prim_map` map are used to compute the cell-centered E-field. + /// `cur_integration_map` map are used to compute the cell-centered E-field. /// - /// @param[in] cur_prim_map Map containing the current values of the - /// cell-centered integrable quantities (before they have been updated + /// @param[in] cur_integration_map Map containing the current values of the + /// cell-centered integration quantities (before they have been updated /// over the current timestep). /// @param[in] xflux_map,yflux_map,zflux_map Maps containing the values of /// the fluxes computed along the x, y, and z directions. The function @@ -164,9 +164,9 @@ class EnzoBfieldMethod : public PUP::able /// @param[in] stale_depth indicates the current stale_depth for the /// supplied quantities. This should nominally be the same as the stale /// depth used to compute the fluxes and that is passed to - /// EnzoIntegrableUpdate::update_quantities. + /// EnzoIntegrationQuanUpdate::update_quantities. virtual void update_all_bfield_components - (EnzoEFltArrayMap &cur_prim_map, EnzoEFltArrayMap &xflux_map, + (EnzoEFltArrayMap &cur_integration_map, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, EnzoEFltArrayMap &out_centered_bfield_map, enzo_float dt, int stale_depth) noexcept=0; diff --git a/src/Enzo/enzo_EnzoBfieldMethodCT.cpp b/src/Enzo/enzo_EnzoBfieldMethodCT.cpp index 67b3912268..d9093d1fe9 100644 --- a/src/Enzo/enzo_EnzoBfieldMethodCT.cpp +++ b/src/Enzo/enzo_EnzoBfieldMethodCT.cpp @@ -36,11 +36,12 @@ void EnzoBfieldMethodCT::register_target_block_ } // get the standard size of cell-centered fields - int mz = bfieldi_l_[0].shape(0); // interface bfield_z cell-centered along z - int my = bfieldi_l_[0].shape(1); // interface bfield_y cell-centered along y + // interface bfield_x is cell-centered along y and z + const int mz = bfieldi_l_[0].shape(0); + const int my = bfieldi_l_[0].shape(1); // interface bfield_x is face-centered along x (and the array includes // values on external faces of the grid) - int mx = bfieldi_l_[0].shape(2) - 1; + const int mx = bfieldi_l_[0].shape(2) - 1; // allocate arrays to hold weight values for each dimension. For a given // dimension, the array tracks the upwind/downwind direction on the cell @@ -202,7 +203,7 @@ void EnzoBfieldMethodCT::identify_upwind(const EnzoEFltArrayMap &flux_map, //---------------------------------------------------------------------- void EnzoBfieldMethodCT::update_all_bfield_components -(EnzoEFltArrayMap &cur_prim_map, EnzoEFltArrayMap &xflux_map, +(EnzoEFltArrayMap &cur_integration_map, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, EnzoEFltArrayMap &out_centered_bfield_map, enzo_float dt, int stale_depth) noexcept @@ -224,12 +225,11 @@ void EnzoBfieldMethodCT::update_all_bfield_components } // First, compute the edge-centered Electric fields (each time, it uses - // the current integrable quantities) - EnzoBfieldMethodCT::compute_all_edge_efields(cur_prim_map, xflux_map, - yflux_map, zflux_map, - center_efield_, - edge_efield_l_, weight_l_, - stale_depth); + // the current integration quantities) + EnzoBfieldMethodCT::compute_all_edge_efields(cur_integration_map, xflux_map, + yflux_map, zflux_map, + center_efield_, edge_efield_l_, + weight_l_, stale_depth); // Update longitudinal B-field (add source terms of constrained transport) for (int dim = 0; dim<3; dim++){ @@ -250,7 +250,7 @@ void EnzoBfieldMethodCT::update_all_bfield_components //---------------------------------------------------------------------- void EnzoBfieldMethodCT::compute_center_efield -(int dim, EFlt3DArray &efield, const EnzoEFltArrayMap &prim_map, +(const int dim, EFlt3DArray &efield, const EnzoEFltArrayMap &integration_map, int stale_depth) { const std::string v_names[3] = {"velocity_x", "velocity_y", "velocity_z"}; @@ -261,10 +261,10 @@ void EnzoBfieldMethodCT::compute_center_efield int k = coord.k_axis(); // Load the jth and kth components of the velocity and cell-centered bfield - EFlt3DArray velocity_j = prim_map.at(v_names[j]); - EFlt3DArray velocity_k = prim_map.at(v_names[k]); - EFlt3DArray bfield_j = prim_map.at(b_names[j]); - EFlt3DArray bfield_k = prim_map.at(b_names[k]); + EFlt3DArray velocity_j = integration_map.at(v_names[j]); + EFlt3DArray velocity_k = integration_map.at(v_names[k]); + EFlt3DArray bfield_j = integration_map.at(b_names[j]); + EFlt3DArray bfield_k = integration_map.at(b_names[k]); for (int iz=stale_depth; iz &edge_efield_l, std::array &weight_l, int stale_depth) { for (int i = 0; i < 3; i++){ - EnzoBfieldMethodCT::compute_center_efield(i, center_efield, prim_map, - stale_depth); + EnzoBfieldMethodCT::compute_center_efield(i, center_efield, + integration_map, stale_depth); EnzoEFltArrayMap *jflux_map; EnzoEFltArrayMap *kflux_map; diff --git a/src/Enzo/enzo_EnzoBfieldMethodCT.hpp b/src/Enzo/enzo_EnzoBfieldMethodCT.hpp index 8c4bc076e8..cd212e5c82 100644 --- a/src/Enzo/enzo_EnzoBfieldMethodCT.hpp +++ b/src/Enzo/enzo_EnzoBfieldMethodCT.hpp @@ -70,8 +70,8 @@ class EnzoBfieldMethodCT : public EnzoBfieldMethod /// Updates all components of the face-centered and the cell-centered bfields /// - /// @param[in] cur_prim_map Map containing the current values of the - /// cell-centered integrable quantities (before they have been updated + /// @param[in] cur_integration_map Map containing the current values of the + /// cell-centered integration quantities (before they have been updated /// over the current timestep). Specifically, the velocity and bfield /// components stored in this mapping are used to compute the /// cell-centered E-field. @@ -80,15 +80,15 @@ class EnzoBfieldMethodCT : public EnzoBfieldMethod /// namely makes use of the various magnetic field fluxes /// @param[out] out_centered_bfield_map Map holding the arrays where the /// updated values for each cell-centered magnetic field component should - /// be stored. These arrays can be aliases of `cur_prim_map` (The + /// be stored. These arrays can be aliases of `cur_integration_map` (The /// arguments can even reference the same object). /// @param[in] dt The (partial) time-step over which to update the magnetic /// fields. /// @param[in] stale_depth indicates the current stale_depth for the /// supplied quantities. This should nominally be the same as the stale /// depth used to compute the fluxes and that is passed to - /// EnzoIntegrableUpdate::update_quantities. - void update_all_bfield_components(EnzoEFltArrayMap &cur_prim_map, + /// EnzoIntegrationQuanUpdate::update_quantities. + void update_all_bfield_components(EnzoEFltArrayMap &cur_integration_map, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, @@ -133,14 +133,14 @@ class EnzoBfieldMethodCT : public EnzoBfieldMethod /// Values of 0, 1 and 2 correspond to the x, y and z directions. /// @param[out] center_efield The array where the computed cell-centered /// values of the E-field are written. - /// @param[in] prim_map Map containing the current values of the - /// cell-centered integrable quantities. Specifically, the velocity and + /// @param[in] integration_map Map containing the current values of the + /// cell-centered integration quantities. Specifically, the velocity and /// bfield entries are used to compute the cell-centered E-field. /// @param[in] stale_depth the stale depth at the time of this function call /// /// @note this function is called in `compute_all_edge_efields` static void compute_center_efield(int dim, EFlt3DArray ¢er_efield, - const EnzoEFltArrayMap &prim_map, + const EnzoEFltArrayMap &integration_map, int stale_depth = 0); /// Computes component i of the edge-centered E-field that sits on the faces @@ -180,10 +180,10 @@ class EnzoBfieldMethodCT : public EnzoBfieldMethod int stale_depth); /// Compute the all of the edge-centered electric fields using the current - /// fluxes and current cell-centered integrable quantities . + /// fluxes and current cell-centered integration quantities . /// - /// @param[in] cur_prim_map Map containing the current values of the - /// cell-centered integrable quantities (before they have been updated + /// @param[in] cur_integration_map Map containing the current values of the + /// cell-centered integration quantities (before they have been updated /// over the current timestep). Specifically, the velocity and bfield /// components stored in this mapping are used to compute the /// cell-centered E-field. @@ -208,7 +208,7 @@ class EnzoBfieldMethodCT : public EnzoBfieldMethod /// the weighting scheme used by Athena++ at a later date). /// @param[in] stale_depth the stale depth at the time of this function call static void compute_all_edge_efields - (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &xflux_map, + (EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, EFlt3DArray ¢er_efield, std::array &edge_efield_l, std::array &weight_l, int stale_depth); diff --git a/src/Enzo/enzo_EnzoCenteredFieldRegistry.cpp b/src/Enzo/enzo_EnzoCenteredFieldRegistry.cpp index 2ed50e7605..f8156005d4 100644 --- a/src/Enzo/enzo_EnzoCenteredFieldRegistry.cpp +++ b/src/Enzo/enzo_EnzoCenteredFieldRegistry.cpp @@ -29,10 +29,14 @@ const ft_map EnzoCenteredFieldRegistry::field_table_ = { //---------------------------------------------------------------------- std::vector EnzoCenteredFieldRegistry::get_registered_quantities -(bool enumerate_components) +(const bool enumerate_components, const bool exclude_non_active_advection) { std::vector out; for (auto const &key_item_pair : EnzoCenteredFieldRegistry::field_table_) { + if (exclude_non_active_advection && + (!(key_item_pair.second).actively_advected)){ + continue; + } if (enumerate_components && (key_item_pair.second).vector_quantity){ out.push_back(key_item_pair.first + std::string("_x")); out.push_back(key_item_pair.first + std::string("_y")); diff --git a/src/Enzo/enzo_EnzoCenteredFieldRegistry.hpp b/src/Enzo/enzo_EnzoCenteredFieldRegistry.hpp index e196d0f371..b9cacb15c1 100644 --- a/src/Enzo/enzo_EnzoCenteredFieldRegistry.hpp +++ b/src/Enzo/enzo_EnzoCenteredFieldRegistry.hpp @@ -149,8 +149,10 @@ class EnzoCenteredFieldRegistry /// When true, the function effectively returns a vector containing all /// registered fields. Otherwise, the returned vector containing the /// names of all quantities listed in FIELD_TABLE. + /// @param exclude_non_active_advection Determines whether quantities that + /// are bit actively advected should be excluded from the list static std::vector get_registered_quantities - (bool enumerate_components); + (const bool enumerate_components, const bool exclude_non_active_advection = false); /// provides the quantity properties listed in FIELD_TABLE (if present) /// diff --git a/src/Enzo/enzo_EnzoEFltArrayMap.hpp b/src/Enzo/enzo_EnzoEFltArrayMap.hpp index 1eba6b0d2e..5ff510acec 100644 --- a/src/Enzo/enzo_EnzoEFltArrayMap.hpp +++ b/src/Enzo/enzo_EnzoEFltArrayMap.hpp @@ -22,7 +22,7 @@ /// ordering is based off the requirements of the Riemann Solver and the class /// offers a method to return a (const?) vector of arrays with that ordering, /// optimizations could be made to the Riemann Solver (as well as -/// EnzoReconstructor and EnzoIntegrableUpdater). This could be most +/// EnzoReconstructor and EnzoIntegrationQuanUpdater). This could be most /// efficiently achieved by internally storing the contained array within a /// vector and also storing some kind of mapping relating the keys to indices /// the vector. diff --git a/src/Enzo/enzo_EnzoEOSIdeal.cpp b/src/Enzo/enzo_EnzoEOSIdeal.cpp index ac7b02cc53..54dfed6fdf 100644 --- a/src/Enzo/enzo_EnzoEOSIdeal.cpp +++ b/src/Enzo/enzo_EnzoEOSIdeal.cpp @@ -39,147 +39,86 @@ bool grackle_variable_gamma_(){ //---------------------------------------------------------------------- -// Helper function that performs a quick check to confirm that certain fields -// are held by both reconstructable_group and integrable_group -void confirm_same_kv_pair_(const EnzoEFltArrayMap &reconstructable, - const EnzoEFltArrayMap &integrable, - bool allowed_to_be_omitted, const std::string &key, - const std::string &func_name) +void EnzoEOSIdeal::primitive_from_integration +(EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &primitive_map, + const int stale_depth, const str_vec_t &passive_list) const { - if (allowed_to_be_omitted && ( (!integrable.contains(key)) || - (!reconstructable.contains(key)) ) ){ - return; - } else { - ASSERT1(func_name.c_str(), - ("The arrays associated with the \"%s\" key in the recontructable " - "and integrable maps are expected to be aliases of each other."), - key.c_str(), integrable.at(key).is_alias(reconstructable.at(key))); - } -} - -//---------------------------------------------------------------------- - -void check_recon_integ_overlap_ -(const EnzoEFltArrayMap &reconstructable, const EnzoEFltArrayMap &integrable, - const std::string &func_name, const str_vec_t &passive_list) -{ - // We assume that the following groups are represented by the same fields in - // integrable and reconstructable. This should probably not be hardcoded - const std::array standard_common_keys = - {"density", "velocity_x", "velocity_y", "velocity_z"}; - const std::array optional_common_keys = - {"bfield_x", "bfield_y", "bfield_z"}; - - for (const std::string& key : standard_common_keys){ - confirm_same_kv_pair_(reconstructable, integrable, false, key, func_name); - } - for (const std::string& key : optional_common_keys){ - confirm_same_kv_pair_(reconstructable, integrable, true, key, func_name); - } - for (const std::string& key : passive_list){ - confirm_same_kv_pair_(reconstructable, integrable, false, key, func_name); + if (grackle_variable_gamma_()){ + ERROR("EnzoEOSIdeal::primitive_from_integration", + "Doesn't currently support spatial variations in gamma"); } -} -//---------------------------------------------------------------------- - -void EnzoEOSIdeal::reconstructable_from_integrable -(EnzoEFltArrayMap &integrable, EnzoEFltArrayMap &reconstructable, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth, - const str_vec_t &passive_list) const -{ + EFlt3DArray density = integration_map.at("density"); + const int mz = density.shape(0); + const int my = density.shape(1); + const int mx = density.shape(2); - // Confirm that the expected fields (e.g. density, vx, vy, vz, bx, by, bz) - // are the same in reconstructable_group and integrable_group - check_recon_integ_overlap_(reconstructable, integrable, - "EnzoEOSIdeal::reconstructable_from_integrable", - passive_list); + // The EOS object doesn't necessarily know what the integration quantities + // are. This means we take something of an exhaustive approach. This could be + // more clever if this operations were made a part of the Hydro Integrator + std::vector quantity_list = + EnzoCenteredFieldRegistry::get_registered_quantities(true, true); - // Simply compute the pressure - pressure_from_integrable(integrable, reconstructable.at("pressure"), - conserved_passive_map, stale_depth); -} - -//---------------------------------------------------------------------- + for (const std::string& key : quantity_list){ + if ((!integration_map.contains(key)) || (!primitive_map.contains(key))){ + continue; + } -void EnzoEOSIdeal::integrable_from_reconstructable -(EnzoEFltArrayMap &reconstructable, EnzoEFltArrayMap &integrable, - int stale_depth, const str_vec_t &passive_list) const -{ - if (grackle_variable_gamma_()){ - // we would need to model a field with the - ERROR("EnzoEOSIdeal::pressure_from_integrable", - "Doesn't currently support spatial variations in gamma"); - } + EFlt3DArray integ_array = integration_map.at(key); + EFlt3DArray prim_array = primitive_map.at(key); - const bool idual = this->uses_dual_energy_formalism(); - const bool mag = (reconstructable.contains("bfield_x") || - reconstructable.contains("bfield_y") || - reconstructable.contains("bfield_z")); - // Confirm that the expected fields (e.g. density, vx, vy, vz, bx, by, bz) - // are the same in reconstructable_group and integrable_group - check_recon_integ_overlap_(reconstructable, integrable, - "EnzoEOSIdeal::integrable_from_reconstructable", - passive_list); - - EFlt3DArray density = reconstructable.get("density", stale_depth); - EFlt3DArray vx = reconstructable.get("velocity_x", stale_depth); - EFlt3DArray vy = reconstructable.get("velocity_y", stale_depth); - EFlt3DArray vz = reconstructable.get("velocity_z", stale_depth); - EFlt3DArray pressure = reconstructable.get("pressure", stale_depth); - - EFlt3DArray bx, by, bz; - if (mag){ - bx = reconstructable.get("bfield_x", stale_depth); - by = reconstructable.get("bfield_y", stale_depth); - bz = reconstructable.get("bfield_z", stale_depth); - } +#ifdef DEBUG_MATCHING_ARRAY_SHAPES + ASSERT6("EnzoEOSIdeal::primitive_from_integration", + ("The array being copied from integration_map has shape " + "(%d,%d,%d), while the destination array has shape (%d,%d,%d). " + "They should be the same."), + mz, my, mx, + prim_array.shape(0), prim_array.shape(1), prim_array.shape(2), + ((prim_array.shape(0) == mz) && + (prim_array.shape(1) == my) && + (prim_array.shape(2) == mx))); +#endif - EFlt3DArray eint, etot; - if (idual){ - eint = integrable.get("internal_energy", stale_depth); + for (int iz = stale_depth; iz < mz - stale_depth; iz++) { + for (int iy = stale_depth; iy < my - stale_depth; iy++) { + for (int ix = stale_depth; ix < mx - stale_depth; ix++) { + prim_array(iz,iy,ix) = integ_array(iz,iy,ix); + } + } + } } - etot = integrable.get("total_energy", stale_depth); - - enzo_float inv_gm1 = 1./(get_gamma()-1.); - for (int iz=0; izuses_dual_energy_formalism(); - const bool mag = (integrable_map.contains("bfield_x") || - integrable_map.contains("bfield_y") || - integrable_map.contains("bfield_z")); + const bool mag = (integration_map.contains("bfield_x") || + integration_map.contains("bfield_y") || + integration_map.contains("bfield_z")); // rather than slicing out the unstaled regions, we may want use the full // array and adjust the iteration limits accordingly. EFlt3DArray density, vx, vy, vz, eint, etot, bx, by, bz; - density = integrable_map.get("density", stale_depth); + density = integration_map.get("density", stale_depth); if (idual){ - eint = integrable_map.get("internal_energy", stale_depth); + eint = integration_map.get("internal_energy", stale_depth); } else { - etot = integrable_map.get("total_energy", stale_depth); - vx = integrable_map.get("velocity_x", stale_depth); - vy = integrable_map.get("velocity_y", stale_depth); - vz = integrable_map.get("velocity_z", stale_depth); + etot = integration_map.get("total_energy", stale_depth); + vx = integration_map.get("velocity_x", stale_depth); + vy = integration_map.get("velocity_y", stale_depth); + vz = integration_map.get("velocity_z", stale_depth); if (mag){ - bx = integrable_map.get("bfield_x", stale_depth); - by = integrable_map.get("bfield_y", stale_depth); - bz = integrable_map.get("bfield_z", stale_depth); + bx = integration_map.get("bfield_x", stale_depth); + by = integration_map.get("bfield_y", stale_depth); + bz = integration_map.get("bfield_z", stale_depth); } } @@ -247,52 +186,10 @@ void EnzoEOSIdeal::pressure_from_integrable //---------------------------------------------------------------------- -void EnzoEOSIdeal::pressure_from_reconstructable -(EnzoEFltArrayMap &reconstructable, EFlt3DArray &pressure, - int stale_depth) const -{ - // This is necessary since other equations of state may not include pressure - // as a reconstructable quantity. - - // We will check if pressure_name the field containing pressure in - // reconstructable_group are the same, if so then do nothing. Otherwise, - // simply copy the values over. - - if (reconstructable.contains("pressure") && - (reconstructable.at("pressure").is_alias(pressure))){ - // The arrays are aliases of each other - return; - } else { - EFlt3DArray old_p = reconstructable.get("pressure", stale_depth); - -#ifdef DEBUG_MATCHING_ARRAY_SHAPES - ASSERT6("EnzoEOSIdeal::pressure_from_reconstructable", - ("The pressure array in reconstructable has shape (%d,%d,%d) " - "while the array passed to this function has shape (%d,%d,%d). " - "They should be the same."), - old_p.shape(0), old_p.shape(1), old_p.shape(2), - pressure.shape(0), pressure.shape(1), pressure.shape(2), - ((old_p.shape(0) == pressure.shape(0)) && - (old_p.shape(1) == pressure.shape(1)) && - (old_p.shape(2) == pressure.shape(2)))); -#endif - - for (int iz=0; iz< old_p.shape(0); iz++) { - for (int iy=0; iy< old_p.shape(1); iy++) { - for (int ix=0; ix< old_p.shape(2); ix++) { - pressure(iz,iy,ix) = old_p(iz,iy,ix); - } - } - } - } -} - -//---------------------------------------------------------------------- - // based on the enzo's hydro_rk implementation of synchronization (found in the // Grid_UpdateMHD.C file) void EnzoEOSIdeal::apply_floor_to_energy_and_sync -(EnzoEFltArrayMap &integrable_map, int stale_depth) const +(EnzoEFltArrayMap &integration_map, const int stale_depth) const { if (grackle_variable_gamma_()){ ERROR("EnzoEOSIdeal::apply_floor_to_energy_and_sync", @@ -300,25 +197,25 @@ void EnzoEOSIdeal::apply_floor_to_energy_and_sync } const bool idual = this->uses_dual_energy_formalism(); - const bool mag = (integrable_map.contains("bfield_x") || - integrable_map.contains("bfield_y") || - integrable_map.contains("bfield_z")); + const bool mag = (integration_map.contains("bfield_x") || + integration_map.contains("bfield_y") || + integration_map.contains("bfield_z")); // in hydro_rk, eta was set equal to eta1 (it didn't use eta2 at all) const double eta = dual_energy_formalism_eta_; EFlt3DArray density, vx, vy, vz, etot, eint, bx, by, bz; - density = integrable_map.get("density", stale_depth); - vx = integrable_map.get("velocity_x", stale_depth); - vy = integrable_map.get("velocity_y", stale_depth); - vz = integrable_map.get("velocity_z", stale_depth); - etot = integrable_map.get("total_energy", stale_depth); + density = integration_map.get("density", stale_depth); + vx = integration_map.get("velocity_x", stale_depth); + vy = integration_map.get("velocity_y", stale_depth); + vz = integration_map.get("velocity_z", stale_depth); + etot = integration_map.get("total_energy", stale_depth); if (idual){ - eint = integrable_map.get("internal_energy", stale_depth); + eint = integration_map.get("internal_energy", stale_depth); } if (mag){ - bx = integrable_map.get("bfield_x", stale_depth); - by = integrable_map.get("bfield_y", stale_depth); - bz = integrable_map.get("bfield_z", stale_depth); + bx = integration_map.get("bfield_x", stale_depth); + by = integration_map.get("bfield_y", stale_depth); + bz = integration_map.get("bfield_z", stale_depth); } float ggm1 = get_gamma()*(get_gamma() - 1.); diff --git a/src/Enzo/enzo_EnzoEOSIdeal.hpp b/src/Enzo/enzo_EnzoEOSIdeal.hpp index 64b9b821d2..71da6a34b9 100644 --- a/src/Enzo/enzo_EnzoEOSIdeal.hpp +++ b/src/Enzo/enzo_EnzoEOSIdeal.hpp @@ -50,29 +50,20 @@ class EnzoEOSIdeal : public EnzoEquationOfState /// CHARM++ Pack / Unpack function void pup (PUP::er &p); - void reconstructable_from_integrable - (EnzoEFltArrayMap &integrable, EnzoEFltArrayMap &reconstructable, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth, - const str_vec_t &passive_list) const; + void primitive_from_integration + (EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &primitive_map, + const int stale_depth, const str_vec_t &passive_list) const; - void integrable_from_reconstructable - (EnzoEFltArrayMap &reconstructable, EnzoEFltArrayMap &integrable, - int stale_depth, const str_vec_t &passive_list) const; - - void pressure_from_integrable - (EnzoEFltArrayMap &integrable_map, const EFlt3DArray &pressure, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth) const; - - void pressure_from_reconstructable(EnzoEFltArrayMap &reconstructable, - EFlt3DArray &pressure, - int stale_depth) const; + void pressure_from_integration + (EnzoEFltArrayMap &integration_map, const EFlt3DArray &pressure, + const int stale_depth) const; inline enzo_float get_density_floor() const { return density_floor_; } enzo_float get_pressure_floor() const { return pressure_floor_; } - void apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integrable_map, - int stale_depth) const; + void apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integration_map, + const int stale_depth) const; bool is_barotropic() const { return false; } diff --git a/src/Enzo/enzo_EnzoEquationOfState.hpp b/src/Enzo/enzo_EnzoEquationOfState.hpp index f3f70e19af..7993159b92 100644 --- a/src/Enzo/enzo_EnzoEquationOfState.hpp +++ b/src/Enzo/enzo_EnzoEquationOfState.hpp @@ -88,101 +88,46 @@ class EnzoEquationOfState : public PUP::able PUP::able::pup(p); } - /// Converts cell-centered integrable primitives to reconstructable primitives + /// Converts integration quantities to primitives /// - /// @param[in] integrable Map holding integrable primitive values that are - /// to be converted - /// @param[out] reconstructable Map holding arrays where the computed - /// reconstructable data is to be stored. Due to the large degree of - /// overlap between integrable and reconstructable quantities, several - /// arrays held in this Map and `integrable` are expected to be aliases. - /// @param[in] conserved_passive_map Map containing the passively advected - /// scalars in conserved form (note that while `integrable` may also - /// contain the same passive scalars, those values are in specific form - - /// which are never used in this calculation). These are provided for - /// Grackle's use. + /// @param[in] integration_map Map holding integration quantities that are + /// to be converted. Passive scalars in this map are expected to be in + /// conserved form (they are densities). + /// @param[out] primitive_map Map holding arrays where the computed + /// primitive data is to be stored. Passive scalars in this map are + /// expected to be in specific form (they are mass fractions). /// @param[in] stale_depth indicates the current stale_depth for the /// supplied cell-centered quantities - /// @param[in] passive_list A list of keys for passive scalars. In this - /// method, it specifies which passive scalars in `reconstructable` and - /// `integrable` should be aliases of each other. + /// @param[in] passive_list A list of keys for passive scalars. These keys + /// will be used to determine which quantities will be copied from the + /// integration_map to the primitive_map. /// - /// For a barotropic EOS, this nominally does nothing - /// For a non-barotropic EOS, this computes pressure + /// Non-passive scalar quantities appearing in both `integration_map` and + /// `primitive_map` are simply deepcopied and passive scalar quantities are + /// converted from conserved-form to specific form. For a non-barotropic EOS, + /// this also computes pressure. /// /// @note - /// This interface is not ideal. `passive_list` is not really used for - /// anything except for checking aliasing. Now that this interface has - /// transitioned from using Groupings of field names to directly handling - /// arrays, it may be better to eliminate this method altogether. - virtual void reconstructable_from_integrable - (EnzoEFltArrayMap &integrable, EnzoEFltArrayMap &reconstructable, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth, - const str_vec_t &passive_list) const =0; - - /// Converts reconstructable primitives to integrable primitives - /// - /// @param[in] reconstructable Map holding reconstructable primitive values - /// that are to be converted - /// @param[out] integrable Map holding arrays where the computed integrable - /// data is to be stored. Due to the large degree of overlap between - /// integrable and reconstructable quantities, several arrays held in - /// this array and reconstructable are expected to be aliases. - /// @param[in] stale_depth indicates the current stale_depth for the - /// supplied cell-centered quantities - /// @param[in] passive_list A list of keys for passive scalars. In this - /// method, it specifies which passive scalars in `reconstructable` and - /// `integrable` should be aliases of each other. - /// - /// For a barotropic EOS, this nominally does nothing - /// For a non-barotropic EOS, this computes specific total energy from - /// pressure. If using the dual energy formalism, it will also compute the - /// internal energy from the pressure - /// - /// @note - /// This interface is not ideal. `passive_list` is not really used for - /// anything except for checking aliasing. Now that this interface has - /// transitioned from using Groupings of field names to directly handling - /// arrays, it may be better to eliminate this method altogether. - virtual void integrable_from_reconstructable - (EnzoEFltArrayMap &reconstructable, EnzoEFltArrayMap &integrable, - int stale_depth, const str_vec_t &passive_list) const =0; - - /// Computes thermal pressure from integrable quantities + /// It's a not obvious to me that the EOS should necessarily be responsible + /// for this operation. + virtual void primitive_from_integration + (EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &primitive_map, + const int stale_depth, const str_vec_t &passive_list) const =0; + + /// Computes thermal pressure from integration quantities /// - /// @param[in] integrable_map Map holding integrable primitive values - /// that are used to compute the pressure + /// @param[in] integration_map Map holding integration quantities that are + /// used to compute the pressure. This should include all necessary + /// passively advected quantities in conserved form. /// @param[out] pressure Array where the thermal pressure is to be stored - /// @param[in] conserved_passive_map Map containing the passively advected - /// scalars in conserved form (note that while `integrable` may also - /// contain the same passive scalars, those values are in specific form - - /// which are never used in this calculation). These are provided for - /// Grackle's use. /// @param[in] stale_depth indicates the current stale_depth for the /// supplied cell-centered quantities /// /// This nominally should wrap EnzoComputePressure. But at the time of /// writing, it doesn't actually wrap EnzoComputePressure - virtual void pressure_from_integrable - (EnzoEFltArrayMap &integrable_map, const EFlt3DArray &pressure, - EnzoEFltArrayMap &conserved_passive_map, int stale_depth) const = 0; - - /// Computes thermal pressure from reconstructable quantities (nominally - /// after reconstruction) - /// - /// @param[in] reconstructable Map holding reconstructable primitive values - /// that are used to compute the pressure - /// @param[out] pressure Array where the thermal pressure is to be stored - /// @param[in] stale_depth indicates the current stale_depth for the - /// supplied cell-centered quantities - /// - /// For a non-barotropic EOS, pressure is considered a reconstructable - /// quantity. In that case, if the pressure array in reconstructable is an - /// alias of the pressure array argument, nothing happens. If they aren't - /// aliases values are copied between arrays. - virtual void pressure_from_reconstructable(EnzoEFltArrayMap &reconstructable, - EFlt3DArray &pressure, - int stale_depth) const = 0; + virtual void pressure_from_integration + (EnzoEFltArrayMap &integration_map, const EFlt3DArray &pressure, + const int stale_depth) const = 0; /// returns the density floor virtual enzo_float get_density_floor() const = 0; @@ -195,9 +140,10 @@ class EnzoEquationOfState : public PUP::able /// and it synchronize the internal energy and total energy fields. If the /// EOS is barotropic, this does nothing. /// - /// @param[in,out] integrable_map Map holding integrable primitives that will - /// be used to apply the floor. It must also include a "total_energy" - /// entry (unless the EOS is barotropic) upon which the floor is applied. + /// @param[in,out] integration_map Map holding integration quantities that + /// will be used to apply the floor. It must also include a + /// "total_energy" entry (unless the EOS is barotropic) upon which the + /// floor is applied. /// @param[in] stale_depth indicates the current stale_depth for the /// supplied cell-centered quantities /// @@ -206,8 +152,8 @@ class EnzoEquationOfState : public PUP::able /// is a local operation that doesn't require data about neighboring cells /// (similar to the implementation of the dual energy formalsim in Enzo's /// Runge Kutta and MHD with Constrained Transport solvers). - virtual void apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integrable_map, - int stale_depth) const = 0; + virtual void apply_floor_to_energy_and_sync(EnzoEFltArrayMap &integration_map, + const int stale_depth) const = 0; /// returns whether the equation of state is barotropic virtual bool is_barotropic() const = 0; diff --git a/src/Enzo/enzo_EnzoIntegrableUpdate.cpp b/src/Enzo/enzo_EnzoIntegrationQuanUpdate.cpp similarity index 67% rename from src/Enzo/enzo_EnzoIntegrableUpdate.cpp rename to src/Enzo/enzo_EnzoIntegrationQuanUpdate.cpp index 259d4bca17..098f16ad26 100644 --- a/src/Enzo/enzo_EnzoIntegrableUpdate.cpp +++ b/src/Enzo/enzo_EnzoIntegrationQuanUpdate.cpp @@ -1,9 +1,9 @@ // See LICENSE_CELLO file for license and copyright information -/// @file enzo_EnzoIntegrableUpdate.cpp +/// @file enzo_EnzoIntegrationQuanUpdate.cpp /// @author Matthew Abruzzo (matthewabruzzo@gmail.com) /// @date Mon June 24 2019 -/// @brief [\ref Enzo] Implementation of EnzoIntegrableUpdate. +/// @brief [\ref Enzo] Implementation of EnzoIntegrationQuanUpdate. #include "cello.hpp" #include "enzo.hpp" @@ -12,10 +12,10 @@ //---------------------------------------------------------------------- static void append_key_to_vec_ -(const str_vec_t &integrable_quantities, FieldCat target_cat, - bool skip_bfield, std::size_t *density_index, str_vec_t &key_vec) +(const str_vec_t &integration_quantities, FieldCat target_cat, + const bool skip_bfield, std::size_t *density_index, str_vec_t &key_vec) { - for (std::string name : integrable_quantities){ + for (const std::string &name : integration_quantities){ bool vector_quantity, actively_advected; FieldCat category; bool success = EnzoCenteredFieldRegistry::quantity_properties @@ -26,14 +26,14 @@ static void append_key_to_vec_ ("\"%s\" is not registered in EnzoCenteredFieldRegistry"), name.c_str(), success); ASSERT1("append_key_to_vec_", - ("\"%s\" should not be listed as an integrable quantity because " + ("\"%s\" should not be listed as an integration quantity because " "it is not actively advected."), name.c_str(), actively_advected); if (category != target_cat){ ASSERT1("append_key_to_vec_", - ("Can't handle the integrable \"%s\" quantity because it has a " - "field category of FieldCat::other"), + ("Can't handle the integration quantity, \"%s\", because it has " + "a field category of FieldCat::other"), name.c_str(), category != FieldCat::other); continue; } else if (skip_bfield && (name == "bfield")){ @@ -54,29 +54,29 @@ static void append_key_to_vec_ //---------------------------------------------------------------------- -EnzoIntegrableUpdate::EnzoIntegrableUpdate -(const str_vec_t& integrable_quantities, - bool skip_B_update) throw() +EnzoIntegrationQuanUpdate::EnzoIntegrationQuanUpdate +(const str_vec_t& integration_quantities, + const bool skip_B_update) throw() { - // Add conserved quantities to integrable_keys_ and identify the index + // Add conserved quantities to integration_keys_ and identify the index // holding the density key density_index_ = std::numeric_limits::max(); - append_key_to_vec_(integrable_quantities, FieldCat::conserved, skip_B_update, - &density_index_, integrable_keys_); + append_key_to_vec_(integration_quantities, FieldCat::conserved, skip_B_update, + &density_index_, integration_keys_); // Confirm that density is in fact a registered quantity - ASSERT("EnzoIntegrableUpdate", - ("\"density\" must be a registered integrable quantity."), + ASSERT("EnzoIntegrationQuanUpdate", + ("\"density\" must be a registered integration quantity."), density_index_ != std::numeric_limits::max()); // Record the first index holding a key for a specific quantity - first_specific_index_ = integrable_keys_.size(); - // Add specific quantities to integrable_keys_ - append_key_to_vec_(integrable_quantities, FieldCat::specific, skip_B_update, - nullptr, integrable_keys_); + first_specific_index_ = integration_keys_.size(); + // Add specific quantities to integration_keys_ + append_key_to_vec_(integration_quantities, FieldCat::specific, skip_B_update, + nullptr, integration_keys_); } //---------------------------------------------------------------------- -void EnzoIntegrableUpdate::clear_dUcons_map +void EnzoIntegrationQuanUpdate::clear_dUcons_map (EnzoEFltArrayMap &dUcons_map, enzo_float value, const str_vec_t &passive_list) const noexcept { @@ -92,13 +92,13 @@ void EnzoIntegrableUpdate::clear_dUcons_map } }; - for (const std::string& key : integrable_keys_){ init_arr(key); } + for (const std::string& key : integration_keys_){ init_arr(key); } for (const std::string& key : passive_list){ init_arr(key); } } //---------------------------------------------------------------------- -void EnzoIntegrableUpdate::accumulate_flux_component +void EnzoIntegrationQuanUpdate::accumulate_flux_component (int dim, double dt, enzo_float cell_width, EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map, int stale_depth, const str_vec_t &passive_list) const noexcept @@ -137,46 +137,44 @@ void EnzoIntegrableUpdate::accumulate_flux_component }; - for (const std::string& key : integrable_keys_){ accumulate(key); } + for (const std::string& key : integration_keys_){ accumulate(key); } for (const std::string& key : passive_list){ accumulate(key); } } //---------------------------------------------------------------------- -EFlt3DArray* EnzoIntegrableUpdate::load_integrable_quantities_ +EFlt3DArray* EnzoIntegrationQuanUpdate::load_integration_quantities_ (EnzoEFltArrayMap &map, int stale_depth) const { - std::size_t nfields = integrable_keys_.size(); + std::size_t nfields = integration_keys_.size(); EFlt3DArray* arr = new EFlt3DArray[nfields]; for (std::size_t i = 0; iget_density_floor(); EFlt3DArray *cur_prim, *dU, *out_prim; - cur_prim = load_integrable_quantities_(initial_integrable_map, stale_depth); - dU = load_integrable_quantities_(dUcons_map, stale_depth); - out_prim = load_integrable_quantities_(out_integrable_map,stale_depth); - std::size_t nfields = integrable_keys_.size(); + cur_prim = load_integration_quantities_(initial_integration_map, stale_depth); + dU = load_integration_quantities_(dUcons_map, stale_depth); + out_prim = load_integration_quantities_(out_integration_map,stale_depth); + std::size_t nfields = integration_keys_.size(); for (int iz = 1; iz < (cur_prim[density_index_].shape(0) - 1); iz++) { for (int iy = 1; iy < (cur_prim[density_index_].shape(1) - 1); iy++) { @@ -185,7 +183,7 @@ void EnzoIntegrableUpdate::update_quantities // get the initial density enzo_float old_rho = cur_prim[density_index_](iz,iy,ix); - // now update the integrable primitives that are conserved + // now update the integration quantities that are conserved for (std::size_t i = 0; i < first_specific_index_; i++){ out_prim[i](iz,iy,ix) = cur_prim[i](iz,iy,ix) + dU[i](iz,iy,ix); } @@ -195,7 +193,7 @@ void EnzoIntegrableUpdate::update_quantities new_rho = EnzoEquationOfState::apply_floor(new_rho, density_floor); out_prim[density_index_](iz,iy,ix) = new_rho; - // update the specific integrable primitives + // update the specific integration primitives enzo_float inv_new_rho = 1./new_rho; for (std::size_t i = first_specific_index_; i < nfields; i++){ out_prim[i](iz,iy,ix) = @@ -208,35 +206,33 @@ void EnzoIntegrableUpdate::update_quantities // apply floor to energy and sync the internal energy with total energy // (the latter only occurs if the dual energy formalism is in use) - eos->apply_floor_to_energy_and_sync(out_integrable_map, stale_depth + 1); + eos->apply_floor_to_energy_and_sync(out_integration_map, stale_depth + 1); delete[] cur_prim; delete[] dU; delete[] out_prim; } //---------------------------------------------------------------------- -void EnzoIntegrableUpdate::update_passive_scalars_ -(EnzoEFltArrayMap &initial_integrable_map, EnzoEFltArrayMap &dUcons_map, - EnzoEFltArrayMap &out_conserved_passive_scalar, int stale_depth, +void EnzoIntegrationQuanUpdate::update_passive_scalars_ +(EnzoEFltArrayMap &initial_integration_map, EnzoEFltArrayMap &dUcons_map, + EnzoEFltArrayMap &out_integration_map, const int stale_depth, const str_vec_t &passive_list) const { - EFlt3DArray cur_rho = initial_integrable_map.get("density", stale_depth); - // cell-centered grid dimensions + EFlt3DArray cur_rho = initial_integration_map.get("density", stale_depth); int mz, my, mx; mz = cur_rho.shape(0); my = cur_rho.shape(1); mx = cur_rho.shape(2); for (const std::string &key : passive_list){ - EFlt3DArray cur_specific, out_conserved, dU; - cur_specific = initial_integrable_map.get(key, stale_depth); - out_conserved = out_conserved_passive_scalar.get(key, stale_depth); + EFlt3DArray cur_conserved, out_conserved, dU; + cur_conserved = initial_integration_map.get(key, stale_depth); + out_conserved = out_integration_map.get(key, stale_depth); dU = dUcons_map.get(key, stale_depth); for (int iz=1; iz -#ifndef ENZO_ENZO_INTEGRABLE_UPDATE_HPP -#define ENZO_ENZO_INTEGRABLE_UPDATE_HPP -class EnzoIntegrableUpdate : public PUP::able +#ifndef ENZO_ENZO_INTEGRATION_QUAN_UPDATE_HPP +#define ENZO_ENZO_INTEGRATION_QUAN_UPDATE_HPP +class EnzoIntegrationQuanUpdate : public PUP::able { - /// @class EnzoIntegrableUpdate + /// @class EnzoIntegrationQuanUpdate /// @ingroup Enzo /// @brief [\ref Enzo] Encapsulates the updating of (actively and - /// passively) advected integrable quantites + /// passively) advected integration quantites public: // interface - /// Create a new EnzoIntegrableUpdate instance - EnzoIntegrableUpdate(const std::vector &integrable_quantities, - bool skip_B_update) throw(); + /// Create a new EnzoIntegrationQuanUpdate instance + EnzoIntegrationQuanUpdate + (const std::vector &integration_quantities, + const bool skip_B_update) throw(); /// Virtual destructor - virtual ~EnzoIntegrableUpdate() + virtual ~EnzoIntegrationQuanUpdate() { } /// CHARM++ PUP::able declaration - PUPable_decl(EnzoIntegrableUpdate); + PUPable_decl(EnzoIntegrationQuanUpdate); /// CHARM++ migration constructor for PUP::able - EnzoIntegrableUpdate (CkMigrateMessage *m) + EnzoIntegrationQuanUpdate (CkMigrateMessage *m) : PUP::able(m) { } @@ -39,25 +41,25 @@ class EnzoIntegrableUpdate : public PUP::able void pup (PUP::er &p) { PUP::able::pup(p); - p|integrable_keys_; + p|integration_keys_; p|first_specific_index_; p|density_index_; } - /// Iterates through all arrays in `dUcons_map` that are registered as - /// integrable quantities or are specified with `passive_list`. All elements + /// Iterates through all arrays in `dUcons_map` that are pre-registered + /// integration quantities or are specified with `passive_list`. All elements /// in these arrays are set to `value`. /// /// @param[in,out] dUcons_map The map holding the arrays that are to be /// modified. These arrays are nominally used to accumulate the changes - /// to all integrable and passively advected quantites. + /// to all integration quantities. /// @param[in] value The value array elements are assigned. /// @param[in] passive_list A list of keys for passive scalars. void clear_dUcons_map(EnzoEFltArrayMap &dUcons_map, enzo_float value, const str_vec_t &passive_list) const noexcept; - /// Computes the change in (the conserved form of) the integrable and - /// passively advected quantites due to the flux divergence along dimension + /// Computes the change in (the conserved form of) the integration quantities + /// (including passive scalars) due to the flux divergence along dimension /// `dim` over the timestep `dt`. These changes are added to the accumulator /// arrays contained by `dUcons_map`. /// @@ -80,39 +82,35 @@ class EnzoIntegrableUpdate : public PUP::able EnzoEFltArrayMap &dUcons_map, int stale_depth, const str_vec_t &passive_list) const noexcept; - /// adds flux divergence (and source terms) to the initial integrable - /// quantities and stores the results in `out_integrable_map` + /// adds flux divergence (and source terms) to the initial integration + /// quantities and stores the results in `out_integration_map` /// - /// @param[in] initial_integrable_map Map of arrays holding the values of - /// integrable quantities from the start of the timestep. The fluxes + /// @param[in] initial_integration_map Map of arrays holding the values of + /// integration quantities from the start of the timestep. The fluxes /// will be added to fields held in this map. This should also contain - /// the passive scalars in specific form (as mass fractions) that are to - /// be integrated. + /// the passive scalars in conserved form (as densities) that are to be + /// integrated. /// @param[in] dUcons_map Map of arrays where the net changes to the - /// integrable quantities and passively advected quantites are stored. + /// integration quantities and passively advected quantites are stored. /// If constrained transport is being used, this will not include arrays /// for the magnetic fields. - /// @param[out] out_integrable_map Map of the fields where the updated - /// integrable quantities will be stored (This can be a reference to the - /// same Map referenced by initial_integrable_map). The updated passively - /// advected scalars will NOT be stored here. - /// @param[out] out_conserved_passive_scalar Map of arrays where the updated - /// passive scalar quantities are stored in conserved form (as densities). + /// @param[out] out_integration_map Map of the fields where the updated + /// integration quantities will be stored (This can be a reference to the + /// same Map referenced by initial_integration_map). The updated + /// passively advected scalars will be stored here. /// @param[in] eos Pointer to the fluid's equation of state object. When /// applicable used for placing a density floor. /// @param[in] stale_depth The stale depth at the time of the function call /// (the stale_depth must be incremented after this function is called) /// @param[in] passive_list A list of keys for passive scalars. void update_quantities - (EnzoEFltArrayMap &initial_integrable_map, EnzoEFltArrayMap &dUcons_map, - EnzoEFltArrayMap &out_integrable_map, - EnzoEFltArrayMap &out_conserved_passive_scalar, - EnzoEquationOfState *eos, int stale_depth, - const str_vec_t &passive_list) const; + (EnzoEFltArrayMap &initial_integration_map, EnzoEFltArrayMap &dUcons_map, + EnzoEFltArrayMap &out_integration_map, EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t &passive_list) const; - /// provides a const vector of all registerred integrable keys - const std::vector integrable_keys() const throw() - { return integrable_keys_; } + /// provides a const vector of all registerred integration keys + const std::vector integration_keys() const throw() + { return integration_keys_; } private: @@ -122,29 +120,29 @@ class EnzoIntegrableUpdate : public PUP::able /// /// (This should called before the density is updated) void update_passive_scalars_ - (EnzoEFltArrayMap &initial_integrable_map, EnzoEFltArrayMap &dUcons_map, - EnzoEFltArrayMap &out_conserved_passive_scalar, int stale_depth, + (EnzoEFltArrayMap &initial_integration_map, EnzoEFltArrayMap &dUcons_map, + EnzoEFltArrayMap &out_integration_map, const int stale_depth, const str_vec_t &passive_list) const; /// Dynamically allocates and constructs an array of instances of EFlt3DArray - /// that are loaded from `map` using the ordering of keys in integrable_keys_ - /// @param[in] map Map of arrays holding data related to each integrable - /// quantities registered in integrable_keys_ + /// that are loaded from `map` using the ordering of keys in integration_keys_ + /// @param[in] map Map of arrays holding data related to each integration + /// quantities registered in integration_keys_ /// @param[in] stale_depth indicates the current stale_depth for the loaded /// quantities. - EFlt3DArray* load_integrable_quantities_(EnzoEFltArrayMap &map, - int stale_depth) const; + EFlt3DArray* load_integration_quantities_(EnzoEFltArrayMap &map, + const int stale_depth) const; private: // attributes - /// Holds key names used to load each integrable quantity component from a + /// Holds key names used to load each integration quantity component from a /// a mapping. Keys for conserved quantities are always listed before the /// specfic quantities. This excludes passively advected scalars. - std::vector integrable_keys_; - /// index of the first specific quantity listed in integrable_keys_ + std::vector integration_keys_; + /// index of the first specific quantity listed in integration_keys_ std::size_t first_specific_index_; - /// index of integrable_keys_ that holds the string "density" + /// index of integration_keys_ that holds the string "density" std::size_t density_index_; }; -#endif /* ENZO_ENZO_INTEGRABLE_UPDATE_HPP */ +#endif /* ENZO_ENZO_INTEGRATION_QUAN_UPDATE_HPP */ diff --git a/src/Enzo/enzo_EnzoLazyPassiveScalarFieldList.cpp b/src/Enzo/enzo_EnzoLazyPassiveScalarFieldList.cpp index f8be864435..a50578b0f1 100644 --- a/src/Enzo/enzo_EnzoLazyPassiveScalarFieldList.cpp +++ b/src/Enzo/enzo_EnzoLazyPassiveScalarFieldList.cpp @@ -35,13 +35,15 @@ EnzoLazyPassiveScalarFieldList::build_passive_list_() noexcept void EnzoLazyPassiveScalarFieldList::pup(PUP::er &p) { p | initialized_; - if (p.isUnpacking()){ - str_vec_t* temp_vals = new std::vector(); - p | *temp_vals; - passive_names_ = std::const_pointer_cast - (std::shared_ptr(temp_vals)); - } else { - str_vec_t temp_copy = *passive_names_; - p | temp_copy; + if (initialized_){ + if (p.isUnpacking()){ + str_vec_t* temp_vals = new std::vector(); + p | *temp_vals; + passive_names_ = std::const_pointer_cast + (std::shared_ptr(temp_vals)); + } else { + str_vec_t temp_copy = *passive_names_; + p | temp_copy; + } } } diff --git a/src/Enzo/enzo_EnzoMethodMHDVlct.cpp b/src/Enzo/enzo_EnzoMethodMHDVlct.cpp index 0d5d36ec1f..24aacae0ce 100644 --- a/src/Enzo/enzo_EnzoMethodMHDVlct.cpp +++ b/src/Enzo/enzo_EnzoMethodMHDVlct.cpp @@ -59,29 +59,42 @@ EnzoMethodMHDVlct::EnzoMethodMHDVlct (std::string rsolver, eos_ = new EnzoEOSIdeal(gamma, density_floor, pressure_floor, dual_energy_formalism, dual_energy_formalism_eta); +#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 */ + // Determine whether magnetic fields are to be used mhd_choice_ = parse_bfield_choice_(mhd_choice); - // determine integrable and reconstructable quantities - std::vector integrable_quantities, reconstructable_quantities; - EnzoMethodMHDVlct::determine_quantities_(eos_, mhd_choice_, - integrable_quantities, - reconstructable_quantities); + 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; + integration_quantities = riemann_solver_->integration_quantities(); + primitive_quantities = riemann_solver_->primitive_quantities(); // Initialize the remaining component objects half_dt_recon_ = EnzoReconstructor::construct_reconstructor - (reconstructable_quantities, half_recon_name, (enzo_float)theta_limiter); + (primitive_quantities, half_recon_name, (enzo_float)theta_limiter); full_dt_recon_ = EnzoReconstructor::construct_reconstructor - (reconstructable_quantities, full_recon_name, (enzo_float)theta_limiter); - riemann_solver_ = EnzoRiemann::construct_riemann(integrable_quantities, - rsolver); - integrable_updater_ = new EnzoIntegrableUpdate(integrable_quantities, - true); + (primitive_quantities, full_recon_name, (enzo_float)theta_limiter); + + integration_quan_updater_ = + new EnzoIntegrationQuanUpdate(integration_quantities, true); - // Determine the list of integrable fields and reconstructable fields and - // then make sure all required fields are defined - build_field_l_(integrable_quantities, integrable_field_list_); - build_field_l_(reconstructable_quantities, reconstructable_field_list_); + // Determine the lists of fields that are required to hold the integration + // quantities and primitives and ensure that they are defined + build_field_l_(integration_quantities, integration_field_list_); + build_field_l_(primitive_quantities, primitive_field_list_); // make sure "pressure" is defined (it's needed to compute the timestep) FieldDescr * field_descr = cello::field_descr(); @@ -131,58 +144,13 @@ EnzoMethodMHDVlct::bfield_choice EnzoMethodMHDVlct::parse_bfield_choice_ //---------------------------------------------------------------------- -void EnzoMethodMHDVlct::determine_quantities_ -(const EnzoEquationOfState *eos, EnzoMethodMHDVlct::bfield_choice mhd_choice, - std::vector &integrable_quantities, - std::vector &reconstructable_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){ - integrable_quantities.push_back(quantity); - reconstructable_quantities.push_back(quantity); - } - - if (mhd_choice != bfield_choice::no_bfield){ - integrable_quantities.push_back("bfield"); - reconstructable_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 integrable quantities - integrable_quantities.push_back("total_energy"); - // add pressure to reconstructable quantities - reconstructable_quantities.push_back("pressure"); - // add specific internal energy to integrable quantities (if using the dual - // engery formalism) - if (eos->uses_dual_energy_formalism()){ - integrable_quantities.push_back("internal_energy"); - } - } -} - -//---------------------------------------------------------------------- - EnzoMethodMHDVlct::~EnzoMethodMHDVlct() { delete eos_; delete half_dt_recon_; delete full_dt_recon_; delete riemann_solver_; - delete integrable_updater_; + delete integration_quan_updater_; if (bfield_method_ != nullptr){ delete bfield_method_; } @@ -202,69 +170,35 @@ void EnzoMethodMHDVlct::pup (PUP::er &p) p|half_dt_recon_; p|full_dt_recon_; p|riemann_solver_; - p|integrable_updater_; + p|integration_quan_updater_; p|mhd_choice_; p|bfield_method_; - p|integrable_field_list_; - p|reconstructable_field_list_; + p|integration_field_list_; + p|primitive_field_list_; p|lazy_passive_list_; } //---------------------------------------------------------------------- -// Returns the unique members of a combination of 2 vectors -std::vector unique_combination_(const std::vector &a, - const std::vector &b) +EnzoEFltArrayMap EnzoMethodMHDVlct::get_integration_map_ +(Block * block, const str_vec_t *passive_list) const noexcept { - std::vector out; - // copy elements from a - for (std::size_t i = 0; i all_prim_fields = - unique_combination_(integrable_field_list_, reconstructable_field_list_); + EnzoEFltArrayMap map("integration"); EnzoFieldArrayFactory array_factory(block,0); - for (const std::string& field_name : all_prim_fields){ - ASSERT1("EnzoMethodMHDVlct::nonpassive_primitive_map_", + for (const std::string& field_name : integration_field_list_){ + ASSERT1("EnzoMethodMHDVlct::get_integration_map_", "EnzoEFltArrayMap can't hold more than one key called \"%s\"", field_name.c_str(), !map.contains(field_name)); map[field_name] = array_factory.from_name(field_name); } - return map; -} - -//---------------------------------------------------------------------- -EnzoEFltArrayMap EnzoMethodMHDVlct::conserved_passive_scalar_map_ -(Block * block) const noexcept -{ - EnzoEFltArrayMap map("conserved_passive_scalar"); - std::shared_ptr passive_list= lazy_passive_list_.get_list(); - EnzoFieldArrayFactory array_factory(block,0); - for (const std::string& field_name : (*passive_list)){ - ASSERT1("EnzoMethodMHDVlct::conserved_passive_scalar_map_", - "EnzoEFltArrayMap can't hold more than one key called \"%s\"", - field_name.c_str(), !map.contains(field_name)); - map[field_name] = array_factory.from_name(field_name); + if (passive_list != nullptr){ + for (const std::string& field_name : (*passive_list)){ + ASSERT1("EnzoMethodMHDVlct::get_integration_map_", + "EnzoEFltArrayMap can't hold more than one key called \"%s\"", + field_name.c_str(), !map.contains(field_name)); + map[field_name] = array_factory.from_name(field_name); + } } return map; } @@ -281,14 +215,19 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() // fields and/or serve as scratch space. // map that holds arrays wrapping the Cello Fields holding each of the - // primitive quantities. Additionally, this also includes temporary arrays - // used to hold the specific form of the passive scalar - EnzoEFltArrayMap primitive_map; // this will be overwritten + // integration quantities. Additionally, this also includes temporary + // arrays used to hold the specific form of the passive scalar + EnzoEFltArrayMap integration_map; // this will be overwritten - // map used for storing primitive values at the half time-step. This - // includes key,array pairs for each entry in primitive_map (there should + // map used for storing integration values at the half time-step. This + // includes key,array pairs for each entry in integration_map (there should // be no aliased fields shared between maps) - EnzoEFltArrayMap temp_primitive_map("temp_primitive"); + EnzoEFltArrayMap temp_integration_map("temp_integration"); + + // Map of arrays used to temporarily store the cell-centered primitive + // quantities that are subsequently reconstructed. This includes arrays for + // storing the specific form of each of the passively advected scalars. + EnzoEFltArrayMap primitive_map("primitive"); // map holding the arrays wrapping each fields corresponding to a passively // advected scalar. The scalar is in conserved form. @@ -298,14 +237,6 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() EnzoEFltArrayMap priml_map("priml"); EnzoEFltArrayMap primr_map("primr"); - // Arrays used to store the pressure computed from the reconstructed left - // and right primitives - // Note: in the case of adiabatic fluids, pressure is a reconstructable - // quantity and entries are included for it in priml_map and - // primr_map. In that case, pressure_l and pressure_r are aliases of - // those arrays. - EFlt3DArray pressure_l, pressure_r; - // maps used to store fluxes (in the future, these will wrap FluxData // entries) EnzoEFltArrayMap xflux_map("xflux"); @@ -313,15 +244,14 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() EnzoEFltArrayMap zflux_map("zflux"); // map of arrays used to accumulate the changes to the conserved forms of - // the integrable quantities and passively advected scalars. In other + // the integration quantities and passively advected scalars. In other // words, at the start of the (partial) timestep, the fields are all set to // zero and are used to accumulate the flux divergence and source terms. If // CT is used, it won't have space to store changes in the magnetic fields. EnzoEFltArrayMap dUcons_map("dUcons"); - setup_arrays_(block, primitive_map, temp_primitive_map, - conserved_passive_scalar_map, priml_map, primr_map, - pressure_l, pressure_r, xflux_map, yflux_map, zflux_map, + setup_arrays_(block, integration_map, temp_integration_map, primitive_map, + priml_map, primr_map, xflux_map, yflux_map, zflux_map, dUcons_map); // Setup a pointer to an array that used to store interface velocity fields @@ -330,7 +260,7 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() // use, don't actually allocate the array and set the pointer to NULL. EFlt3DArray interface_velocity_arr, *interface_velocity_arr_ptr; if (eos_->uses_dual_energy_formalism()){ - EFlt3DArray density = primitive_map.at("density"); + EFlt3DArray density = integration_map.at("density"); interface_velocity_arr = EFlt3DArray(density.shape(0), density.shape(1), density.shape(2)); interface_velocity_arr_ptr = &interface_velocity_arr; @@ -352,27 +282,13 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() // refreshed) extends over. int stale_depth = 0; - // convert the passive scalars from conserved form to specific form - // (outside the integrator, they are treated like conserved densities) - compute_specific_passive_scalars_(*(lazy_passive_list_.get_list()), - primitive_map["density"], - conserved_passive_scalar_map, - primitive_map, stale_depth); - // repeat the following loop twice (for half time-step and full time-step) - for (int i=0;i<2;i++){ double cur_dt = (i == 0) ? dt/2. : dt; - EnzoEFltArrayMap& cur_integrable_map = - (i == 0) ? primitive_map : temp_primitive_map; - EnzoEFltArrayMap& out_integrable_map = - (i == 0) ? temp_primitive_map : primitive_map; - // For the purposes of making the calculation procedure slightly more - // more explicit, we distinguish between cur_integrable_group and - // cur_reconstructable_group. Due to the high level of overlap between - // these, they are simply aliases of the same underlying Grouping that - // holds groups for both of them - EnzoEFltArrayMap& cur_reconstructable_map = cur_integrable_map; + EnzoEFltArrayMap& cur_integration_map = + (i == 0) ? integration_map : temp_integration_map; + EnzoEFltArrayMap& out_integration_map = + (i == 0) ? temp_integration_map : integration_map; EnzoReconstructor *reconstructor; @@ -380,56 +296,36 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() reconstructor = half_dt_recon_; } else { reconstructor = full_dt_recon_; - - // After the fluxes were added to the passive scalar in the first half - // timestep, the values were stored in conserved form in the fields - // held by conserved_passive_scalar_map. - // Need to convert them to specific form - compute_specific_passive_scalars_(*(lazy_passive_list_.get_list()), - cur_integrable_map["density"], - conserved_passive_scalar_map, - cur_integrable_map, stale_depth); } // set all elements of the arrays in dUcons_map to 0 (throughout the rest // of the current loop, flux divergence and source terms will be // accumulated in these arrays) - integrable_updater_->clear_dUcons_map(dUcons_map, 0., - *(lazy_passive_list_.get_list())); - - // Compute the reconstructable quantities from the integrable quantites - // Although cur_integrable_map holds the passive scalars in integrable - // form, the conserved form of the values is required in case Grackle is - // being used. - // - // Note: cur_integrable_map and cur_reconstructable_map are aliases of - // the same map since there is such a large degree of overlap between - // reconstructable and integrable quantities - // - // For a barotropic gas, the following nominally does nothing - // For a non-barotropic gas, the following nominally computes pressure - eos_->reconstructable_from_integrable(cur_integrable_map, - cur_reconstructable_map, - conserved_passive_scalar_map, - stale_depth, - *(lazy_passive_list_.get_list())); + integration_quan_updater_->clear_dUcons_map + (dUcons_map, 0., *(lazy_passive_list_.get_list())); + + // Compute the primitive quantities from the integration quantites + // This basically copies all quantities that are both and an integration + // quantity and a primitive and converts the passsive scalars from + // conserved-form to specific-form (i.e. from density to mass fraction). + // For a non-barotropic gas, this also computes pressure + eos_->primitive_from_integration(cur_integration_map, primitive_map, + stale_depth, + *(lazy_passive_list_.get_list())); // Compute flux along each dimension - compute_flux_(0, cur_dt, cell_widths[0], cur_reconstructable_map, - priml_map, primr_map, pressure_l, pressure_r, - xflux_map, dUcons_map, interface_velocity_arr_ptr, - *reconstructor, bfield_method_, stale_depth, - *(lazy_passive_list_.get_list())); - compute_flux_(1, cur_dt, cell_widths[1], cur_reconstructable_map, - priml_map, primr_map, pressure_l, pressure_r, - yflux_map, dUcons_map, interface_velocity_arr_ptr, - *reconstructor, bfield_method_, stale_depth, - *(lazy_passive_list_.get_list())); - compute_flux_(2, cur_dt, cell_widths[2], cur_reconstructable_map, - priml_map, primr_map, pressure_l, pressure_r, - zflux_map, dUcons_map, interface_velocity_arr_ptr, - *reconstructor, bfield_method_, stale_depth, - *(lazy_passive_list_.get_list())); + compute_flux_(0, cur_dt, cell_widths[0], primitive_map, + priml_map, primr_map, xflux_map, dUcons_map, + interface_velocity_arr_ptr, *reconstructor, bfield_method_, + stale_depth, *(lazy_passive_list_.get_list())); + compute_flux_(1, cur_dt, cell_widths[1], primitive_map, + priml_map, primr_map, yflux_map, dUcons_map, + interface_velocity_arr_ptr, *reconstructor, bfield_method_, + stale_depth, *(lazy_passive_list_.get_list())); + compute_flux_(2, cur_dt, cell_widths[2], primitive_map, + priml_map, primr_map, zflux_map, dUcons_map, + interface_velocity_arr_ptr, *reconstructor, bfield_method_, + stale_depth, *(lazy_passive_list_.get_list())); // increment the stale_depth stale_depth+=reconstructor->immediate_staling_rate(); @@ -438,25 +334,22 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw() // Update Bfields if (bfield_method_ != nullptr) { - bfield_method_->update_all_bfield_components(cur_integrable_map, + bfield_method_->update_all_bfield_components(cur_integration_map, xflux_map, yflux_map, zflux_map, - out_integrable_map, cur_dt, + out_integration_map, + cur_dt, stale_depth); bfield_method_->increment_partial_timestep(); } - // Update quantities (includes flux divergence and source terms) - // This needs to happen after updating the cell-centered B-field so that - // the pressure floor can be applied to the total energy (and if - // necessary the total energy can be synchronized with internal energy) - // - // Note: updated passive scalars are NOT saved in out_integrable_group in - // specific form. Instead they are saved in conserved_passive_scalars - // in conserved form. - integrable_updater_->update_quantities - (primitive_map, dUcons_map, out_integrable_map, - conserved_passive_scalar_map, eos_, stale_depth, + // Update the integration quantities (includes flux divergence and source + // terms). This currently needs to happen after updating the + // cell-centered B-field so that the pressure floor can be applied to the + // total energy (and if necessary the total energy can be synchronized + // with the internal energy) + integration_quan_updater_->update_quantities + (integration_map, dUcons_map, out_integration_map, eos_, stale_depth, *(lazy_passive_list_.get_list())); // increment stale_depth since the inner values have been updated @@ -493,53 +386,18 @@ void EnzoMethodMHDVlct::check_mesh_and_ghost_size_(Block *block) const noexcept //---------------------------------------------------------------------- -void EnzoMethodMHDVlct::compute_specific_passive_scalars_ -(const str_vec_t &passive_list, EFlt3DArray& density, - EnzoEFltArrayMap& conserved_passive_scalar_map, - EnzoEFltArrayMap& specific_passive_scalar_map, int stale_depth) const noexcept -{ - int mz = density.shape(0); - int my = density.shape(1); - int mx = density.shape(2); - - for (const std::string& key : passive_list){ - EFlt3DArray cur_conserved = conserved_passive_scalar_map.at(key); - EFlt3DArray out_specific = specific_passive_scalar_map.at(key); - - for (int iz = stale_depth; iz < mz - stale_depth; iz++) { - for (int iy = stale_depth; iy < my - stale_depth; iy++) { - for (int ix = stale_depth; ix < mx - stale_depth; ix++) { - out_specific(iz,iy,ix) = cur_conserved(iz,iy,ix)/density(iz,iy,ix); - } - } - } - - } -} - -//---------------------------------------------------------------------- - void EnzoMethodMHDVlct::compute_flux_ -(int dim, double cur_dt, enzo_float cell_width, - EnzoEFltArrayMap &reconstructable_map, +(const int dim, const double cur_dt, const enzo_float cell_width, + EnzoEFltArrayMap &primitive_map, EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map, - EFlt3DArray &pressure_l, EFlt3DArray &pressure_r, EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map, EFlt3DArray *interface_velocity_arr_ptr, EnzoReconstructor &reconstructor, - EnzoBfieldMethod *bfield_method, int stale_depth, + EnzoBfieldMethod *bfield_method, const int stale_depth, const str_vec_t& passive_list) const noexcept { - // purely for the purposes of making the caluclation more explicit, we define - // the following aliases for priml_map and primr_map - EnzoEFltArrayMap &reconstructable_l = priml_map; - EnzoEFltArrayMap &reconstructable_r = primr_map; - EnzoEFltArrayMap &integrable_l = priml_map; - EnzoEFltArrayMap &integrable_r = primr_map; - // First, reconstruct the left and right interface values - reconstructor.reconstruct_interface(reconstructable_map, - reconstructable_l, reconstructable_r, + reconstructor.reconstruct_interface(primitive_map, priml_map, primr_map, dim, eos_, stale_depth, passive_list); // We temporarily increment the stale_depth for the rest of this calculation @@ -549,45 +407,36 @@ void EnzoMethodMHDVlct::compute_flux_ // Need to set the component of reconstructed B-field along dim, equal to // the corresponding longitudinal component of the B-field tracked at cell - // interfaces (should potentially be handled internally by reconstructor) + // interfaces if (bfield_method != nullptr) { - bfield_method->correct_reconstructed_bfield(reconstructable_l, - reconstructable_r, + bfield_method->correct_reconstructed_bfield(priml_map, primr_map, dim, cur_stale_depth); } - // Calculate integrable values on left and right faces: - eos_->integrable_from_reconstructable(reconstructable_l, integrable_l, - cur_stale_depth, passive_list); - eos_->integrable_from_reconstructable(reconstructable_r, integrable_r, - cur_stale_depth, passive_list); - - // Calculate pressure on left and right faces: - eos_->pressure_from_reconstructable(reconstructable_l, pressure_l, - cur_stale_depth); - eos_->pressure_from_reconstructable(reconstructable_r, pressure_r, - cur_stale_depth); - // Next, compute the fluxes - riemann_solver_->solve(integrable_l, integrable_r, pressure_l, pressure_r, - flux_map, dim, eos_, cur_stale_depth, passive_list, + riemann_solver_->solve(priml_map, primr_map, flux_map, dim, eos_, + cur_stale_depth, passive_list, interface_velocity_arr_ptr); - // Accumulate the change in integrable quantities from these flux_map in - // dUcons_map - integrable_updater_->accumulate_flux_component(dim, cur_dt, cell_width, - flux_map, dUcons_map, - cur_stale_depth, passive_list); - - // if using dual energy formalism, compute the dual energy formalism + // Compute the changes in the conserved form of the integration quantities + // from the fluxes and use these values to update dUcons_map (which is used + // to accumulate the total change in these quantities over the current + // [partial] timestep) + integration_quan_updater_->accumulate_flux_component(dim, cur_dt, cell_width, + flux_map, dUcons_map, + cur_stale_depth, + passive_list); + + // if using dual energy formalism, compute the component of the internal + // energy source term for this dim (and update dUcons_map). if (eos_->uses_dual_energy_formalism()){ EnzoSourceInternalEnergy eint_src; - eint_src.calculate_source(dim, cur_dt, cell_width, reconstructable_map, + eint_src.calculate_source(dim, cur_dt, cell_width, primitive_map, dUcons_map, *interface_velocity_arr_ptr, eos_, cur_stale_depth); } - // Finally, have the record the upwind direction (for handling CT) + // Finally, have bfield_method record the upwind direction (for handling CT) if (bfield_method != nullptr){ bfield_method->identify_upwind(flux_map, dim, cur_stale_depth); } @@ -596,7 +445,7 @@ void EnzoMethodMHDVlct::compute_flux_ //---------------------------------------------------------------------- static void add_temporary_arrays_to_map_ -(EnzoEFltArrayMap &map, std::array &shape, +(EnzoEFltArrayMap &map, const std::array &shape, const std::vector* const nonpassive_names, const str_vec_t* const passive_lists) { @@ -617,11 +466,9 @@ static void add_temporary_arrays_to_map_ //---------------------------------------------------------------------- void EnzoMethodMHDVlct::setup_arrays_ -(Block *block, EnzoEFltArrayMap &primitive_map, - EnzoEFltArrayMap &temp_primitive_map, - EnzoEFltArrayMap &conserved_passive_scalar_map, +(Block *block, EnzoEFltArrayMap &integration_map, + EnzoEFltArrayMap &temp_integration_map, EnzoEFltArrayMap &primitive_map, EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map, - EFlt3DArray &pressure_l, EFlt3DArray &pressure_r, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, EnzoEFltArrayMap &dUcons_map) noexcept { @@ -629,68 +476,53 @@ void EnzoMethodMHDVlct::setup_arrays_ // allocate stuff! Make sure to do it in a way such that we don't have to // separately deallocate it! - // To assist with setting up arrays, let's create a list of ALL primitive - // keys (including passive scalars) and all integrable keys. These are the - // same thing except the latter excludes quantities (like pressure) that we - // don't compute pressure for. - std::vector combined_key_list = unique_combination_ - (integrable_field_list_, reconstructable_field_list_); - - // First, setup conserved_passive_scalar_map - conserved_passive_scalar_map = conserved_passive_scalar_map_(block); - - // Next, setup nonpassive components of primitive_map - primitive_map = nonpassive_primitive_map_(block); - std::array shape = {primitive_map.at("density").shape(0), - primitive_map.at("density").shape(1), - primitive_map.at("density").shape(2)}; - add_temporary_arrays_to_map_(primitive_map, shape, nullptr, - (lazy_passive_list_.get_list()).get()); + // First, setup integration_map + integration_map = get_integration_map_ + (block, (lazy_passive_list_.get_list()).get()); + + const std::array shape = {integration_map.at("density").shape(0), + integration_map.at("density").shape(1), + integration_map.at("density").shape(2)}; - // Then, setup temp_primitive_map - add_temporary_arrays_to_map_(temp_primitive_map, shape, &combined_key_list, + // Next, setup temp_integration_map + add_temporary_arrays_to_map_(temp_integration_map, shape, + &integration_field_list_, (lazy_passive_list_.get_list()).get()); - // Prepare arrays to hold fluxes (it should include groups for all actively - // and passively advected quantities) + // Prepare arrays to hold fluxes. It should include keys for all integration + // quantities actively (including passively advected scalars) EnzoEFltArrayMap* flux_maps[3] = {&zflux_map, &yflux_map, &xflux_map}; for (std::size_t i = 0; i < 3; i++){ std::array cur_shape = shape; // makes a deep copy cur_shape[i] -= 1; add_temporary_arrays_to_map_(*(flux_maps[i]), cur_shape, - &combined_key_list, + &integration_field_list_, (lazy_passive_list_.get_list()).get()); } - // Prepare fields used to accumulate all changes to the actively advected and - // passively advected quantities. If CT is in use, dUcons_group should not - // have storage for magnetic fields since CT independently updates magnetic - // fields (this exclusion is implicitly handled integrable_updater_) - std::vector tmp = integrable_updater_->integrable_keys(); + // Prepare fields used to accumulate all changes to the integration + // quantities (including passively advected scalars). If CT is in use, + // dUcons_group should not have storage for magnetic fields since CT + // independently updates magnetic fields (this exclusion is implicitly + // handled by integration_quan_updater_) + std::vector tmp = integration_quan_updater_->integration_keys(); add_temporary_arrays_to_map_(dUcons_map, shape, &tmp, (lazy_passive_list_.get_list()).get()); - // Prepare temporary fields for priml and primr + // Setup primitive_map + add_temporary_arrays_to_map_(primitive_map, shape, + &primitive_field_list_, + (lazy_passive_list_.get_list()).get()); + + // Prepare maps for holding the left and right reconstructed primitives. // As necessary, we pretend that these are centered along: // - z and have shape (mz-1, my, mx) // - y and have shape ( mz,my-1, mx) // - x and have shape ( mz, my,mx-1) - add_temporary_arrays_to_map_(priml_map, shape, &combined_key_list, + add_temporary_arrays_to_map_(priml_map, shape, &primitive_field_list_, (lazy_passive_list_.get_list()).get()); - add_temporary_arrays_to_map_(primr_map, shape, &combined_key_list, + add_temporary_arrays_to_map_(primr_map, shape, &primitive_field_list_, (lazy_passive_list_.get_list()).get()); - - // If there are pressure entries in priml_map and primr_map (depends on the - // EOS), set pressure_l and pressure_name_r equal to - // those field names. Otherwise, reserve/allocate left/right pressure fields - if (priml_map.contains("pressure")) { - pressure_l = priml_map.at("pressure"); - pressure_r = primr_map.at("pressure"); - } else { - pressure_l = EFlt3DArray(shape[0], shape[1], shape[2]); - pressure_r = EFlt3DArray(shape[0], shape[1], shape[2]); - } - } //---------------------------------------------------------------------- @@ -700,39 +532,37 @@ double EnzoMethodMHDVlct::timestep ( Block * block ) const throw() // analogous to ppm timestep calulation, probably want to require that cfast // is no smaller than some tiny positive number. - // Construct a map holding the field data for each of the (non-passive) - // primitive quantities. - EnzoEFltArrayMap primitive_map = nonpassive_primitive_map_(block); + // Constructs a map containing the field data for each integration quantity + // This includes each passively advected scalar (as densities) + EnzoEFltArrayMap integration_map = get_integration_map_ + (block, (lazy_passive_list_.get_list()).get()); if (eos_->uses_dual_energy_formalism()){ // synchronize eint and etot. // This is only strictly necessary after problem initialization and when // there is an inflow boundary condition - eos_->apply_floor_to_energy_and_sync(primitive_map, 0); + eos_->apply_floor_to_energy_and_sync(integration_map, 0); } // Compute thermal pressure (this presently requires that "pressure" is a // permanent field) EnzoFieldArrayFactory array_factory(block); EFlt3DArray pressure = array_factory.from_name("pressure"); - EnzoEFltArrayMap conserved_passive_scalar_map = - conserved_passive_scalar_map_(block); - eos_->pressure_from_integrable(primitive_map, pressure, - conserved_passive_scalar_map, 0); + eos_->pressure_from_integration(integration_map, pressure, 0); // Now load other necessary quantities - enzo_float gamma = eos_->get_gamma(); - EFlt3DArray density = primitive_map.at("density"); - EFlt3DArray velocity_x = primitive_map.at("velocity_x"); - EFlt3DArray velocity_y = primitive_map.at("velocity_y"); - EFlt3DArray velocity_z = primitive_map.at("velocity_z"); + const enzo_float gamma = eos_->get_gamma(); + EFlt3DArray density = integration_map.at("density"); + EFlt3DArray velocity_x = integration_map.at("velocity_x"); + EFlt3DArray velocity_y = integration_map.at("velocity_y"); + EFlt3DArray velocity_z = integration_map.at("velocity_z"); const bool mhd = (mhd_choice_ != bfield_choice::no_bfield); EFlt3DArray bfieldc_x, bfieldc_y, bfieldc_z; if (mhd) { - bfieldc_x = primitive_map.at("bfield_x"); - bfieldc_y = primitive_map.at("bfield_y"); - bfieldc_z = primitive_map.at("bfield_z"); + bfieldc_x = integration_map.at("bfield_x"); + bfieldc_y = integration_map.at("bfield_y"); + bfieldc_z = integration_map.at("bfield_z"); } // widths of cells diff --git a/src/Enzo/enzo_EnzoMethodMHDVlct.hpp b/src/Enzo/enzo_EnzoMethodMHDVlct.hpp index 2d86c2634b..3896385902 100644 --- a/src/Enzo/enzo_EnzoMethodMHDVlct.hpp +++ b/src/Enzo/enzo_EnzoMethodMHDVlct.hpp @@ -11,7 +11,7 @@ /// EnzoEquationOfState: Equation of State for Gas /// EnzoReconstructor: Reconstructs primitive variables /// EnzoRiemann: Solves the Riemann Problem -/// EnzoIntegrableUpdate: handles the updating of advected quantites +/// EnzoIntegrationQuanUpdate: Handles updates to integration quantites /// /// This class can be run with and without magnetic fields. When run with /// magnetic fields, it makes use of the following component: @@ -19,17 +19,20 @@ /// /// Some notes on implementation /// - this Method tracks specific total energy (referred to as total_energy) -/// - We categorize quantities as reconstructable and integrable primitives. -/// Nearly every field overlaps between the two cases. Examples of the -/// primitives for adiabatic, ideal gas (without dual energy formalism): -/// - density (both integrable and reconstructable) -/// - velocity (both integrable and reconstructable) -/// - pressure (just reconstructable) -/// - total_energy (just integrable) +/// - We categorize quantities as integration quantities quantities and +/// primitives. Nearly every field overlaps between the two cases. Each +/// integration quantity is either a conserved quantity or a conserved +/// quanity divided by density (i.e. specific) Examples of the integration +/// quantities and primitives for adiabatic, ideal gas (without dual +/// energy formalism) include: +/// - density (both integration and primitive) +/// - velocity (both integration and primitive) +/// - pressure (just primitive) +/// - total_energy (just integration) /// When supporting magnetic fields, there is also: -/// - bfield (both integrable and reconstructable) +/// - bfield (both integration and primitive) /// When using the dual energy formalism there is also: -/// - internal_energy (just integrable) +/// - internal_energy (just integration) /// /// EnzoEFltArrayMap Objects /// ------------------------ @@ -38,26 +41,27 @@ /// class implements a map/dictionary that holds instances of EFlt3DArray. /// All arrays in a given map are assumed to have the same shape. /// -/// Currently, 8 different maps are used: -/// 1. primitive_map: holds the integrable and reconstructed quantities -/// that are stored in Cello fields. It also holds the values of all -/// passively advected scalars, in specific form, which are stored -/// in temporary arrays -/// 2. temp_primive_map: This holds temporary arrays for each of the -/// quantities in primitive_map. These arrays are used to store the -/// estimated state at the partial timestep. -/// 3. priml_map: holds left reconstructed primitive fields (has the +/// Currently, 9 different maps are used: +/// 1. integration_map: Map of arrays wrapping the Cello Fields holding +/// each of the integration quantities. This includes each of the +/// passive scalars (as densities). +/// 2. temp_integration_map: Map of arrays used to hold temporary values +/// of for each of the quantities in integration_map. These arrays +/// are used to store the estimated values at the partial timestep. +/// 3. primitive_map: Map of arrays used to temporarily store the +/// primitive quantities +/// 4. priml_map: holds left reconstructed primitive quantities (has the /// same keys as primitive_map) -/// 4. primr_map: holds right reconstructed primitive fields (has the +/// 5. primr_map: holds right reconstructed primitive quantities (has the /// same keys as primitive_map) -/// 5. xflux_map: holds fluxes in the x-direction -/// 6. yflux_map: holds fluxes in the y-direction -/// 7. zflux_map: holds fluxes in the z-direction -/// 8. dUcons_map: holds arrays which are used to accumulate the total -/// change in the conserved versions of all integrable quantities +/// 6. xflux_map: holds fluxes in the x-direction +/// 7. yflux_map: holds fluxes in the y-direction +/// 8. zflux_map: holds fluxes in the z-direction +/// 9. dUcons_map: holds arrays which are used to accumulate the total +/// change in the conserved versions of all integration quantities /// (includes both flux divergence and source terms) /// Note: -/// - All arrays in maps 1, 2, and 8 have the same shapes as +/// - All arrays in maps 1, 2, 3, and 8 have the same shapes as /// cell-centered Cello Fields /// - All arrays in priml_map and primr_map technically have the shape /// of a cell-centered field, but they are treated as though they have @@ -66,11 +70,6 @@ /// - For the purposes of these enumerated maps, we assume that the /// length of a face-centered array along the dimension with /// face-centering is 1 less than that of a cell-centered array -/// - All reconstructed fields are not technically registered as -/// cell-centered fields. Consequently, they are formally registered -/// as cell-centered fields (to guarantee that they have enough space). -/// - All reconstructable and integrable primitive quantities have -/// key-array pairs named for them in maps number 1, 2, 3, and 4. #ifndef ENZO_ENZO_METHOD_VLCT_HPP #define ENZO_ENZO_METHOD_VLCT_HPP @@ -118,11 +117,11 @@ class EnzoMethodMHDVlct : public Method { half_dt_recon_(nullptr), full_dt_recon_(nullptr), riemann_solver_(nullptr), - integrable_updater_(nullptr), + integration_quan_updater_(nullptr), mhd_choice_(bfield_choice::no_bfield), bfield_method_(nullptr), - integrable_field_list_(), - reconstructable_field_list_(), + integration_field_list_(), + primitive_field_list_(), lazy_passive_list_() { } @@ -146,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 quantities from (FIELD_TABLE) to be reconstructed and - /// integrated and that will be integrated. - /// - /// @param[in] eos Pointer to the fluid's EquationOfState object. - /// @param[in] mhd_choice Encodes how the integrator will handle B-fields - /// @param[out] integrable_quantities Reference to a vector that get's filled - /// by this function with the integrable quantities (matching names in - /// FIELD_TABLE) used by the integrator - /// @param[out] reconstructable_quantities Reference to a vector that get's - /// filled by this function with the names of 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 &integrable_quantities, - std::vector &reconstructable_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 @@ -170,32 +152,15 @@ class EnzoMethodMHDVlct : public Method { /// @param[in] block used to determine the current mesh size and ghost depth void check_mesh_and_ghost_size_(Block *block) const noexcept; - /// Converts conservative passive scalars (which are originally densities) - /// to specific form (basically just divide by density) - /// - /// @param[in] passive_list A list of keys for passive scalars. - /// @param[in] density Array holding the current density values - /// @param[in] conserved_passive_scalar_map Map of the arrays containing the - /// current values of the passively advected scalars (in conserved form) - /// @param[out] specific_passive_scalar_map Map of arrays where the specific - /// form of the scalars will be stored. - /// @param[in] stale_depth The current stale depth - void compute_specific_passive_scalars_ - (const str_vec_t &passive_list, EFlt3DArray& density, - EnzoEFltArrayMap& conserved_passive_scalar_map, - EnzoEFltArrayMap& specific_passive_scalar_map, - int stale_depth) const noexcept; - - /// Constructs a map containing the field data for each primitive (except for - /// the passively advected scalars). - EnzoEFltArrayMap nonpassive_primitive_map_(Block * block) const noexcept; - - /// Constructs a map containing the field data for each passively advected - /// scalar (in conserved form). - EnzoEFltArrayMap conserved_passive_scalar_map_(Block * block) const noexcept; + /// Constructs a map containing the field data for each integration quantity + /// This includes all passively advected scalars (as densities) included in + /// passive_list + EnzoEFltArrayMap get_integration_map_(Block * block, + const str_vec_t *passive_list) + const noexcept; /// Computes the fluxes along a given dimension, `dim`, and accumulate the - /// changes to the integrable quantities in `dUcons_map` + /// changes to the integration quantities in `dUcons_map` /// /// If using the dual energy formalism, this also computes a part of the /// internal energy density source term, @@ -215,26 +180,22 @@ class EnzoMethodMHDVlct : public Method { /// 1, and 2 correspond to the x, y, and z directions, respectively. /// @param[in] dt The current timestep. /// @param[in] cell_width The cell width along dimension `dim`. - /// @param[in] reconstructable_map Map of arrays holding cell-centered + /// @param[in] primitive_map Map of arrays holding cell-centered /// primitive quantities that are to be reconstructed (This includes /// specific passive scalars). /// @param[in] priml_map,primr_map Maps of arrays used to temporarily - /// hold the left/right reconstructed face-centered reconstructable and - /// integrable quantities. These arrays should have the shape of a - /// cell-centered field, but are treated as though they have the shape of - /// a that is face-centered along `dim` (If a cell-centered field holds - /// `N` elements along `dim`, then such a face-centered field should only - /// have `N-1` elements along `dim`). - /// @param[in] pressure_l,pressure_r Arrays used to temporarily store the - /// left/right pressure values. The shape of these arrays should be the - /// same as in priml_map,primr_map. Note, pressure_l (pressure_r) is - /// allowed to be a shallow copy of an array in priml_map (primr_map). + /// hold the left/right reconstructed face-centered primitives. These + /// arrays should have the shape of a cell-centered field, but are + /// treated as though they have the shape of a that is face-centered + /// along `dim` (If a cell-centered field holds `N` elements along `dim`, + /// then such a face-centered field should only have `N-1` elements along + /// `dim`). /// @param[in] flux_map Holds arrays where the calculated fluxes /// will be stored. The arrays should be face-centered along `dim`. /// If a cell-centered field holds `N` elements along `dim`, then this /// should only hold `N-1` elements along `dim`. /// @param[in,out] dUcons_map Map of arrays where the changes to the - /// integrable quantities are accumulated. If constrained transport is + /// integration quantities are accumulated. If constrained transport is /// being used, this won't include arrays for the magnetic fields. /// @param[in] interface_velocity_arr_ptr Pointer to an array to /// temporarily hold the computed component of the velocity at the cell @@ -244,7 +205,7 @@ class EnzoMethodMHDVlct : public Method { /// under the dual energy formalism). If the value is `nullptr`, then the /// interface velocity is not stored in the array. /// @param[in] reconstructor the instance of EnzoReconstructor to use to - /// update reconstruct the face centered values + /// update reconstruct the face-centered primitives /// @param[in,out] bfield_method When using running with bfield handling, this /// is a pointer to an instance of EnzoBfieldMethod. During the function /// call, the internal state is updated. If not handling bfields, this @@ -252,60 +213,49 @@ class EnzoMethodMHDVlct : public Method { /// @param[in] stale_depth indicates the current stale depth (before /// performing reconstruction) /// @param[in] passive_list A list of keys for passively advected scalars. - /// - /// @par Note - /// It might be worth breaking this into 2 functions (where one of them - /// handles the reconstruction of the fields and calculation of the flux and - /// the other additionally handles the accumulation of values in flux_map - /// and calculates any relevant source terms. void compute_flux_ - (int dim, double cur_dt, enzo_float cell_width, - EnzoEFltArrayMap &reconstructable_map, + (const int dim, const double cur_dt, const enzo_float cell_width, + EnzoEFltArrayMap &primitive_map, EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map, - EFlt3DArray &pressure_l, EFlt3DArray &pressure_r, EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map, EFlt3DArray *interface_velocity_arr_ptr, EnzoReconstructor &reconstructor, - EnzoBfieldMethod *bfield_method, int stale_depth, + EnzoBfieldMethod *bfield_method, const int stale_depth, const str_vec_t& passive_list) const noexcept; /// Setup arrays used throughout `compute`. This includes both arrays that /// wrap Cello fields AND temporary arrays used as scratch space. /// /// @param[in] block holds data to be processed - /// @param[out] primitive_map Map of arrays wrapping the Cello Fields holding - /// each of the integrable and reconstructable quantity. This also holds - /// temporary arrays where the specific form of the passively advected - /// scalars will be stored. - /// @param[out] temp_primitive_map Map for storing the integrable and - /// reconstructable quantities at the half timestep. This should have all + /// @param[out] integration_map Map of arrays wrapping the Cello Fields + /// holding each of the integration quantities. This includes each of the + /// passive scalars (as densities). + /// @param[out] temp_integration_map Map for storing the integration + /// quantities at the half timestep. This should have all /// the same entries as primitive_map. However, all arrays in this map /// are temporary. + /// @param[out] primitive_map Map of arrays used to temporarily store the + /// cell-centered primitive quantities that are subsequently + /// reconstructed. This includes arrays for storing the specific form of + /// each of the passively advected scalars. /// @param[out] priml_map,primr_map Maps of arrays used to temporarily - /// hold the left/right reconstructed face-centered reconstructable and - /// integrable quantities. These arrays should have the shape of a - /// cell-centered field so that they can be reused for multiple - /// dimensions. These have the same keys as primitive_map. - /// @param[out] pressure_l,pressure_r Arrays used to temporarily store the - /// left/right pressure values. The shape of these arrays should be the - /// same as in priml_map,primr_map. Note, for adiabatic equations of - /// state pressure_l (pressure_r) is a shallow copy of an array in - /// priml_map (primr_map). + /// hold the left/right reconstructed face-centered primitive quantities. + /// These arrays should have the shape of a cell-centered field so that + /// they can be reused for multiple dimensions. These have the same keys + /// as primitive_map. /// @param[out] xflux_map, yflux_map, zflux_map Maps of temporary arrays that /// are used to store the x, y, and z fluxes. A given map of arrays will /// hold values along at the face-centers along the direction of the /// fluxes. Note, if a cell-centered field holds `N` elements along /// `dim`, then this should only hold `N-1` elements along `dim`. /// @param[out] dUcons_map Map of temporary arrays used to accumulate the - /// changes to the conserved forms of the integrable quantities and + /// changes to the conserved forms of the integration quantities and /// passively advected scalars. If CT is used, this grouping won't have /// space to store changes in the magnetic fields (that update is handled /// separately). void setup_arrays_ - (Block *block, EnzoEFltArrayMap &primitive_map, - EnzoEFltArrayMap &temp_primitive_map, - EnzoEFltArrayMap &conserved_passive_scalar_map, + (Block *block, EnzoEFltArrayMap &integration_map, + EnzoEFltArrayMap &temp_integration_map, EnzoEFltArrayMap &primitive_map, EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map, - EFlt3DArray &pressure_l, EFlt3DArray &pressure_r, EnzoEFltArrayMap &xflux_map, EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map, EnzoEFltArrayMap &dUcons_map) noexcept; @@ -321,8 +271,8 @@ class EnzoMethodMHDVlct : public Method { EnzoReconstructor *full_dt_recon_; /// Pointer to the Riemann solver EnzoRiemann *riemann_solver_; - /// Pointer to the integrable quantity updater - EnzoIntegrableUpdate *integrable_updater_; + /// Pointer to the integration quantity updater + EnzoIntegrationQuanUpdate *integration_quan_updater_; /// Indicates how magnetic fields are handled bfield_choice mhd_choice_; @@ -330,13 +280,14 @@ class EnzoMethodMHDVlct : public Method { /// Pointer to the BfieldMethod handler EnzoBfieldMethod *bfield_method_; - /// Names of the integrable fields (only includes the field names for + /// Names of the integration fields (only includes the field names for /// actively advected quantities). These also serve as the keys to the /// mappings of arrays used in the calculation - std::vector integrable_field_list_; - /// Names of the reconstructable primitive fields. These also serve as the - /// keys to the mappings of arrays used in the calculation - std::vector reconstructable_field_list_; + std::vector integration_field_list_; + /// Names of the primitive fields (this should exclude passively advected + /// scalars). These also serve as the keys to the mappings of arrays used in + /// the calculation + std::vector primitive_field_list_; /// Lazy initializer of the list of fields holding passive scalars EnzoLazyPassiveScalarFieldList lazy_passive_list_; diff --git a/src/Enzo/enzo_EnzoReconstructor.hpp b/src/Enzo/enzo_EnzoReconstructor.hpp index c68a15e915..82fa76cfdc 100644 --- a/src/Enzo/enzo_EnzoReconstructor.hpp +++ b/src/Enzo/enzo_EnzoReconstructor.hpp @@ -61,9 +61,10 @@ class EnzoReconstructor : public PUP::able /// Reconstructs the interface values /// /// @param[in] prim_map Map holding the data for the cell-centered - /// reconstructable primitives. This is expected to have keys for all - /// of the active reconstructed quantities registered with the factory - /// method (plus all of the keys listed in `passive lists`) + /// primitives that are to be reconstructed. This is expected to have + /// a key-array pair for each entry in the list passed as the + /// ``active_reconstructed_quantities`` argument of the factory method + /// (plus all of the keys listed in `passive lists`) /// @param[out] priml_map,primr_map Holds existing arrays where the /// left/right reconstructed, face-centered primitives are written. /// These must supply the same keys that are expected for prim_map. @@ -85,8 +86,8 @@ class EnzoReconstructor : public PUP::able /// for other axes. virtual void reconstruct_interface (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &priml_map, - EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t& passive_list)=0; + EnzoEFltArrayMap &primr_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t& passive_list)=0; /// The rate amount by which the stale_depth increases after the current /// reconstructor is used to update the fluid over a (partial or full) @@ -125,8 +126,8 @@ class EnzoReconstructor : public PUP::able { return total_staling_rate() - immediate_staling_rate(); } protected: - /// list of the key names for all components of (non-passively advected - /// quantities) that are to be reconstructed. + /// list of the key names for all components of (non-passively advected) + /// primitives that are to be reconstructed. std::vector active_key_names_; }; diff --git a/src/Enzo/enzo_EnzoReconstructorNN.cpp b/src/Enzo/enzo_EnzoReconstructorNN.cpp index e814d538c2..df03e252cb 100644 --- a/src/Enzo/enzo_EnzoReconstructorNN.cpp +++ b/src/Enzo/enzo_EnzoReconstructorNN.cpp @@ -12,8 +12,8 @@ void EnzoReconstructorNN::reconstruct_interface (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &priml_map, - EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t& passive_list) + EnzoEFltArrayMap &primr_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t& passive_list) { // determine components of i unit vector diff --git a/src/Enzo/enzo_EnzoReconstructorNN.hpp b/src/Enzo/enzo_EnzoReconstructorNN.hpp index 17a937afca..3410b95bda 100644 --- a/src/Enzo/enzo_EnzoReconstructorNN.hpp +++ b/src/Enzo/enzo_EnzoReconstructorNN.hpp @@ -39,8 +39,8 @@ class EnzoReconstructorNN : public EnzoReconstructor void reconstruct_interface (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &priml_map, - EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t& passive_list); + EnzoEFltArrayMap &primr_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t& passive_list); int total_staling_rate() { return 1; } diff --git a/src/Enzo/enzo_EnzoReconstructorPLM.hpp b/src/Enzo/enzo_EnzoReconstructorPLM.hpp index 8114be1ada..4f1a109ba4 100644 --- a/src/Enzo/enzo_EnzoReconstructorPLM.hpp +++ b/src/Enzo/enzo_EnzoReconstructorPLM.hpp @@ -165,8 +165,8 @@ class EnzoReconstructorPLM : public EnzoReconstructor void reconstruct_interface (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &priml_map, - EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t& passive_list); + EnzoEFltArrayMap &primr_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t& passive_list); int total_staling_rate() { return 2; } @@ -184,8 +184,8 @@ class EnzoReconstructorPLM : public EnzoReconstructor template void EnzoReconstructorPLM::reconstruct_interface (EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &priml_map, - EnzoEFltArrayMap &primr_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t& passive_list) + EnzoEFltArrayMap &primr_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t& passive_list) { EnzoPermutedCoordinates coord(dim); Limiter limiter_func = Limiter(); @@ -193,8 +193,8 @@ void EnzoReconstructorPLM::reconstruct_interface auto fn = [coord, limiter_func, theta_limiter, stale_depth, &prim_map, &priml_map, &primr_map](const std::string &key, - bool use_floor, - enzo_float prim_floor) + const bool use_floor, + const enzo_float prim_floor) { // Cast the problem as reconstructing values at: // wl(k, j, i+3/2) and wr(k,j,i+1/2) diff --git a/src/Enzo/enzo_EnzoRiemann.cpp b/src/Enzo/enzo_EnzoRiemann.cpp index 1263f0b089..77b55d44da 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(const std::string& solver, const bool mhd, + const 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 8a0f90504f..c4265b5c3b 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 integrable quantities (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(const std::string& solver, const bool mhd, + const bool internal_energy); - EnzoRiemann() throw() + EnzoRiemann() noexcept {} /// Virtual destructor @@ -54,15 +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 integrable 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[in] pressure_array_l,pressure_array_r Arrays holding the - /// precomputed left/right reconstructed pressure values. These should - /// have the same face-centering as the arrays in priml_map/primr_map. - /// @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 @@ -73,7 +74,7 @@ class EnzoRiemann : public PUP::able /// `dim` (the array should not include exterior faces of the block and /// should be cell-centered along other dimensions). This quantity is /// used to compute the internal energy source term (needed under the - /// dual energy formalism). If the value is `NULL`, then the interface + /// dual energy formalism). If the value is `nullptr`, then the interface /// velocity is not stored in the array. /// /// @note It's alright for arrays in `priml_map` and `primr_map` to have the @@ -84,11 +85,19 @@ class EnzoRiemann : public PUP::able /// applies to the other arrays passed as arguments. virtual 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 str_vec_t &passive_list, - EFlt3DArray *interface_velocity) const = 0; + EnzoEFltArrayMap &flux_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t &passive_list, + const 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 b91897bc0b..271c4d988b 100644 --- a/src/Enzo/enzo_EnzoRiemannImpl.hpp +++ b/src/Enzo/enzo_EnzoRiemannImpl.hpp @@ -74,15 +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 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 @@ -151,24 +151,14 @@ class EnzoRiemannImpl : public EnzoRiemann "ImplFunctor's operator() method doesn't have the correct " "function signature"); - /// The following vector holds advection quantities for which passive - /// advection fluxes can be computed when they are not directly handled by a - /// Riemann Solver. Unlike normal passively advected scalars, quantities - /// listed here nominally play a more direct role in hydrodynamical - /// integration. The canonical example is "specific internal energy" - /// - /// Note: This is unrelated to Riemann Solver Fallback. - static const std::vector PassiveFallbackAdvectionQuantities; public: // interface /// 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(){ }; @@ -185,11 +175,16 @@ class EnzoRiemannImpl : public EnzoRiemann void pup (PUP::er &p); 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 str_vec_t &passive_list, - EFlt3DArray *interface_velocity) const; + EnzoEFltArrayMap &flux_map, const int dim, + const EnzoEquationOfState *eos, const int stale_depth, + const str_vec_t &passive_list, + const 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 @@ -198,88 +193,43 @@ protected : //methods EnzoEFltArrayMap &prim_map_r, EnzoEFltArrayMap &flux_map, const EFlt3DArray &density_flux, - int stale_depth, + const int stale_depth, const str_vec_t &passive_list) const throw(); 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 integrable quantities that require passive fluxes - std::vector passive_integrable_quantities_; - - /// Field categories of each quantity in passive_integrable_quantities_ - std::vector passive_integrable_categories_; + /// 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_; }; //---------------------------------------------------------------------- -// TODO: move this to EnzoRiemann so that it isn't initialized for every -// template - -template -const std::vector -EnzoRiemannImpl::PassiveFallbackAdvectionQuantities = - {"internal_energy"}; - -//---------------------------------------------------------------------- - template -EnzoRiemannImpl::EnzoRiemannImpl -(std::vector integrable_groups) +EnzoRiemannImpl::EnzoRiemannImpl(const 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()); - } +{ + integration_quantities_ = LUT::integration_quantity_names(); - // Check that all quantities in integrable_groups are in lut_groups or in - // EnzoRiemannImpl::PassiveFallbackAdvectionQuantities - 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 (std::find(PassiveFallbackAdvectionQuantities.begin(), - PassiveFallbackAdvectionQuantities.end(), - group) != PassiveFallbackAdvectionQuantities.end()){ - passive_integrable_quantities_.push_back(group); - } 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); } } - // determine the categories of each passive integrable group - for (std::string group : passive_integrable_quantities_){ - FieldCat category; - EnzoCenteredFieldRegistry::quantity_properties(group, nullptr, - &category, nullptr); - passive_integrable_categories_.push_back(category); + calculate_internal_energy_flux_ = internal_energy; + if (calculate_internal_energy_flux_){ + integration_quantities_.push_back("internal_energy"); } } @@ -290,19 +240,50 @@ void EnzoRiemannImpl::pup (PUP::er &p) { EnzoRiemann::pup(p); - p|integrable_groups_; - p|passive_integrable_quantities_; - p|passive_integrable_categories_; + 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_(const enzo_float left, + const enzo_float right, + const 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_(const enzo_float density_l, + const enzo_float pressure_l, + const enzo_float density_r, + const enzo_float pressure_r, + const enzo_float gamma, + const 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, - const EFlt3DArray &pressure_array_l, const EFlt3DArray &pressure_array_r, - EnzoEFltArrayMap &flux_map, int dim, EnzoEquationOfState *eos, - int stale_depth, const str_vec_t &passive_list, - EFlt3DArray *interface_velocity) const + EnzoEFltArrayMap &flux_map, const int dim, const EnzoEquationOfState *eos, + const int stale_depth, const str_vec_t &passive_list, + const EFlt3DArray *interface_velocity) const { const bool barotropic = eos->is_barotropic(); @@ -317,16 +298,26 @@ void EnzoRiemannImpl::solve ASSERT("EnzoRiemannImpl::solve", "currently no support for barotropic eos", !barotropic); - // 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; } + // TODO: Add special handling for barotropic equations of state + const EFlt3DArray pressure_array_l = prim_map_l.at("pressure"); + const EFlt3DArray pressure_array_r = prim_map_r.at("pressure"); + + // 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; - wl_arrays = load_array_of_fields(prim_map_l, dim); - wr_arrays = load_array_of_fields(prim_map_r, dim); - flux_arrays = load_array_of_fields(flux_map, dim); + wl_arrays = load_array_of_fields(prim_map_l, dim, true); + wr_arrays = load_array_of_fields(prim_map_r, dim, true); + flux_arrays = load_array_of_fields(flux_map, dim, false); ImplFunctor func; @@ -348,8 +339,10 @@ 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); - lutarray Ur = enzo_riemann_utils::compute_conserved(wr); + lutarray Ul = enzo_riemann_utils::compute_conserved(wl, + gamma); + lutarray Ur = enzo_riemann_utils::compute_conserved(wr, + gamma); // compute the interface fluxes lutarray Fl = enzo_riemann_utils::active_fluxes(wl, Ul, @@ -368,8 +361,11 @@ 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]); } } } @@ -381,21 +377,18 @@ void EnzoRiemannImpl::solve 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. } } //---------------------------------------------------------------------- -inline void passive_advection_helper_ -(const str_vec_t &passive_keys, - EnzoEFltArrayMap& prim_map_l, EnzoEFltArrayMap& prim_map_r, - EnzoEFltArrayMap& flux_map, const EFlt3DArray &density_flux, - int stale_depth) noexcept +template +void EnzoRiemannImpl::solve_passive_advection_ +(EnzoEFltArrayMap &prim_map_l, EnzoEFltArrayMap &prim_map_r, + EnzoEFltArrayMap &flux_map, const EFlt3DArray &density_flux, + const int stale_depth, const str_vec_t &passive_list) const throw() { - std::size_t num_keys = passive_keys.size(); + const std::size_t num_keys = passive_list.size(); if (num_keys == 0) {return;} // This was essentially transcribed from hydro_rk in Enzo: @@ -406,9 +399,9 @@ inline void passive_advection_helper_ EFlt3DArray *flux_arrays = new EFlt3DArray[num_keys]; for (std::size_t ind=0; ind0) ? wl : wr; - flux_arrays[key_ind](iz,iy,ix) = w_upwind * dens_flux; + const enzo_float dens_flux = density_flux(iz,iy,ix); + const enzo_float wl = wl_arrays[key_ind](iz,iy,ix); + const enzo_float wr = wr_arrays[key_ind](iz,iy,ix); + flux_arrays[key_ind](iz,iy,ix) = + calc_passive_scalar_flux_(wl, wr, dens_flux); } } @@ -429,35 +422,4 @@ inline void passive_advection_helper_ delete[] wl_arrays; delete[] wr_arrays; delete[] flux_arrays; } -//---------------------------------------------------------------------- - -template -void EnzoRiemannImpl::solve_passive_advection_ -(EnzoEFltArrayMap &prim_map_l, EnzoEFltArrayMap &prim_map_r, - EnzoEFltArrayMap &flux_map, const EFlt3DArray &density_flux, - int stale_depth, const str_vec_t &passive_list) const throw() -{ - // First address the passive scalars (they are listed in passive_list) - passive_advection_helper_(passive_list, prim_map_l, prim_map_r, flux_map, - density_flux, stale_depth); - - // Next address integrable quantites that have passively advected fluxes - for (std::size_t i=0; i>`. 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 2717be509a..f72eab94ac 100644 --- a/src/Enzo/enzo_EnzoRiemannUtils.hpp +++ b/src/Enzo/enzo_EnzoRiemannUtils.hpp @@ -18,20 +18,58 @@ //---------------------------------------------------------------------- +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(const enzo_float density, + const enzo_float pressure, + const enzo_float gamma) noexcept + { return pressure / ( (gamma - 1.0) * density); } + + static inline enzo_float eint_dens(const enzo_float density, const enzo_float pressure, + const enzo_float gamma) noexcept + { return pressure / (gamma - 1.0); } + + static inline enzo_float sound_speed(const enzo_float density, const enzo_float pressure, + const enzo_float gamma) noexcept + { return std::sqrt(gamma * pressure / density); } + +}; + +//---------------------------------------------------------------------- + namespace enzo_riemann_utils{ - /// Computes the conserved counterpart for every integrable primitive. - /// Integrable primitives are categorized as conserved, specific, and other + template + inline enzo_float mag_pressure(const lutarray prim) noexcept + { + enzo_float bi = (LUT::bfield_i >= 0) ? prim[LUT::bfield_i] : 0; + enzo_float bj = (LUT::bfield_j >= 0) ? prim[LUT::bfield_j] : 0; + enzo_float bk = (LUT::bfield_k >= 0) ? prim[LUT::bfield_k] : 0; + return 0.5 * (bi*bi + bj*bj + bk *bk); + } + + //---------------------------------------------------------------------- + + /// Computes the conserved counterpart for every primitive. /// - /// quantities classified as conserved are copied, quantities classified as - /// primitive are multiplied by density, and for simplicity, quantities - /// classified as other are copied (There are no obvious cases where there - /// should ever be a quanitity classified as other + /// Primitives are generally categorized as conserved or specific. A + /// "conserved" primitive is directly copied, and a "specific" primitive + /// is multiplied by the density. The main exception is pressure, which is + /// stored in prim[LUT::total_energy]. /// - /// This has been factored out EnzoRiemannImpl to reduce the number of - /// template functions created when the same LUT is reused + /// This is currently assumed to only be implemented for a non-barotropic EOS + /// with a calorically perfect EOS + /// + /// We might want to consolidate this with active_fluxes template - inline lutarray compute_conserved(const lutarray prim) noexcept + inline lutarray compute_conserved(const lutarray prim, + const enzo_float gamma) noexcept { lutarray cons; @@ -45,26 +83,25 @@ namespace enzo_riemann_utils{ for (std::size_t i = LUT::specific_start; i < LUT::NEQ; i++) { cons[i] = prim[i] * prim[LUT::density]; } - return cons; - } - //---------------------------------------------------------------------- + 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 = EOSStructIdeal::eint_dens(density, pressure, + gamma); - template - inline enzo_float mag_pressure(const lutarray prim) noexcept - { - enzo_float bi = (LUT::bfield_i >= 0) ? prim[LUT::bfield_i] : 0; - enzo_float bj = (LUT::bfield_j >= 0) ? prim[LUT::bfield_j] : 0; - enzo_float bk = (LUT::bfield_k >= 0) ? prim[LUT::bfield_k] : 0; - return 0.5 * (bi*bi + bj*bj + bk *bk); - } + const enzo_float vi = prim[LUT::velocity_i]; + const enzo_float vj = prim[LUT::velocity_j]; + const enzo_float vk = prim[LUT::velocity_k]; + const enzo_float kinetic_edens = 0.5 * density * (vi*vi + vj*vj + vk*vk); - //---------------------------------------------------------------------- + enzo_float magnetic_edens = mag_pressure(prim); - 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]); } + cons[LUT::total_energy] = internal_edens + kinetic_edens + magnetic_edens; + } + + return cons; + } //---------------------------------------------------------------------- @@ -80,8 +117,10 @@ 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 B2 = (bi*bi + bj*bj + bk *bk); + const enzo_float cs = EOSStructIdeal::sound_speed(prim_vals[LUT::density], + pressure, gamma); + const enzo_float cs2 = std::pow(cs,2); + const enzo_float B2 = (bi*bi + bj*bj + bk *bk); if (B2 == 0){ return std::sqrt(cs2); } @@ -109,8 +148,10 @@ 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 B2 = (bi*bi + bj*bj + bk *bk); + const enzo_float cs = EOSStructIdeal::sound_speed(prim_vals[LUT::density], + pressure, gamma); + const enzo_float cs2 = std::pow(cs,2); + const enzo_float B2 = (bi*bi + bj*bj + bk *bk); if (B2 == 0){ return std::sqrt(cs2); } @@ -207,7 +248,7 @@ namespace enzo_riemann_utils{ /// the mesh). This allows for appropriate loading of reconstructed fields. template inline std::array load_array_of_fields - (EnzoEFltArrayMap& map, int dim) noexcept + (EnzoEFltArrayMap& map, const int dim, const bool prim) noexcept { std::array arr; // in the case where we don't have reconstructed values (dim = -1) we assume @@ -217,9 +258,16 @@ namespace enzo_riemann_utils{ // define a lambda function to execute for every member of lut. For each // member in lut, its passed: 1. the member's name // 2. the associated index - auto fn = [coord, &arr, &map](std::string name, int index) + auto fn = [coord, prim, &arr, &map](const std::string& name, + const int index) { - if (index != -1){ arr[index] = map.at(parse_mem_name_(name, coord)); } + if (index != -1){ + if (prim && (index == LUT::total_energy)){ + arr[index] = map.at("pressure"); + } else { + arr[index] = map.at(parse_mem_name_(name, coord)); + } + } }; LUT::for_each_entry(fn); diff --git a/src/Enzo/enzo_EnzoSourceInternalEnergy.cpp b/src/Enzo/enzo_EnzoSourceInternalEnergy.cpp index d000273aa4..5cc1132ab0 100644 --- a/src/Enzo/enzo_EnzoSourceInternalEnergy.cpp +++ b/src/Enzo/enzo_EnzoSourceInternalEnergy.cpp @@ -11,9 +11,10 @@ //---------------------------------------------------------------------- void EnzoSourceInternalEnergy::calculate_source -(int dim, double dt, enzo_float cell_width, EnzoEFltArrayMap &prim_map, - EnzoEFltArrayMap &dUcons_map, EFlt3DArray &interface_velocity, - EnzoEquationOfState *eos, int stale_depth) const throw() +(const int dim, const double dt, const enzo_float cell_width, + EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &dUcons_map, + EFlt3DArray &interface_velocity, const EnzoEquationOfState *eos, + const int stale_depth) const throw() { // SANITY CHECKS: ASSERT("EnzoSourceInternalEnergy::calculate_source", @@ -24,7 +25,7 @@ void EnzoSourceInternalEnergy::calculate_source "The EOS can't be barotropic and use the dual energy formalism.", !(eos->is_barotropic()) ); - enzo_float dtdx = dt/cell_width; + const enzo_float dtdx = dt/cell_width; EnzoPermutedCoordinates coord(dim); @@ -42,14 +43,11 @@ void EnzoSourceInternalEnergy::calculate_source CSlice full_ax(nullptr, nullptr); // load cell-centered quantities. - // define: rho_center(k,j,i) -> rho(k,j,i+1) - // eint_center(k,j,i) -> eint(k,j,i+1) + // define: pressure_center(k,j,i) -> pressure(k,j,i+1) // deint_dens_center(k,j,i) -> deint_dens(k,j,i+1) - EFlt3DArray rho, eint, rho_center, eint_center; - rho = prim_map.get("density", stale_depth); - eint = prim_map.get("internal_energy", stale_depth); - rho_center = coord.get_subarray(rho, full_ax, full_ax, CSlice(1, -1)); - eint_center = coord.get_subarray(eint, full_ax, full_ax, CSlice(1, -1)); + EFlt3DArray pressure = prim_map.get("pressure", stale_depth); + EFlt3DArray pressure_center = coord.get_subarray(pressure, full_ax, full_ax, + CSlice(1, -1)); EFlt3DArray deint_dens = dUcons_map.get("internal_energy", stale_depth); EFlt3DArray deint_dens_center = coord.get_subarray(deint_dens, full_ax, @@ -71,12 +69,13 @@ void EnzoSourceInternalEnergy::calculate_source enzo_float gm1 = eos->get_gamma() - 1.; - for (int iz=0; iz