diff --git a/input/Gravity/jeans_wave/initial_jeans.incl b/input/Gravity/jeans_wave/initial_jeans.incl index 0d897dd4a6..d94a488e74 100644 --- a/input/Gravity/jeans_wave/initial_jeans.incl +++ b/input/Gravity/jeans_wave/initial_jeans.incl @@ -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. diff --git a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp index 20105a49ce..f9e0157cb6 100644 --- a/src/Enzo/enzo_EnzoMethodPmDeposit.cpp +++ b/src/Enzo/enzo_EnzoMethodPmDeposit.cpp @@ -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& 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(); @@ -239,7 +214,7 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw() #endif for (int ip=0; ip& 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 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 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 temp(4*m, 0.0); + std::vector rfield(m, 0.0); int gxi=gx; int gyi=gy; @@ -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 density_tot_az = density_tot_arr.subarray + (CSlice(gz, mz - gz), CSlice(gy, my - gy), CSlice(gx, mx - gx)); + + for (int iz=0; izcycle() == 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; idata()->field()); - if (rank >= 2) for (int i=0; i density_tot_arr = + field.view("density_total"); + CelloArray density_particle_arr = + field.view("density_particle"); + CelloArray density_particle_accum_arr = + field.view("density_particle_accumulate"); - if (rank >= 3) for (int i=0; icompute_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; izcell_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();