Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

VL+CT QoL Improvements #94

Merged
merged 15 commits into from
Nov 30, 2021
Merged

VL+CT QoL Improvements #94

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
00214d0
First step in shifting the Riemann Solver to compute the fluxes from …
mabruzzo Apr 22, 2021
41070ad
Refactored the internal energy source term to make use of pressure (i…
mabruzzo Apr 25, 2021
4b23eca
Started updating EnzoMethodMHDVlct and EnzoEquationOfState to replace…
mabruzzo Apr 25, 2021
63ad9b7
Refactored EnzoMethodMHDVlct so that the integration map holds the pa…
mabruzzo Apr 30, 2021
a8fa814
Renamed EnzoIntegrableUpdate as EnzoIntegrationQuanUpdate and changed…
mabruzzo May 1, 2021
9f3ab8c
Tweaked EnzoBfieldMethod to reflect the change in the variable names
mabruzzo May 8, 2021
270d9e9
Updated the website documentation for the hydro/mhd toolkit (except t…
mabruzzo May 4, 2021
de32f09
Updates to EnzoRiemann related to the shift to using primitives to co…
mabruzzo May 5, 2021
9ffe97d
Merge branch 'master' into vlct-varnames
mabruzzo May 20, 2021
ee5ee15
Tweaked the tolerance on the VL+CT dual-energy cloud test for the HLL…
mabruzzo May 27, 2021
e0c43df
Apply most suggestions from code review
mabruzzo May 28, 2021
cba3d9d
Fixed compilation bug introduced in prior commit.
mabruzzo May 28, 2021
ba45e50
Fixed a minor bug in EnzoLazyPassiveScalarFieldList
mabruzzo May 28, 2021
860d2a6
tweaked description of EnzoCenteredFieldRegistry in website docs (pos…
mabruzzo Nov 13, 2021
722f041
Fixing typos in documentation (to address @gregbryan's PR Review comm…
mabruzzo Nov 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
593 changes: 282 additions & 311 deletions doc/source/devel/hydro_infrastructure.rst
View file
Open in desktop

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions input/vlct/run_dual_energy_cloud_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
2 changes: 1 addition & 1 deletion src/Enzo/_enzo.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
2 changes: 1 addition & 1 deletion src/Enzo/enzo.CI
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ module enzo {
PUPable EnzoMethodMHDVlct;

PUPable EnzoEOSIdeal;
PUPable EnzoIntegrableUpdate;
PUPable EnzoIntegrationQuanUpdate;
PUPable EnzoReconstructorNN;
// EnzoReconstructorPLMEnzoRKLim
PUPable EnzoReconstructorPLM<PLM_EnzoRKLimiter>;
Expand Down
10 changes: 5 additions & 5 deletions src/Enzo/enzo_EnzoBfieldMethod.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
Expand Down
36 changes: 18 additions & 18 deletions src/Enzo/enzo_EnzoBfieldMethodCT.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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++){
Expand All @@ -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"};
Expand All @@ -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]);
Comment on lines +264 to +267
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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]);
const EFlt3DArray &velocity_j = integration_map.at(v_names[j]);
const EFlt3DArray &velocity_k = integration_map.at(v_names[k]);
const EFlt3DArray &bfield_j = integration_map.at(b_names[j]);
const EFlt3DArray &bfield_k = integration_map.at(b_names[k]);

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've thought about this exact sort of change a lot, and I'm not sure whether this is an improvement or not.

The thing is that EFlt3DArray can be thought of as a shared pointer. By making a copy of it (what we have now), we pay the overhead of incrementing/decrementing the reference count (although the current case probably uses move semantics). By storing a reference to it, we avoid increments/decrements, but now we effectively have a reference to a pointer (which is similar to a pointer to a pointer). I'm not sure what's better

An alternative idea is just to mark velocity_j, velocity_k, bfield_j, and bfield_k as const.

The same sort of semantics apply to Kokkos::Views, so I'm willing to go with your instinct here. But before I change it, do you have anything else to add?.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The thing is that EFlt3DArray can be thought of as a shared pointer. By making a copy of it (what we have now), we pay the overhead of incrementing/decrementing the reference count (although the current case probably uses move semantics). By storing a reference to it, we avoid increments/decrements, but now we effectively have a reference to a pointer (which is similar to a pointer to a pointer). I'm not sure what's better

We've been using essentially const Kokkos:View& recently, I appreciate that we keep the same object. Sometimes it can make debugging easier since the address of the container doesn't change between contexts.

The other confusing piece to me is that all of our arrays can probably be made const &, since we even though we're changing data in the array we're not changing the data itself.

The same sort of semantics apply to Kokkos::Views, so I'm willing to go with your instinct here. But before I change it, do you have anything else to add?.

I'm not planning on adding any more suggestions. I would suggest to try batching all of my suggestions together into one commit via github, then see if it compiles and runs. I've probably missed a few things that would also need to be constant but that should be straightfoward to find via the compiler. Worst case you could undo the commit entirely.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other confusing piece to me is that all of our arrays can probably be made const &, since we even though we're changing data in the array we're not changing the data itself

That's part of the reason why I've largely avoided referring to EFlt3DArray as a const in most places yet.

When I originally developed CelloArray, a lot of it's design was strongly motivated by my familiarity with NumPy. I wasn't fond of the way that traditional "container-semantics" will produce a deepcopy when you perform an assignment (it also makes subarray creation more confusing). It was originally my intention to prevent modification of the contents for const CelloArray, but that led to some strange quirks in the Copy-Constructor and Copy-Assignment operation. (Some stuff with temporaries break. It's confusing to new users why certain things break, but I don't think it actually limits what we want to do in the main codebase)

At around the time that I was struggling with this, I started reading about Kokkos and I realized that CelloArray has pointer-semantics just like Kokkos::View (a better name for CelloArray is probably CelloView), so I decided to lean into it.

I'm working on something that might help the situation. I think I mentioned this before, but I've slowly been working on a new PR, in which I introduce support for CelloArray<const T, D> which prevents the contents of the array from being mutated (I needed to make sure that this case element access works right support and add code to facilitate casts from CelloArray<T, D> and CelloArray<const T, D>). I got a little side-tracked over the past 2 weeks, but I hope to get back to it soon-ish. Once I'm done, it's my intention to alias CelloArray<const enzo_float, D> as something like RdOnlyEFlt3DArray (I'm definitely open to better names) and start introducing that throughout the VL+CT integrator.

But, It might be worthwhile to discuss changing the CelloArray semantics:

  • If you and @pgrete, don't think it would be confusing, it might be worth revisiting the idea of making const CelloArray totally immutable. I'm not convinced that it's completely impossible. It might be easier now that I have a working test suite and I'm planning to remove the bulk assignment stuff from the CelloArray API (which will drastically simplify the implementation).
  • I'd also be willing to consider a shift to container semantics (I don't think this would break anything). But if we did that, I would want some help with improving

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just as an aside, the aforementioned support for CelloArray<const T, D> is being implemented as part of PR #113

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the EFlt3DArray arrays that are currently used, I do think using const EFlt3DArray& is preferable to EFlt3DArray since the array shouldn't be changed even if the data is changed, and using the compiler to enforce that is preferable.

I'm working on something that might help the situation. I think I mentioned this before, but I've slowly been working on a new PR, in which I introduce support for CelloArray<const T, D> which prevents the contents of the array from being mutated (I needed to make sure that this case element access works right support and add code to facilitate casts from CelloArray<T, D> and CelloArray<const T, D>).

This would be the preferable syntax so that it's clear the array points to const data.

I got a little side-tracked over the past 2 weeks, but I hope to get back to it soon-ish. Once I'm done, it's my intention to alias CelloArray<const enzo_float, D> as something like RdOnlyEFlt3DArray (I'm definitely open to better names) and start introducing that throughout the VL+CT integrator.

I'd lean against aliasing the type. The code would be more explicit using CelloArray<const enzo_float, D> and CelloArray<enzo_float, D> everywhere so that it's clear that it's a CelloArray and references const or non-const data.

* If you and @pgrete, don't think it would be confusing, it might be worth revisiting the idea of making `const CelloArray` totally immutable. I'm not convinced that it's completely impossible. It might be easier now that I have a working test suite and I'm planning to remove the bulk assignment stuff from the `CelloArray` API (which will drastically simplify the implementation).

Make it so const CelloArray means the data it points to can't be modified? I'm against that approach, I think it makes sense to preserve CelloArray<const T, D> vs. const CelloArray<T, D> type.


for (int iz=stale_depth; iz<efield.shape(0)-stale_depth; iz++) {
for (int iy=stale_depth; iy<efield.shape(1)-stale_depth; iy++) {
Expand Down Expand Up @@ -530,15 +530,15 @@ void EnzoBfieldMethodCT::compute_edge_efield
//----------------------------------------------------------------------

void EnzoBfieldMethodCT::compute_all_edge_efields
(EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &xflux_map,
(EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &xflux_map,
EnzoEFltArrayMap &yflux_map, EnzoEFltArrayMap &zflux_map,
EFlt3DArray &center_efield, std::array<EFlt3DArray,3> &edge_efield_l,
std::array<EFlt3DArray,3> &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;
Expand Down
24 changes: 12 additions & 12 deletions src/Enzo/enzo_EnzoBfieldMethodCT.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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,
Expand Down Expand Up @@ -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 &center_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
Expand Down Expand Up @@ -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.
Expand All @@ -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 &center_efield, std::array<EFlt3DArray,3> &edge_efield_l,
std::array<EFlt3DArray,3> &weight_l, int stale_depth);
Expand Down
6 changes: 5 additions & 1 deletion src/Enzo/enzo_EnzoCenteredFieldRegistry.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,14 @@ const ft_map EnzoCenteredFieldRegistry::field_table_ = {
//----------------------------------------------------------------------

std::vector<std::string> EnzoCenteredFieldRegistry::get_registered_quantities
(bool enumerate_components)
(const bool enumerate_components, const bool exclude_non_active_advection)
{
std::vector<std::string> 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"));
Expand Down
4 changes: 3 additions & 1 deletion src/Enzo/enzo_EnzoCenteredFieldRegistry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::string> 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)
///
Expand Down
2 changes: 1 addition & 1 deletion src/Enzo/enzo_EnzoEFltArrayMap.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
Loading