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

Light refactor of EnzoMethodPmDeposit::compute and gas-deposition drift-time bugfix #239

Merged
merged 4 commits into from
Jul 1, 2022
Merged

Light refactor of EnzoMethodPmDeposit::compute and gas-deposition drift-time bugfix #239

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fea421f
Light refactoring of EnzoMethodPmDeposit.
mabruzzo Jun 25, 2022
193200e
Additional refactoring of EnzoMethodPmDeposit.cpp.
mabruzzo Jun 27, 2022
eb7698d
Fixed bug relating to drift-time for gas-deposition in EnzoMethodPmDe…
mabruzzo Jun 27, 2022
2b63d7d
clarifying comment.
mabruzzo Jun 27, 2022
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
8 changes: 0 additions & 8 deletions input/Gravity/jeans_wave/initial_jeans.incl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,6 @@

list = ["pm_deposit", "gravity"];

pm_deposit {
# every cycle, this method tries to "drift" the density field forward
# (like it's a grid of particles) by alpha*timestep.

alpha = 0.0; # this test-problem encounters issues when this is
# non-zero
}

gravity {
solver = "cg";
# choose mass units to adjust the value of the jeans wavelength is 0.5.
Expand Down
303 changes: 188 additions & 115 deletions src/Enzo/enzo_EnzoMethodPmDeposit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,71 +95,46 @@ void EnzoMethodPmDeposit::pup (PUP::er &p)

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

void EnzoMethodPmDeposit::compute ( Block * block) throw()
{

if (enzo::simulation()->cycle() == enzo::config()->initial_cycle) {
// Check if the gravity method is being used and that pm_deposit
// precedes the gravity method.
ASSERT("EnzoMethodPmDeposit",
"Error: pm_deposit method must precede gravity method.",
enzo::problem()->method_precedes("pm_deposit", "gravity"));
}

if (block->is_leaf()) {

namespace { // define local helper functions in anonymous namespace

/// deposits mass from all gravitating particles onto density_particle_arr
///
/// @param[out] density_particle_arr The array where the deposited mass
/// density is stored
/// @param[in] Block Contains the particle data to use for accumulation
/// @param[in] dt_div_cosmoa Length of time to "drift" the particles before
/// before deposition divided by the scale factor (computed for the time
/// after particles have been drifted)
/// @param[in] inv_vol One divided by the volume of a cell. When
/// `cello::rank()` is 2, this is really "inverse-area" and when it's 1,
/// this is really "inverse-cell-width". In cosmological simulations this
/// should be the comoving quantity.
/// @param[in] mx,my,mz Specifies the number of cells along each
/// dimension of an array (including ghost cells)
/// @param[in] gx,gy,gz Specifies the number of cells in the ghost zone
/// for each dimensions
void deposit_particles_(const CelloArray<enzo_float,3>& density_particle_arr,
Block* block, double dt_div_cosmoa, double inv_vol,
int mx, int my, int mz,
int gx, int gy, int gz)
{
Particle particle (block->data()->particle());
Field field (block->data()->field());

int rank = cello::rank();
enzo_float * de_t = (enzo_float *)
field.values("density_total");
enzo_float * de_p = (enzo_float *)
field.values("density_particle");
enzo_float * de_pa = (enzo_float *)
field.values("density_particle_accumulate");

int mx,my,mz;
field.dimensions(0,&mx,&my,&mz);
int nx,ny,nz;
field.size(&nx,&ny,&nz);
int gx,gy,gz;
field.ghost_depth(0,&gx,&gy,&gz);

// NOTE: density_total is now cleared in EnzoMethodGravity to
// instead of here to possible race conditions with refresh. This
// means EnzoMethodPmDeposit ("pm_deposit") currently CANNOT be
// used without EnzoMethodGravity ("gravity")
// compute extent of the active zone
const int nx = mx - 2 * gx;
const int ny = (rank >=2) ? my - 2 * gy : 1;
const int nz = (rank >=3) ? mz - 2 * gz : 1;

const int m = mx*my*mz;
std::fill_n(de_p,m,0.0);
std::fill_n(de_pa,m,0.0);
enzo_float * de_p = density_particle_arr.data();

// Get block extents and cell widths
// Get block extents
double xm,ym,zm;
double xp,yp,zp;
double hx,hy,hz;
block->lower(&xm,&ym,&zm);
block->upper(&xp,&yp,&zp);
block->cell_width(&hx,&hy,&hz);

// To calculate densities from particles with "mass" attributes
// or constants, we need the inverse volume of cells in this block.
double inv_vol = 1.0;
if (rank >= 1) inv_vol /= hx;
if (rank >= 2) inv_vol /= hy;
if (rank >= 3) inv_vol /= hz;

// Get cosmological scale factors, if cosmology is turned on
enzo_float cosmo_a=1.0;
enzo_float cosmo_dadt=0.0;
EnzoPhysicsCosmology * cosmology = enzo::cosmology();
if (cosmology) {
cosmology->compute_expansion_factor(&cosmo_a,&cosmo_dadt,
block->time() + alpha_*block->dt());
}

const double dt = alpha_ * block->dt() / cosmo_a;

// Get the number of particle types in the "is_gravitating" group
ParticleDescr * particle_descr = cello::particle_descr();
Expand Down Expand Up @@ -239,7 +214,7 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
#endif

for (int ip=0; ip<np; ip++) {
double x = xa[ip*dp] + vxa[ip*dv]*dt;
double x = xa[ip*dp] + vxa[ip*dv]*dt_div_cosmoa;

double tx = nx*(x - xm) / (xp - xm) - 0.5;

Expand Down Expand Up @@ -284,8 +259,8 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()

for (int ip=0; ip<np; ip++) {

double x = xa[ip*dp] + vxa[ip*dv]*dt;
double y = ya[ip*dp] + vya[ip*dv]*dt;
double x = xa[ip*dp] + vxa[ip*dv]*dt_div_cosmoa;
double y = ya[ip*dp] + vya[ip*dv]*dt_div_cosmoa;

double tx = nx*(x - xm) / (xp - xm) - 0.5;
double ty = ny*(y - ym) / (yp - ym) - 0.5;
Expand Down Expand Up @@ -358,9 +333,9 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()

// Copy batch particle velocities to temporary block field velocities

double x = xa[ip*dp] + vxa[ip*dv]*dt;
double y = ya[ip*dp] + vya[ip*dv]*dt;
double z = za[ip*dp] + vza[ip*dv]*dt;
double x = xa[ip*dp] + vxa[ip*dv]*dt_div_cosmoa;
double y = ya[ip*dp] + vya[ip*dv]*dt_div_cosmoa;
double z = za[ip*dp] + vza[ip*dv]*dt_div_cosmoa;

double tx = nx*(x - xm) / (xp - xm) - 0.5;
double ty = ny*(y - ym) / (yp - ym) - 0.5;
Expand Down Expand Up @@ -450,18 +425,58 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()

} // Loop over particle types in "is_gravitating" group

//--------------------------------------------------
// Add gas density
//--------------------------------------------------
}

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

/// deposits mas density from gas onto density_tot_arr
///
/// @param[in, out] density_tot_arr The array where density gets accumulated
/// @param[in] field Contains the field data to use for accumulation
/// @param[in] dt Length of time to "drift" the density field before
/// deposition divided by the scale factor at the deposition time
/// @param[in] hx_prop,hy_prop,hz_prop The width of cell along
/// each axis. These specify the proper lengths at the time that we
/// deposit the density (after any drift).
/// @param[in] mx,my,mz Specifies the number of cells along each
/// dimension of an array (including ghost cells)
/// @param[in] gx,gy,gz Specifies the number of cells in the ghost zone
/// for each dimensions
void deposit_gas_(const CelloArray<enzo_float, 3>& density_tot_arr,
Field& field, enzo_float dt_div_cosmoa,
enzo_float hx_prop, enzo_float hy_prop, enzo_float hz_prop,
int mx, int my, int mz,
int gx, int gy, int gz){

int rank = cello::rank();
const int m = mx*my*mz;

// compute extent of the active zone
int nx = mx - 2 * gx;
int ny = (rank >=2) ? my - 2 * gy : 1;
int nz = (rank >=3) ? mz - 2 * gz : 1;

// retrieve primary fields needed for depositing gas density
enzo_float * de = (enzo_float *) field.values("density");
enzo_float * vxf = (enzo_float *) field.values("velocity_x");
enzo_float * vyf = (enzo_float *) field.values("velocity_y");
enzo_float * vzf = (enzo_float *) field.values("velocity_z");

// allocate and zero-initialize scratch arrays for missing velocity
// components.
std::vector<enzo_float> vel_scratch(m*(3 - rank), 0.0);
if (rank < 2) vyf = vel_scratch.data();
if (rank < 3) vzf = vel_scratch.data() + m;

enzo_float * temp = new enzo_float [4*m];
enzo_float * de_gas = new enzo_float [m];
enzo_float * rfield = new enzo_float [m];
// deposited_gas_density is a temporary array that just includes cells in
// the active zone
const CelloArray<enzo_float, 3> deposited_gas_density(nz,ny,nx);
// CelloArray sets elements to zero by default (making next line redundant)
std::fill_n(deposited_gas_density.data(), nx*ny*nz, 0.0);

std::fill_n(temp, 4*m, 0.0);
std::fill_n(de_gas, m, 0.0);
std::fill_n(rfield, m, 0.0);
// allocate temporary arrays
std::vector<enzo_float> temp(4*m, 0.0);
std::vector<enzo_float> rfield(m, 0.0);

int gxi=gx;
int gyi=gy;
Expand All @@ -472,66 +487,124 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
int i0 = 0;
int i1 = 1;

/* Stefan: Not sure if these next four lines are correct.
I include the factors of cosmo_a for consistency with
previous version. Also, not sure why dtf is
initialised with the value of alpha_.
*/
FORTRAN_NAME(dep_grid_cic)(de, deposited_gas_density.data(), temp.data(),
vxf, vyf, vzf,
&dt_div_cosmoa, rfield.data(), &rank,
&hx_prop,&hy_prop,&hz_prop,
&mx,&my,&mz,
&gxi,&gyi,&gzi,
&nxi,&nyi,&nzi,
&i0,&i0,&i0,
&nx,&ny,&nz,
&i1,&i1,&i1);

// build a slice of density_tot that just includes the active zone
CelloArray<enzo_float,3> density_tot_az = density_tot_arr.subarray
(CSlice(gz, mz - gz), CSlice(gy, my - gy), CSlice(gx, mx - gx));

for (int iz=0; iz<nz; iz++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=0; ix<nx; ix++) {
density_tot_az(iz,iy,ix) += deposited_gas_density(iz,iy,ix);
}
}
}
}

enzo_float hxf = hx * cosmo_a;
enzo_float hyf = hy * cosmo_a;
enzo_float hzf = hz * cosmo_a;
enzo_float dtf = alpha_;
}

enzo_float * vxf = (enzo_float *) field.values("velocity_x");
enzo_float * vyf = (enzo_float *) field.values("velocity_y");
enzo_float * vzf = (enzo_float *) field.values("velocity_z");
//----------------------------------------------------------------------

void EnzoMethodPmDeposit::compute ( Block * block) throw()
{

if (enzo::simulation()->cycle() == enzo::config()->initial_cycle) {
// Check if the gravity method is being used and that pm_deposit
// precedes the gravity method.
ASSERT("EnzoMethodPmDeposit",
"Error: pm_deposit method must precede gravity method.",
enzo::problem()->method_precedes("pm_deposit", "gravity"));
}

enzo_float * vx = new enzo_float [m];
enzo_float * vy = new enzo_float [m];
enzo_float * vz = new enzo_float [m];
if (block->is_leaf()) {

for (int i=0; i<m; i++) vx[i] = vxf[i];
Field field (block->data()->field());

if (rank >= 2) for (int i=0; i<m; i++) vy[i] = vyf[i];
else for (int i=0; i<m; i++) vy[i] = 0.0;
CelloArray<enzo_float,3> density_tot_arr =
field.view<enzo_float>("density_total");
CelloArray<enzo_float,3> density_particle_arr =
field.view<enzo_float>("density_particle");
CelloArray<enzo_float,3> density_particle_accum_arr =
field.view<enzo_float>("density_particle_accumulate");

if (rank >= 3) for (int i=0; i<m; i++) vz[i] = vzf[i];
else for (int i=0; i<m; i++) vz[i] = 0.0;
int mx,my,mz;
field.dimensions(0,&mx,&my,&mz);
int gx,gy,gz;
field.ghost_depth(0,&gx,&gy,&gz);

const int m = mx*my*mz;
std::fill_n(density_tot_arr.data(), mx*my*mz, 0.0);

FORTRAN_NAME(dep_grid_cic)(de,de_gas,temp,
vx, vy, vz,
&dtf, rfield, &rank,
&hxf,&hyf,&hzf,
&mx,&my,&mz,
&gxi,&gyi,&gzi,
&nxi,&nyi,&nzi,
&i0,&i0,&i0,
&nx,&ny,&nz,
&i1,&i1,&i1);
// NOTE 2022-06-24: previously, we filled density_particle_accum_arr with
// zeros at around this location and included the following note:
// NOTE: density_total is now cleared in EnzoMethodGravity to
// instead of here to possible race conditions with refresh. This
// means EnzoMethodPmDeposit ("pm_deposit") currently CANNOT be
// used without EnzoMethodGravity ("gravity")
// This operation & comment didn't sense since we completely overwrite
// values of density_total & density_particle_accum_arr later in this method

for (int i=0; i<mx*my*mz; i++) {
de_t[i] = de_pa[i] = de_p[i];
// Get cosmological scale factors, if cosmology is turned on
enzo_float cosmo_a=1.0;
enzo_float cosmo_dadt=0.0;
EnzoPhysicsCosmology * cosmology = enzo::cosmology();
if (cosmology) {
cosmology->compute_expansion_factor(&cosmo_a,&cosmo_dadt,
block->time() + alpha_*block->dt());
}
const double dt_div_cosmoa = alpha_ * block->dt() / cosmo_a;

for (int iz=gz; iz<mz-gz; iz++) {
for (int iy=gy; iy<my-gy; iy++) {
for (int ix=gx; ix<mx-gx; ix++) {
int i = ix + mx*(iy + my*iz);
int ig = (ix-gx) + nx*((iy-gy) + ny*(iz-gz));
de_t[i] += de_gas[ig];
}
}
double hx,hy,hz;
block->cell_width(&hx,&hy,&hz);

// add mass from particles
{
int rank = cello::rank();
// To calculate densities from particles with "mass" attributes
// or constants, we need the inverse volume of cells in this block.
double inv_vol = 1.0;
if (rank >= 1) inv_vol /= hx;
if (rank >= 2) inv_vol /= hy;
if (rank >= 3) inv_vol /= hz;

deposit_particles_(density_particle_arr, block, dt_div_cosmoa, inv_vol,
mx, my, mz,
gx, gy, gz);

// update density_tot_arr
density_particle_arr.copy_to(density_tot_arr);
density_particle_arr.copy_to(density_particle_accum_arr);
}

delete [] rfield;
delete [] temp;
delete [] vx;
delete [] vy;
delete [] vz;

delete [] de_gas;
// add mass from gas
{
// Grid_DepositBaryons.C from enzo-dev overrides the drift time for the
// gas density to be zero when using the PPM and Zeus solvers.
//
// The way that the Gravity source terms are currently included with
// EnzoMethodPpm and EnzoMethodMHDVlct is currently very similar (in
// terms of temporal integration), so we also set this to zero.
const double gas_dt_div_cosmoa = 0.0;

// The use of "proper" cell-widths was carried over for consistency with
// earlier versions of the code. Based on Grid_DepositBaryons.C from
// enzo-dev, it seems like this may not be correct.
deposit_gas_(density_tot_arr, field, gas_dt_div_cosmoa,
hx*cosmo_a, hy*cosmo_a, hz*cosmo_a,
mx, my, mz,
gx, gy, gz);
}
}

block->compute_done();
Expand Down