add accel derived field output

inputs/precipitator_lowres.in

@@ -10,20 +10,15 @@ dt         = 50.0       # time increment between outputs
 
 <parthenon/output2>
 file_type  = hdf5      # HDF5 data dump
-variables  = prim, drho_over_rho, dP_over_P, dK_over_K, dT_over_T, entropy, temperature, tcool_over_tff, dv_x, dv_y, dv_z
-dt         = 500.0        # Time increment between outputs
+variables  = prim, drho_over_rho, dP_over_P, dK_over_K, dT_over_T, entropy, temperature, tcool_over_tff, dv_x, dv_y, dv_z, accel_x, accel_y, accel_z
+dt         = 100.0        # Time increment between outputs
 id         = prim      # Name to append to output
 
 <parthenon/output3>
 file_type  = rst        # Binary data dump
 dt         = 2000.0        # time increment between outputs
 id         = restart
 
-<parthenon/output4>
-file_type    = ascent
-dt           = 20.0
-actions_file = inputs/precipitator_no_binning.yaml
-
 <parthenon/time>
 cfl        = 0.2        # The Courant, Friedrichs, & Lewy (CFL) Number
 nlim       = -1         # cycle limit
@@ -112,5 +107,5 @@ num_modes = 12
 vertical_driving_only = true
 sol_weight = 1.0
 rseed = 42
-t_corr = 500.0   # Myr
+t_corr = 50.0   # Myr
 max_height = 80.0 # kpc

src/pgen/precipitator.cpp

@@ -326,6 +326,9 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re
     auto v_norm = std::sqrt(v2_sum / (Lx * Ly * Lz) / (SQR(sigma_v)));
 
     auto turbHeat_pack = md->PackVariables(std::vector<std::string>{"turbulent_heating"});
+    auto accel_x_pack = md->PackVariables(std::vector<std::string>{"accel_x"});
+    auto accel_y_pack = md->PackVariables(std::vector<std::string>{"accel_y"});
+    auto accel_z_pack = md->PackVariables(std::vector<std::string>{"accel_z"});
 
     pmb->par_for(
         "apply_perturb_v", 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
@@ -339,6 +342,14 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re
           }
           const Real dv_z = perturb_pack(b, 2, k, j, i) / v_norm;
 
+          // save normalized acceleration field
+          const auto &accel_x = accel_x_pack(b);
+          const auto &accel_y = accel_y_pack(b);
+          const auto &accel_z = accel_z_pack(b);
+          accel_x(0, k, j, i) = dv_x;
+          accel_y(0, k, j, i) = dv_y;
+          accel_z(0, k, j, i) = dv_z;
+
           // compute old kinetic energy
           const auto &u = cons(b);
           const Real rho = u(IDN, k, j, i);
@@ -637,6 +648,17 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
   m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
   pkg->AddField("dv_z", m);
 
+  // add accel_x field
+  m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
+  pkg->AddField("accel_x", m);
+  // add accel_y field
+  m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
+  pkg->AddField("accel_y", m);
+  // add accel_z field
+  m = Metadata({Metadata::Cell, Metadata::OneCopy}, std::vector<int>({1}));
+  pkg->AddField("accel_z", m);
+
+
   const Units units(pin);
   Kokkos::Random_XorShift64_Pool<> random_pool(/*seed=*/12345);