From 9bcced159a4e4a5720d6758fca9f8084574e6a1c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 7 Jul 2022 13:26:22 -0400 Subject: [PATCH 1/2] Fix bug from PR#239 in EnzoMethodPmDeposit PR #239 introduced a minor bug in which array used for storing the output of dep_grid_cic was not large enough. This caused dep_grid_cic to write values to invalid locations in memory. This bug manifested in weird ways when using compiling with intel compilers (under special conditions) In addition to fixing that bug, some variable names were updated so that their significance became more apparent --- src/Enzo/enzo_EnzoMethodPmDeposit.cpp | 61 +++++++++++++++++---------- 1 file changed, 38 insertions(+), 23 deletions(-) diff --git a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp index f9e0157cb6..2c44b86c5a 100644 --- a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp +++ b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp @@ -468,44 +468,59 @@ namespace { // define local helper functions in anonymous namespace if (rank < 2) vyf = vel_scratch.data(); if (rank < 3) vzf = vel_scratch.data() + m; - // deposited_gas_density is a temporary array that just includes cells in - // the active zone - const CelloArray 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); + // allocate memory for storing the deposited gas density. The treatment of + // this array is a little odd: + // - After calling dep_grid_cic, we treat this memory as though just + // enough space has been allocated for holding data in the active + // zone (as though it holds just nz*ny*nx entries). Note: this is + // mainly known from earlier implemenations of this code. + // - It seems as though this array MUST have an allocation for extra + // array allocations. If we allocate just nz*ny*nx entries, then + // dep_grid_cic ends up writing to invalid memory locations. Earlier + // versions of the code allocate mz*my*mx entries which seems + // sufficient. + std::vector deposited_gas_density(m, 0.0); // allocate temporary arrays - std::vector temp(4*m, 0.0); - std::vector rfield(m, 0.0); - - int gxi=gx; - int gyi=gy; - int gzi=gz; - int nxi=mx-gx-1; - int nyi=my-gy-1; - int nzi=mz-gz-1; - int i0 = 0; - int i1 = 1; + std::vector temp(4*m, 0.0); // scratch space + std::vector rfield(m, 0.0); // legacy argument relevant for + // patch-based AMR + + // specify the indices corresponding to the start and end (inclusive) of + // the source fields (namely de, vxf, vyf, vzf) + int src_starti_x=gx; + int src_starti_y=gy; + int src_starti_z=gz; + int src_endi_x=mx-gx-1; + int src_endi_y=my-gy-1; + int src_endi_z=mz-gz-1; + + // dummy args + int offset_val = 0; + int refine_factor = 1; 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, + &src_starti_x, &src_starti_y, &src_starti_z, + &src_endi_x, &src_endi_y, &src_endi_z, + &offset_val, &offset_val, &offset_val, &nx,&ny,&nz, - &i1,&i1,&i1); + &refine_factor, &refine_factor, &refine_factor); - // build a slice of density_tot that just includes the active zone - CelloArray density_tot_az = density_tot_arr.subarray + // construct a subarray density_tot that just includes the active zone + CelloArray density_tot_active_zone = density_tot_arr.subarray (CSlice(gz, mz - gz), CSlice(gy, my - gy), CSlice(gx, mx - gx)); for (int iz=0; iz Date: Thu, 7 Jul 2022 19:44:42 -0400 Subject: [PATCH 2/2] fixing another bug that arose because I zeroed out the wrong array. --- src/Enzo/enzo_EnzoMethodPmDeposit.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp index 2c44b86c5a..962c27ecca 100644 --- a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp +++ b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp @@ -557,8 +557,7 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw() 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); + std::fill_n(density_particle_arr.data(), mx*my*mz, 0.0); // NOTE 2022-06-24: previously, we filled density_particle_accum_arr with // zeros at around this location and included the following note: