Skip to content

Commit

Permalink
This commit is from Will Hicks. It undoes changes introduced in PR en…
Browse files Browse the repository at this point in the history 
…zo-project#173
  • Loading branch information
@WillHicks96 @mabruzzo
WillHicks96 authored and mabruzzo committed Jan 26, 2023
1 parent ebb3ab0 commit c01ebb8
Show file tree
Hide file tree
Showing 12 changed files with 348 additions and 498 deletions.
4 changes: 2 additions & 2 deletions src/Enzo/enzo_EnzoComputeCoolingTime.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,8 +73,8 @@ void EnzoComputeCoolingTime::compute_(Block * block,
"Grackle must be enabled in order to compute the cooling time",
enzo::config()->method_grackle_use_grackle );
const EnzoMethodGrackle* grackle_method = enzo::grackle_method();
grackle_method->calculate_cooling_time(EnzoFieldAdaptor(block, i_hist_), ct,
0, grackle_units, grackle_fields);
grackle_method->calculate_cooling_time(block, ct,
grackle_units, grackle_fields, i_hist_);
}

#endif
220 changes: 64 additions & 156 deletions src/Enzo/enzo_EnzoComputePressure.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,209 +43,117 @@ void EnzoComputePressure::compute ( Block * block) throw()

if (!block->is_leaf()) return;

Field field = block->data()->field();
EnzoBlock * enzo_block = enzo::block(block);
Field field = enzo_block->data()->field();

ASSERT("EnzoComputePressure::compute()",
"'pressure' must be defined as a permanent field",
field.field_id("pressure") >= 0);
// TODO: possibly check that pressure is cell-centered
enzo_float * p = field.is_field("pressure") ?
(enzo_float*) field.values("pressure", i_hist_) : NULL;

compute(block, (enzo_float*)field.values("pressure", i_hist_));
if (!p){
ERROR("EnzoComputePressure::compute()",
" 'pressure' is not defined as a permanent field");
}

compute(block, p);
}

//----------------------------------------------------------------------

void EnzoComputePressure::compute (Block * block, enzo_float *p,
int stale_depth /* 0 */) throw()
void EnzoComputePressure::compute (Block * block, enzo_float * p) throw()
{

if (!block->is_leaf()) return;

compute_(block, p, stale_depth);
compute_(block, p);
}

//----------------------------------------------------------------------

void EnzoComputePressure::compute_(Block * block,
enzo_float * p,
int stale_depth /* 0 */
enzo_float * p
#ifdef CONFIG_USE_GRACKLE
, code_units * grackle_units, /*NULL*/
grackle_field_data * grackle_fields /*NULL*/
#endif
)
{
// make a CelloArray that wraps p
Field field = block->data()->field();
int gx,gy,gz;
field.ghost_depth (field.field_id("density"),&gx,&gy,&gz);
int nx,ny,nz;
field.size (&nx,&ny,&nz);
CelloArray<enzo_float,3> p_arr(p, nz + 2*gz, ny + 2*gy, nx + 2*gx);

EnzoFieldAdaptor f_adaptor(block, i_hist_);

bool dual_energy = !enzo::fluid_props()->dual_energy_config().is_disabled();
bool mhd = field.is_field("bfield_x");

EnzoComputePressure::compute_pressure(f_adaptor, p_arr, mhd, dual_energy,
gamma_, stale_depth);
}

//----------------------------------------------------------------------

namespace { // local function

template<class K>
void exec_loop_(int mz, int my, int mx, int stale_depth, K& kernel){

const int rank = cello::rank();

const int ix_start = stale_depth;
const int ix_stop = mx - stale_depth;

const int iy_start = (rank > 1) ? stale_depth : 0;
const int iy_stop = (rank > 1) ? my - stale_depth : my;

const int iz_start = (rank > 2) ? stale_depth : 0;
const int iz_stop = (rank > 2) ? mz - stale_depth : mz;

if ((ix_start >= ix_stop) | (iy_start >= iy_stop) | (iz_start >= iz_stop)){
ERROR("exec_loop_", "stale_depth is too large");
} else if (stale_depth < 0){
ERROR("exec_loop_", "stale_depth is negative");
}

for (int iz = iz_start; iz < iz_stop; iz++){
for (int iy = iy_start; iy < iy_stop; iy++){
for (int ix = ix_start; ix < ix_stop; ix++){
kernel(iz, iy, ix);
}
}
}
}
if (!block->is_leaf()) return;

}
EnzoBlock * enzo_block = enzo::block(block);

//----------------------------------------------------------------------
Field field = enzo_block->data()->field();

void EnzoComputePressure::compute_pressure
(const EnzoFieldAdaptor& fadaptor,
const CelloArray<enzo_float, 3>& p,
bool mhd,
bool dual_energy,
double gamma,
int stale_depth, /* 0 */
bool ignore_grackle /*false*/
#ifdef CONFIG_USE_GRACKLE
, code_units * grackle_units, /*nullptr*/
grackle_field_data * grackle_fields /*nullptr*/
#endif
) throw()
{
bool mhd = field.is_field("bfield_x");

if (enzo::config()->method_grackle_use_grackle & !ignore_grackle){
if (enzo::config()->method_grackle_use_grackle){
#ifdef CONFIG_USE_GRACKLE
// the following assertion is not strictly necessary (the problem will be
// caught later in EnzoMethodGrackle), but this is more informative...
ASSERT("EnzoMethodGrackle::calculate_pressure",
"until PR #106 is merged into grackle, stale_depth must be 0 since "
"grackle's local_calculate_pressure ignores the grid_start and "
"grid_stop members of grackle_field_data",
stale_depth == 0);

if (!fadaptor.consistent_with_field_strides(p)){
ERROR("EnzoMethodGrackle::calculate_pressure",
"When using grackle to compute pressure, the output array must "
"have identical strides to the fields.");
}
const EnzoMethodGrackle* grackle_method = enzo::grackle_method();
grackle_method->calculate_pressure(fadaptor, p.data(), stale_depth,
grackle_units, grackle_fields);
grackle_method->calculate_pressure(block, p, grackle_units,
grackle_fields, i_hist_);
#else
ERROR("EnzoComputePressure::compute_()",
"Attempting to compute pressure with method Grackle "
"but Enzo-E has not been compiled with Grackle (set use_grackle = 1) \n");
#endif
} else {

using RdOnlyEFltArr = CelloArray<const enzo_float, 3>;

const RdOnlyEFltArr d = fadaptor.view("density");
} else {

enzo_float gm1 = gamma - 1.0;
const int rank = cello::rank();

const int mz = d.shape(0);
const int my = d.shape(1);
const int mx = d.shape(2);
const enzo_float * d = (enzo_float*) field.values("density", i_hist_);

ASSERT("EnzoComputePressure::compute_pressure",
"Shapes must be the same",
((mz == p.shape(0)) & (my == p.shape(1)) & (mx == p.shape(2))));
const enzo_float * vx = (enzo_float*) field.values("velocity_x", i_hist_);
const enzo_float * vy = (enzo_float*) field.values("velocity_y", i_hist_);
const enzo_float * vz = (enzo_float*) field.values("velocity_z", i_hist_);

if (dual_energy) {
const enzo_float * bx = mhd ?
(enzo_float*) field.values("bfield_x", i_hist_) : nullptr;
const enzo_float * by = mhd ?
(enzo_float*) field.values("bfield_y", i_hist_) : nullptr;
const enzo_float * bz = mhd ?
(enzo_float*) field.values("bfield_z", i_hist_) : nullptr;

const RdOnlyEFltArr ie = fadaptor.view("internal_energy");
const enzo_float * te = (enzo_float*) field.values("total_energy", i_hist_);
const enzo_float * ie = (enzo_float*) field.values("internal_energy", i_hist_);

auto loop_body = [=](int iz, int iy, int ix)
{ p(iz,iy,ix) = gm1 * d(iz,iy,ix) * ie(iz,iy,ix); };
exec_loop_(mz, my, mx, stale_depth, loop_body);
int nx,ny,nz;
field.size(&nx,&ny,&nz);

} else { // not using dual energy formalism
int gx,gy,gz;
field.ghost_depth (0,&gx,&gy,&gz);
if (rank < 2) gy = 0;
if (rank < 3) gz = 0;

const int rank = cello::rank();
const RdOnlyEFltArr te = fadaptor.view("total_energy");
int m = (nx+2*gx) * (ny+2*gy) * (nz+2*gz);
enzo_float gm1 = gamma_ - 1.0;

// fetch velocity arrays
const RdOnlyEFltArr vx = fadaptor.view("velocity_x");
const RdOnlyEFltArr vy = (rank >= 2)
? fadaptor.view("velocity_y") : RdOnlyEFltArr();
const RdOnlyEFltArr vz = (rank >= 3)
? fadaptor.view("velocity_z") : RdOnlyEFltArr();
if (true) { //(enzo::config()->ppm_dual_energy) {

// fetch bfield arrays
const RdOnlyEFltArr bx = (mhd)
? fadaptor.view("bfield_x") : RdOnlyEFltArr();
const RdOnlyEFltArr by = (mhd & (rank >= 2))
? fadaptor.view("bfield_y") : RdOnlyEFltArr();
const RdOnlyEFltArr bz = (mhd & (rank >= 3))
? fadaptor.view("bfield_z") : RdOnlyEFltArr();
for (int i=0; i<m; i++) {
p[i] = gm1 * d[i] * ie[i];
}

} else {
if (rank == 1) {

auto loop_body = [=](int iz, int iy, int ix)
{
enzo_float ke = 0.5 * (vx(iz,iy,ix) * vx(iz,iy,ix));
enzo_float me_den = mhd ? 0.5*(bx(iz,iy,ix) * bx(iz,iy,ix)) : 0.;
p(iz,iy,ix) = gm1 * (d(iz,iy,ix) * (te(iz,iy,ix) - ke) - me_den);
};
exec_loop_(mz, my, mx, stale_depth, loop_body);

for (int i=0; i<m; i++) {
enzo_float ke = 0.5*vx[i]*vx[i];
enzo_float me_den = mhd ? 0.5*bx[i]*bx[i] : 0.;
p[i] = gm1 * (d[i] * (te[i] - ke) - me_den);
}
} else if (rank == 2) {

auto loop_body = [=](int iz, int iy, int ix)
{
enzo_float ke = 0.5*(vx(iz,iy,ix) * vx(iz,iy,ix) +
vy(iz,iy,ix) * vy(iz,iy,ix));
enzo_float me_den = mhd ? 0.5*(bx(iz,iy,ix) * bx(iz,iy,ix) +
by(iz,iy,ix) * by(iz,iy,ix)) : 0.;
p(iz,iy,ix) = gm1 * (d(iz,iy,ix) * (te(iz,iy,ix) - ke) - me_den);
};
exec_loop_(mz, my, mx, stale_depth, loop_body);

for (int i=0; i<m; i++) {
enzo_float ke = 0.5*(vx[i]*vx[i] + vy[i]*vy[i]);
enzo_float me_den = mhd ? 0.5*(bx[i]*bx[i] + by[i]*by[i]) : 0.;
p[i] = gm1 * (d[i] * (te[i] - ke) - me_den);
}
} else if (rank == 3) {

auto loop_body = [=](int iz, int iy, int ix)
{
enzo_float ke = 0.5*(vx(iz,iy,ix) * vx(iz,iy,ix) +
vy(iz,iy,ix) * vy(iz,iy,ix) +
vz(iz,iy,ix) * vz(iz,iy,ix));
enzo_float me_den = mhd ? 0.5*(bx(iz,iy,ix) * bx(iz,iy,ix) +
by(iz,iy,ix) * by(iz,iy,ix) +
bz(iz,iy,ix) * bz(iz,iy,ix)) : 0.;
p(iz,iy,ix) = gm1 * (d(iz,iy,ix) * (te(iz,iy,ix) - ke) - me_den);
};
exec_loop_(mz, my, mx, stale_depth, loop_body);

for (int i=0; i<m; i++) {
enzo_float ke = 0.5*(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
enzo_float me_den = mhd ?
0.5*(bx[i]*bx[i] + by[i]*by[i] + bz[i]*bz[i]) : 0.;
p[i] = gm1 * (d[i] * (te[i] - ke) - me_den);
}
}
}
}
Expand Down
70 changes: 9 additions & 61 deletions src/Enzo/enzo_EnzoComputePressure.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,72 +39,20 @@ class EnzoComputePressure : public Compute {
return "pressure";
}

/// Perform the computation on the block and store the results in the
/// "pressure" field
/// Perform the computation on the block
void compute( Block * block) throw();

/// Perform the computation on the block and store the result in the provided
/// array.
///
/// @param[in] block Where input arrays are loaded from
/// @param[out] p Array that should hold the result.
/// @param[in] stale_depth The number of the outermost zones to exclude from
/// the calculation
///
/// This function will rely on the Grackle-supplied function for computing
/// the pressure if the simulation is configured to use `EnzoMethodGrackle`
/// (in this case the `gamma` value passed to the constructor is ignored).
/// Otherwise, it falls back to an alternate version of the calculation.
void compute( Block * block,
enzo_float* p,
int stale_depth = 0) throw();

/// Perform the calculation with the provided `grackle_units` or
/// `grackle_fields` object.
///
/// @note
/// This is primarily used by `EnzoInitialGrackleTest`
void compute( Block * block, enzo_float * p) throw();

void compute_(Block * block,
enzo_float * p,
int stale_depth = 0
#ifdef CONFIG_USE_GRACKLE
, code_units * grackle_units = nullptr,
grackle_field_data * grackle_fields = nullptr
#endif
);

/// static method to compute thermal pressure
///
/// @param[in] fadaptor Contains arrays of quantities used to compute the
/// thermal pressure.
/// @param[out] p Array where the thermal pressure is to be stored
/// @param[in] mhd Whether the simulation includes magnetic fields
/// @param[in] dual_energy Whether the simulation uses the dual-energy
/// formalism
/// @param[in] gamma Specifies the adiabatic index. This is entirely ignored
/// if the Grackle-supplied function for computing pressure is used.
/// @param[in] stale_depth The number of the outermost zones to exclude from
/// the calculation.
/// @param[in] ignore_grackle indicates whether to ignore the Grackle
/// routine for computing pressure. This exists to force the usage of the
/// `gamma` parameter (and can be used to ignore the effects of molecular
/// hydrogen on the adiabatic index)
///
/// By default, this function will rely on the Grackle-supplied function for
/// computing the pressure if the simulation is configured to use
/// `EnzoMethodGrackle`. The `ignore_grackle` parameter only has meaning in
/// this scenario (if the simulation does not use `EnzoMethodGrackle`, then
/// the value of `ignore_grackle` is meaningless).
static void compute_pressure(const EnzoFieldAdaptor& fadaptor,
const CelloArray<enzo_float, 3>& p,
bool mhd, bool dual_energy, double gamma,
int stale_depth = 0,
bool ignore_grackle = false
enzo_float * p
#ifdef CONFIG_USE_GRACKLE
, code_units * grackle_units = nullptr,
grackle_field_data * grackle_fields = nullptr
, code_units * grackle_units = NULL,
grackle_field_data * grackle_fields = NULL
#endif
) throw();
);



protected: // attributes

Expand Down
6 changes: 3 additions & 3 deletions src/Enzo/enzo_EnzoComputeTemperature.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,14 +89,14 @@ void EnzoComputeTemperature::compute_(Block * block,

Field field = enzo_block->data()->field();

const enzo_float gamma = enzo::fluid_props()->gamma();
const enzo_float gamma = enzo::fluid_props()->gamma(); //EnzoBlock::Gamma[cello::index_static()];

if (enzo::config()->method_grackle_use_grackle){

#ifdef CONFIG_USE_GRACKLE
const EnzoMethodGrackle* grackle_method = enzo::grackle_method();
grackle_method->calculate_temperature(EnzoFieldAdaptor(block, i_hist_), t,
0, grackle_units, grackle_fields);
grackle_method->calculate_temperature(block, t, grackle_units,
grackle_fields, i_hist_);
#else
ERROR("EnzoComputeTemperature::compute_()",
"Attempting to compute temperature with method Grackle "
Expand Down
Loading

0 comments on commit c01ebb8

Please sign in to comment.