print injected cell details

print only data for injected cells add user control for modifying sigma, print injected cell data add Barrier and new var for sigma0 remove vismf write add comments, cleanup, and add user-defined var for ubound of rel_diff, and lbound of ndens for magnetization comp remove print statements ifdef pulsar for isp_nd in fulldiags fix tab fix eol add pulsar input fix comments
RevathiJambunathan · May 19, 2022 · 48d0e4a · 48d0e4a
1 parent 991e0b9
commit 48d0e4a
Show file tree

Hide file tree

Showing 8 changed files with 431 additions and 99 deletions.
diff --git a/Source/Diagnostics/ComputeDiagFunctors/NumberDensityFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/NumberDensityFunctor.H
@@ -19,7 +19,7 @@ public:
      * \param[in] crse_ratio for interpolating field values from the simulation MultiFab, src_mf,
                   to the output diagnostic MultiFab, mf_dst.
      * \param[in] species_index index of species for which number density is computed
-     * \param[in] scomp starting component of mf_src 
+     * \param[in] scomp starting component of mf_src
      * \param[in] ncomp Number of component of mf_src to cell-center in dst multifab.
      */
     NumberDensityFunctor(const amrex::MultiFab * const mf_src, const int lev,

diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp
@@ -154,7 +154,7 @@ FlushFormatCheckpoint::WriteToFile (
                     amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "pml_rz"));
             }
 #endif
-        } 
+        }
 #ifdef PULSAR
         // write magnetization multifab
         VisMF::Write(warpx.getPulsar().get_magnetization(lev), amrex::MultiFabFileFullPrefix(lev, checkpointname, default_level_prefix, "magnetization"));

diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
@@ -585,9 +585,10 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
 
     // Species index to loop over species that dump rho per species
     int i = 0;
+#ifdef PULSAR
     // Species index to loop over species that dump number density per species
     int isp_nd = 0;
-
+#endif
     const auto nvar = static_cast<int>(m_varnames_fields.size());
     const auto nspec = static_cast<int>(m_pfield_species.size());
     const auto ntot = static_cast<int>(nvar + m_pfield_varnames.size() * nspec);

diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
@@ -122,7 +122,6 @@ WarpX::Evolve (int numsteps)
         }
 
 #ifdef PULSAR
-        m_pulsar->TotalParticles();
         m_pulsar->ComputePlasmaNumberDensity();
         m_pulsar->ComputePlasmaMagnetization();
         // inject particles for pulsar simulation

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
@@ -755,6 +755,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
     amrex::Real Sigma0_threshold = Pulsar::m_Sigma0_threshold;
     const MultiFab& magnetization_mf = WarpX::GetInstance().getPulsar().get_magnetization(lev);
     amrex::MultiFab* injection_flag_mf = WarpX::GetInstance().getPulsar().get_pointer_injection_flag(lev);
+    const int modify_sigma_threshold = Pulsar::modify_sigma_threshold;
 #endif
 
     const auto dx = geom.CellSizeArray();
@@ -903,8 +904,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
                                                     pulsar_dR_star*buffer_factor) )
             {
                 // compute threshold magnetization Sigma = Sigma0 * (Rstar/r)^3
-                amrex::Real Sigma_threshold = Sigma0_threshold * (Rstar/rad) * (Rstar/rad) * (Rstar/rad);
-
+                amrex::Real Sigma_threshold = Sigma0_threshold;
+                if (modify_sigma_threshold == 1) {
+                    Sigma_threshold = Sigma0_threshold * (Rstar/rad) * (Rstar/rad) * (Rstar/rad);
+                }
+
                 // inject particles if magnetization sigma > threshold magnetization Sigma
                 if (mag(lo_tile_index[0] + i, lo_tile_index[1] + j, lo_tile_index[2] + k) > Sigma_threshold )
                 {
@@ -913,7 +917,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
                         // instead of modiying number of particles, the weight is changed
                         pcounts[index] = num_ppc;
                     } else if (pulsar_modifyParticleWtAtInjection == 0) {
-		        const amrex::XDim3 ppc_per_dim = inj_pos->getppcInEachDim();
+                const amrex::XDim3 ppc_per_dim = inj_pos->getppcInEachDim();
                         // Modiying number of particles injected
                         // (could lead to round-off errors)
                         pcounts[index] = static_cast<int>(ppc_per_dim.x*std::cbrt(pulsar_injection_fraction))

diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H
@@ -42,6 +42,7 @@ public:
     void TuneSigma0Threshold (const int step);
     void TotalParticles ();
     amrex::Real SumInjectionFlag ();
+    void PrintInjectedCellValues ();
 
 
     AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -518,6 +519,16 @@ public:
     std::list<amrex::Real> ROI_list;
     int list_size = 0;
     static int ROI_avg_window_size;
+    static int modify_sigma_threshold;
+    static int m_print_injected_celldata;
+    static int m_print_celldata_starttime;
+    /** Lower bound for number density to compute magnetization (Default is 1.e-16)*/
+    static amrex::Real m_lbound_ndens_magnetization;
+    /** Upeer bound for relative difference sigma tune method, such that,
+     ** if relative_difference > upper bound, then sigma is modified by upper bound.
+     ** Upper bound specified must be less than 1, and default is 0.1 (10%)
+     */
+    static amrex::Real m_ubound_reldiff_sigma0;
 private:
 };