From 48d0e4aaa2c00e90c29f7c218d4f5829111fa5ea Mon Sep 17 00:00:00 2001 From: RevathiJambunathan Date: Mon, 9 May 2022 15:46:26 -0700 Subject: [PATCH] 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 --- .../NumberDensityFunctor.H | 2 +- .../FlushFormats/FlushFormatCheckpoint.cpp | 2 +- Source/Diagnostics/FullDiagnostics.cpp | 3 +- Source/Evolve/WarpXEvolve.cpp | 1 - .../Particles/PhysicalParticleContainer.cpp | 10 +- Source/Particles/PulsarParameters.H | 11 + Source/Particles/PulsarParameters.cpp | 238 ++++++++++------ Tools/pulsar_input/inputs.corotating.3d.PM | 263 ++++++++++++++++++ 8 files changed, 431 insertions(+), 99 deletions(-) create mode 100644 Tools/pulsar_input/inputs.corotating.3d.PM diff --git a/Source/Diagnostics/ComputeDiagFunctors/NumberDensityFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/NumberDensityFunctor.H index dff533810c9..b7369f43d47 100644 --- 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 index 6f44cf1a756..3cea6d0aa75 100644 --- 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 index fabbf15fb72..0dbd2c50e3b 100644 --- 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(m_varnames_fields.size()); const auto nspec = static_cast(m_pfield_species.size()); const auto ntot = static_cast(nvar + m_pfield_varnames.size() * nspec); diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b25a7108cc1..c2d00b5de66 100644 --- 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 index fbe7b0da405..cfba3cb7541 100644 --- 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(ppc_per_dim.x*std::cbrt(pulsar_injection_fraction)) diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H index f1fbc281091..8c94ff9b130 100644 --- 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 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: }; diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index 834841f417e..83f17c29791 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -86,6 +86,11 @@ amrex::Real Pulsar::m_max_Sigma0; amrex::Real Pulsar::m_sum_injection_rate = 0.; std::string Pulsar::m_sigma_tune_method; int Pulsar::ROI_avg_window_size = 50; +int Pulsar::modify_sigma_threshold = 0.; +int Pulsar::m_print_injected_celldata; +int Pulsar::m_print_celldata_starttime; +amrex::Real Pulsar::m_lbound_ndens_magnetization = 1.e-16; +amrex::Real Pulsar::m_ubound_reldiff_sigma0 = 0.1; Pulsar::Pulsar () @@ -95,7 +100,6 @@ Pulsar::Pulsar () void Pulsar::ReadParameters () { - amrex::Print() << " pulsar read data \n"; amrex::ParmParse pp("pulsar"); pp.query("pulsarType",m_pulsar_type); @@ -239,6 +243,13 @@ Pulsar::ReadParameters () { pp.get("sigma_tune_method",m_sigma_tune_method); amrex::Print() << " sigma tune method " << m_sigma_tune_method << "\n"; pp.get("ROI_avg_size",ROI_avg_window_size); + pp.get("modify_sigma_threshold", modify_sigma_threshold); + pp.get("print_injected_celldata", m_print_injected_celldata); + pp.get("print_celldata_starttime", m_print_celldata_starttime); + pp.query("lowerBound_ndens_magnetization", m_lbound_ndens_magnetization); + if (m_sigma_tune_method == "relative_difference") { + pp.query("upperBound_reldiff_sigma0", m_ubound_reldiff_sigma0); + } } @@ -298,7 +309,6 @@ Pulsar::InitDataAtRestart () void Pulsar::InitData () { - amrex::Print() << " pulsar init data \n"; auto & warpx = WarpX::GetInstance(); const int nlevs_max = warpx.finestLevel() + 1; //allocate number density multifab @@ -1061,9 +1071,6 @@ Pulsar::ComputePlasmaNumberDensity () get_particle_position(p, xw, yw, zw); // Get position in AMReX convention to calculate corresponding index. - // Ideally this will be replaced with the AMReX NGP interpolator - // Always do x direction. No RZ case because it's not implemented, and code - // will have aborted int ii = 0, jj = 0, kk = 0; amrex::ParticleReal x = p.pos(0); amrex::Real lx = (x - plo[0]) * dxi[0]; @@ -1085,7 +1092,7 @@ Pulsar::ComputePlasmaNumberDensity () ); } - + const Geometry& geom = warpx.Geom(lev); const auto dx = geom.CellSizeArray(); #if defined WARPX_DIM_3D @@ -1117,14 +1124,14 @@ Pulsar::ComputePlasmaMagnetization () const int nlevs_max = warpx.finestLevel() + 1; std::vector species_names = warpx.GetPartContainer().GetSpeciesNames(); const int nspecies = species_names.size(); - const amrex::Real min_ndens = 1.e-16; // make this user-defined + const amrex::Real min_ndens = m_lbound_ndens_magnetization; constexpr amrex::Real mu0_m_c2_inv = 1._rt / (PhysConst::mu0 * PhysConst::m_e * PhysConst::c * PhysConst::c); for (int lev = 0; lev < nlevs_max; ++lev) { m_magnetization[lev]->setVal(0._rt); - const amrex::MultiFab& Bx_mf = warpx.GetInstance().getBfield(lev, 0); - const amrex::MultiFab& By_mf = warpx.GetInstance().getBfield(lev, 1); - const amrex::MultiFab& Bz_mf = warpx.GetInstance().getBfield(lev, 2); + const amrex::MultiFab& Bx_mf = warpx.getBfield(lev, 0); + const amrex::MultiFab& By_mf = warpx.getBfield(lev, 1); + const amrex::MultiFab& Bz_mf = warpx.getBfield(lev, 2); for (MFIter mfi(*m_magnetization[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { const amrex::Box& bx = mfi.tilebox(); @@ -1161,10 +1168,6 @@ Pulsar::ComputePlasmaMagnetization () ); } // mfiter loop } // loop over levels - std::stringstream ss; - int myproc = ParallelDescriptor::MyProc(); - ss << "/diags/sigma" << Concatenate("_",warpx.getistep(0),5)<< Concatenate("_",myproc,5); - VisMF::Write( *m_magnetization[0], ss.str()); } void @@ -1177,6 +1180,12 @@ Pulsar::TuneSigma0Threshold (const int step) amrex::Real dt = warpx.getdt(0); // Total number of cells that have injected particles amrex::Real total_injected_cells = SumInjectionFlag(); + // for debugging print injected cell data + if (m_print_injected_celldata == 1) { + if (warpx.getistep(0) >= m_print_celldata_starttime) { + PrintInjectedCellValues(); + } + } using PTDType = typename WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; for (int isp = 0; isp < nspecies; ++isp) { @@ -1202,73 +1211,71 @@ Pulsar::TuneSigma0Threshold (const int step) reduce_ops ); ws_total = amrex::get<0>(ws_r); - amrex::ParallelDescriptor::ReduceRealSum(ws_total, ParallelDescriptor::IOProcessorNumber()); + amrex::ParallelDescriptor::ReduceRealSum(ws_total); total_weight_allspecies += ws_total; } + // injection rate is sum of particle weight over all species per timestep amrex::Real current_injection_rate = total_weight_allspecies / dt; if (list_size < ROI_avg_window_size) { ROI_list.push_back(current_injection_rate); m_sum_injection_rate += current_injection_rate; list_size++; } else { - amrex::Print() << " front : " << ROI_list.front() << "\n"; m_sum_injection_rate -= ROI_list.front(); ROI_list.pop_front(); ROI_list.push_back(current_injection_rate); - amrex::Print() << " current injection rate : " << current_injection_rate << " list back " << ROI_list.back() << "\n"; m_sum_injection_rate += ROI_list.back(); } amrex::Print() << " current_injection rate " << current_injection_rate << " sum : " << m_sum_injection_rate << "\n"; amrex::Real specified_injection_rate = m_GJ_injection_rate * m_injection_rate; if (m_injection_tuning_interval.contains(step+1) ) { - amrex::Print() << " period for injection tuning : " << m_injection_tuning_interval.localPeriod(step+1) << "\n"; - amrex::Print() << " species rate is : " << specified_injection_rate << " current rate : " << current_injection_rate << "\n"; - amrex::Print() << " Sigma0 before mod : " << m_Sigma0_threshold << "\n"; + // Sigma0 before modification amrex::Real m_Sigma0_pre = m_Sigma0_threshold; -// amrex::Real avg_injection_rate = m_sum_injection_rate/(m_injection_tuning_interval.localPeriod(step+1) * 1.); - amrex::Print() << " list size : " << list_size << "\n"; - amrex::Real avg_injection_rate = m_sum_injection_rate/(list_size * 1.); - amrex::Print() << " avg inj rate " << avg_injection_rate << "\n"; + // running average of injection rate over average window + amrex::Real avg_injection_rate = m_sum_injection_rate/(list_size * 1._rt); + amrex::Real new_sigma0_threshold = m_Sigma0_threshold; if (avg_injection_rate < specified_injection_rate) { // reduce sigma0 so more particles can be injected - //m_Sigma0_threshold *= total_weight_allspecies/specified_injection_rate; if (m_sigma_tune_method == "10percent") { - m_Sigma0_threshold = m_Sigma0_threshold - 0.1 * m_Sigma0_threshold; + // Modify threshold magnetization threshold by 10% + new_sigma0_threshold = m_Sigma0_threshold - 0.1 * m_Sigma0_threshold; } if (m_sigma_tune_method == "relative_difference") { if (avg_injection_rate == 0.) { - // if no particles are injection only change the threshold sigma by 10% - m_Sigma0_threshold = m_Sigma0_threshold - 0.1 * m_Sigma0_threshold; + // If no particles are injected only change the threshold sigma by 10% + new_sigma0_threshold = m_Sigma0_threshold - 0.1 * m_Sigma0_threshold; } else { amrex::Real rel_diff = (specified_injection_rate - avg_injection_rate)/specified_injection_rate; - if (rel_diff < 0.1) { - amrex::Print() << " rel_diff " << rel_diff << "\n"; - m_Sigma0_threshold = m_Sigma0_threshold - rel_diff * m_Sigma0_threshold; + if (rel_diff < m_ubound_reldiff_sigma0) { + // If relative difference is less than user-defined upper bound, reduce magnetization by rel_diff + new_sigma0_threshold = m_Sigma0_threshold - rel_diff * m_Sigma0_threshold; } else { - amrex::Print() << " rel_diff " << rel_diff << " using upper bound 10%"<< "\n"; - m_Sigma0_threshold = m_Sigma0_threshold - 0.1 * m_Sigma0_threshold; + // Maximum relative difference is set by user-defined upper bound + new_sigma0_threshold = m_Sigma0_threshold - m_ubound_reldiff_sigma0 * m_Sigma0_threshold; } } } } else if (avg_injection_rate > specified_injection_rate ) { - //m_Sigma0_threshold *= specified_injection_rate/total_weight_allspecies; + // Increase sigma0 so fewer particles are injected if (m_sigma_tune_method == "10percent") { - m_Sigma0_threshold = m_Sigma0_threshold + 0.1 * m_Sigma0_threshold; + // Modify threshold magnetization threshold by 10% + new_sigma0_threshold = m_Sigma0_threshold + 0.1 * m_Sigma0_threshold; } if (m_sigma_tune_method == "relative_difference") { amrex::Real rel_diff = (avg_injection_rate - specified_injection_rate)/specified_injection_rate; - if (rel_diff < 0.1) { - amrex::Print() << " rel_diff " << rel_diff << "\n"; - m_Sigma0_threshold = m_Sigma0_threshold + rel_diff * m_Sigma0_threshold; + if (rel_diff < m_ubound_reldiff_sigma0) { + // If relative difference is less than user-defined upper bound, increase magnetization by rel_diff + new_sigma0_threshold = m_Sigma0_threshold + rel_diff * m_Sigma0_threshold; } else { - amrex::Print() << " rel_diff " << rel_diff << " using upper bound 10%"<< "\n"; - m_Sigma0_threshold = m_Sigma0_threshold + 0.1 * m_Sigma0_threshold; + // If rel_diff > upper bound, then increase sigma0 by m_ubound_reldiff_sigma0 + new_sigma0_threshold = m_Sigma0_threshold + m_ubound_reldiff_sigma0 * m_Sigma0_threshold; } } } - if (m_Sigma0_threshold < m_min_Sigma0) m_Sigma0_threshold = m_min_Sigma0; - if (m_Sigma0_threshold > m_max_Sigma0) m_Sigma0_threshold = m_max_Sigma0; - amrex::Print() << " Simg0 modified to : " << m_Sigma0_threshold << "\n"; + if (new_sigma0_threshold < m_min_Sigma0) new_sigma0_threshold = m_min_Sigma0; + if (new_sigma0_threshold > m_max_Sigma0) new_sigma0_threshold = m_max_Sigma0; + // Store modified new sigma0 in member variable, m_Sigma0_threshold + m_Sigma0_threshold = new_sigma0_threshold; amrex::AllPrintToFile("RateOfInjection") << warpx.getistep(0) << " " << warpx.gett_new(0) << " " << dt << " " << specified_injection_rate << " " << avg_injection_rate << " " << m_Sigma0_pre << " "<< m_Sigma0_threshold << " " << m_min_Sigma0 << " " << m_max_Sigma0 << " " << m_Sigma0_baseline<< "\n"; } amrex::AllPrintToFile("ROI") << warpx.getistep(0) << " " << warpx.gett_new(0) << " " << dt << " " << specified_injection_rate << " " << current_injection_rate << " "<< m_Sigma0_threshold << " " << total_injected_cells << "\n"; @@ -1280,7 +1287,7 @@ Pulsar::TotalParticles () auto& warpx = WarpX::GetInstance(); std::vector species_names = warpx.GetPartContainer().GetSpeciesNames(); const int nspecies = species_names.size(); - amrex::Long total_particles = 0; + amrex::Long total_particles = 0; for (int isp = 0; isp < nspecies; ++isp) { auto& pc = warpx.GetPartContainer().GetParticleContainer(isp); @@ -1288,7 +1295,6 @@ Pulsar::TotalParticles () amrex::ParallelDescriptor::ReduceLongSum(np_total, ParallelDescriptor::IOProcessorNumber()); total_particles += np_total; } - amrex::Print() << " total particles " << total_particles << "\n"; } amrex::Real @@ -1297,50 +1303,98 @@ Pulsar::SumInjectionFlag () return m_injection_flag[0]->sum(); } -//void -//Pulsar::PrintInjectedCellValues () -//{ -// auto& warpx = WarpX::GetInstance(); -// std::vector species_names = warpx.GetPartContainer();GetSpeciesNames(); -// const int nspecies = species_names.size(); -// -// int total_injected_cells = SumInjectionFlag(); -// // x, y, z, r, theta, phi, injection_flag, magnetization, ndens_p, ndens_e, Bx, By, Bz, Bmag -// int total_diags = 14; -// amrex::Gpu::DeviceVector InjectedCellDiag(total_diags,0.0); -// const int lev = 0; -// const auto dx_lev = warpx.Geom(lev).CellSizeArray(); -// const amrex::RealBox* real_box = warpx.Geom(lev).ProbDomain(); -// const amrex::MultiFab& Bx_mf = warpx.getBfield(lev,0); -// const amrex::MultiFab& By_mf = warpx.getBfield(lev,1); -// const amrex::MultiFab& Bz_mf = warpx.getBfield(lev,2); -// const amrex::MultiFab& injectionflag_mf = *m_injection_flag[lev]; -// const amrex::MultiFab& magnetization_mf = *m_magnetization[lev]; -// const amrex::MultiFab& ndens_mf = *m_plasma_number_density[lev]; -// int cell_counter = 0; -// for (MFIter mfi(injectionflag_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) -// { -// const amrex::Box & bx = mfi.tilebox(); -// amrex::Array4 const& Bx = Bx_mf[mfi].array(); -// amrex::Array4 const& By = By_mf[mfi].array(); -// amrex::Array4 const& Bz = Bz_mf[mfi].array(); -// amrex::Array4 const& ndens = ndens_mf[mfi].array(); -// amrex::Array4 const& injection = injectionflag_mf[mfi].array(); -// amrex::Array4 const& mag = magnetization_mf[mfi].array(); -// amrex::LoopOnCpu( bx, -// [=] (int i, int j, int k) -// { -// if (injection(i,j,k) == 1) { -// amrex::Real Bx_cc = (Bx(i,j,k) + Bx(i+1, j, k)) / 2.0; -// amrex::Real By_cc = (By(i,j,k) + By(i, j+1, k)) / 2.0; -// amrex::Real Bz_cc = (Bz(i,j,k) + Bz(i, j, k+1)) / 2.0; -// // magnitude of magnetic field -// amrex::Real B_mag = std::sqrt( Bx_cc * Bx_cc -// + By_cc * By_cc -// + Bz_cc * Bz_cc ); -// InjectedCellDiag[0] -// counter++; -// } -// }); -// } -//} +void +Pulsar::PrintInjectedCellValues () +{ + auto& warpx = WarpX::GetInstance(); + std::vector species_names = warpx.GetPartContainer().GetSpeciesNames(); + const int nspecies = species_names.size(); + int total_injected_cells = static_cast(SumInjectionFlag()); + // x, y, z, r, theta, phi, injection_flag, magnetization, ndens_p, ndens_e, Bx, By, Bz, Bmag, rho + int total_diags = 15; + amrex::Gpu::DeviceVector InjectedCellDiag(total_injected_cells*total_diags,0.0); + amrex::Real * InjectedCellDiagData = InjectedCellDiag.data(); + const int lev = 0; + const auto dx_lev = warpx.Geom(lev).CellSizeArray(); + const amrex::RealBox real_box = warpx.Geom(lev).ProbDomain(); + const amrex::MultiFab& Bx_mf = warpx.getBfield(lev,0); + const amrex::MultiFab& By_mf = warpx.getBfield(lev,1); + const amrex::MultiFab& Bz_mf = warpx.getBfield(lev,2); + const amrex::MultiFab& injectionflag_mf = *m_injection_flag[lev]; + const amrex::MultiFab& magnetization_mf = *m_magnetization[lev]; + const amrex::MultiFab& ndens_mf = *m_plasma_number_density[lev]; + GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + center_star_arr[idim] = m_center_star[idim]; + } + std::unique_ptr rho; + auto& mypc = warpx.GetPartContainer(); + rho = mypc.GetChargeDensity(lev, true); + amrex::MultiFab & rho_mf = *rho; + + Gpu::DeviceScalar cell_counter(0); + int* cell_counter_d = cell_counter.dataPtr(); + for (MFIter mfi(injectionflag_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box & bx = mfi.tilebox(); + amrex::Array4 const& Bx = Bx_mf[mfi].array(); + amrex::Array4 const& By = By_mf[mfi].array(); + amrex::Array4 const& Bz = Bz_mf[mfi].array(); + amrex::Array4 const& ndens = ndens_mf[mfi].array(); + amrex::Array4 const& injection = injectionflag_mf[mfi].array(); + amrex::Array4 const& mag = magnetization_mf[mfi].array(); + amrex::Array4 const& rho_arr = rho_mf[mfi].array(); + amrex::LoopOnCpu( bx, + [=] (int i, int j, int k) + { + if (injection(i,j,k) == 1) { + // Cartesian Coordinates + amrex::Real x = i * dx_lev[0] + real_box.lo(0) + 0.5_rt * dx_lev[0]; + amrex::Real y = j * dx_lev[1] + real_box.lo(1) + 0.5_rt * dx_lev[1]; + amrex::Real z = k * dx_lev[2] + real_box.lo(2) + 0.5_rt * dx_lev[2]; + // convert cartesian to spherical coordinates + amrex::Real r, theta, phi; + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + amrex::Real Bx_cc = (Bx(i,j,k) + Bx(i+1, j, k)) / 2.0; + amrex::Real By_cc = (By(i,j,k) + By(i, j+1, k)) / 2.0; + amrex::Real Bz_cc = (Bz(i,j,k) + Bz(i, j, k+1)) / 2.0; + // magnitude of magnetic field + amrex::Real B_mag = std::sqrt( Bx_cc * Bx_cc + + By_cc * By_cc + + Bz_cc * Bz_cc ); + int counter = *cell_counter_d; + InjectedCellDiagData[counter*total_diags+0] = x; + InjectedCellDiagData[counter*total_diags+1] = y; + InjectedCellDiagData[counter*total_diags+2] = z; + InjectedCellDiagData[counter*total_diags+3] = r; + InjectedCellDiagData[counter*total_diags+4] = theta; + InjectedCellDiagData[counter*total_diags+5] = phi; + InjectedCellDiagData[counter*total_diags+6] = injection(i, j, k); + InjectedCellDiagData[counter*total_diags+7] = mag(i, j, k); + InjectedCellDiagData[counter*total_diags+8] = ndens(i, j, k, 0); + InjectedCellDiagData[counter*total_diags+9] = ndens(i, j, k, 1); + InjectedCellDiagData[counter*total_diags+10] = Bx_cc; + InjectedCellDiagData[counter*total_diags+11] = By_cc; + InjectedCellDiagData[counter*total_diags+12] = Bz_cc; + InjectedCellDiagData[counter*total_diags+13] = B_mag; + InjectedCellDiagData[counter*total_diags+14] = rho_arr(i, j, k); + const int unity = 1; + amrex::HostDevice::Atomic::Add(cell_counter_d, unity); + } + }); + } + amrex::Print() << " counter : " << cell_counter.dataValue() << " total cells injected " << total_injected_cells << "\n"; + std::stringstream ss; + ss << Concatenate("InjectionCellData", warpx.getistep(0), 5); + amrex::AllPrintToFile(ss.str()) << " cell_index x y z r theta phi injection magnetization ndens_p ndens_e Bx By Bz Bmag rho \n" ; + for (int icell = 0; icell < total_injected_cells; ++icell ) { + if (InjectedCellDiagData[icell*total_diags + 6] == 1) { + amrex::AllPrintToFile(ss.str()) << icell << " "; + for (int idata = 0; idata < total_diags; ++idata) { + amrex::AllPrintToFile(ss.str()) << InjectedCellDiagData[icell * total_diags + idata] << " "; + } + } + amrex::AllPrintToFile(ss.str()) << "\n"; + } +} diff --git a/Tools/pulsar_input/inputs.corotating.3d.PM b/Tools/pulsar_input/inputs.corotating.3d.PM new file mode 100644 index 00000000000..f0867c2dc77 --- /dev/null +++ b/Tools/pulsar_input/inputs.corotating.3d.PM @@ -0,0 +1,263 @@ +# Description: +# +# This inputs file sets up a NS with these properties: +# - E = -(omega*R)[cross]B inside the NS +# - external E and B outside the star +# +# See the "Pulsar Setup" section at the end for the options +# +# This initializes the electrons and positrons with a corotating momentum function. +# Based on the pulsar_type = "active" or "dead", particles are injected continuously or +# until rho_GJ is reached + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 20000 +amr.n_cell = 768 768 768 +amr.max_grid_size = 64 +amr.blocking_factor = 64 +amr.max_level = 0 +geometry.coord_sys = 0 # 0: Cartesian +geometry.dims = 3 +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 180000 180000 180000 +boundary.field_lo = pml pml pml +boundary.field_hi = pml pml pml +#amr.restart = ./diags/chk004500 + +################################# +############ NUMERICS ########### +################################# +algo.maxwell_solver = yee +warpx.verbose = 1 +warpx.do_dive_cleaning = 0 +warpx.use_filter = 1 +warpx.cfl = .95 +warpx.do_pml_in_domain = 1 +warpx.pml_has_particles = 1 +warpx.do_pml_j_damping = 1 +my_constants.pi = 3.141592653589793 +my_constants.dens = 5.544e6 +my_constants.xc = 90000 +my_constants.yc = 90000 +my_constants.zc = 90000 +my_constants.r_star = 12000 +my_constants.omega = 6245.676 +my_constants.c = 299792458. +my_constants.B_star = 8.0323e-6 # magnetic field of NS (T) +my_constants.dR = 234.375 +my_constants.to = 2.e-4 +my_constants.skin_depth = 2000 +algo.particle_shape = 3 + +######################################## +########### load balancing ############# +######################################## +algo.load_balance_intervals = 100 +algo.load_balance_efficiency_ratio_threshold = 1.05 +algo.load_balance_knapsack_factor = 1.2 +algo.costs_heuristic_particles_wt = 0.65 +algo.costs_heuristic_cells_wt = 0.35 +algo.load_balance_costs_update = heuristic +algo.load_balance_with_sfc = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = plasma_e plasma_p + +plasma_e.charge = -q_e +plasma_e.mass = m_e +plasma_e.injection_style = "NUniformPerCell" +plasma_e.profile = parse_density_function +plasma_e.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star + skin_depth)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star+skin_depth-dR)) )*dens" +plasma_e.num_particles_per_cell_each_dim = 2 2 2 +plasma_e.momentum_distribution_type = parse_momentum_function + +## Inject stationary pairs +plasma_e.momentum_function_ux(x,y,z) = "0.0" +plasma_e.momentum_function_uy(x,y,z) = "0.0" +## uz is always 0 for the injection methods above +plasma_e.momentum_function_uz(x,y,z) = "0.0" +plasma_e.do_continuous_injection = 0 +plasma_e.density_min = 0.1 + +plasma_p.charge = q_e +plasma_p.mass = m_e +plasma_p.injection_style = "NUniformPerCell" +plasma_p.profile = parse_density_function +plasma_p.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star+skin_depth)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star+skin_depth-dR)) )*dens" +plasma_p.num_particles_per_cell_each_dim = 2 2 2 +plasma_p.momentum_distribution_type = parse_momentum_function +## Inject stationary pairs +plasma_p.momentum_function_ux(x,y,z) = "0.0" +plasma_p.momentum_function_uy(x,y,z) = "0.0" +## uz ip always 0 for the injection methods above +plasma_p.momentum_function_uz(x,y,z) = "0.0" +plasma_p.do_continuous_injection = 0 +plasma_p.density_min = 0.1 + +#### Adding particle attributes ### +plasma_e.addRealAttributes = injectiontime +plasma_e.attribute.injectiontime(x,y,z,ux,uy,uz,t) = t +plasma_p.addRealAttributes = injectiontime +plasma_p.attribute.injectiontime(x,y,z,ux,uy,uz,t) = t + +diagnostics.diags_names = plt chk +plt.diag_type = Full +plt.intervals = 500 +plt.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz divE conductor rho part_per_cell rho_plasma_e rho_plasma_p magnetization ndens_plasma_p ndens_plasma_e injectionflag + +chk.format = checkpoint +chk.diag_type = Full +chk.intervals = 500 + +############ Reduced Diags ######### +warpx.reduced_diags_names = total_rho_domain +#warpx.reduced_diags_names = divE_conductor divE_domain Ex_xmin Ey_ymin Ez_zmin Ex_xmax Ey_ymax Ez_zmax total_rho_conductor total_rho_domain LB_cost LB_eff PartNum + +LB_cost.type = LoadBalanceCosts +LB_cost.intervals = 50 + +LB_eff.type = LoadBalanceEfficiency +LB_eff.intervals = 50 + +PartNum.type = ParticleNumber +PartNum.intervals = 10 + +total_rho_conductor.type = ParticleNumber +total_rho_conductor.filter_function = " ( (x >= xc-r_star) * (x <= xc+r_star) * (y >= yc-r_star) * (y <= yc+r_star) * (z >= zc-r_star) * (z <= zc+r_star) ) * 1 + 0" +total_rho_conductor.intervals = 1 + +total_rho_domain.type = ParticleNumber +total_rho_domain.filter_function = " (x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc) )^0.5 >= (r_star + skin_depth)" +total_rho_domain.intervals = 500 + +divE_conductor.type = divEReduction +divE_conductor.reduction_type = integral +divE_conductor.integration_type = volume +divE_conductor.reduced_function(x,y,z) = " ( (x >= xc-r_star) * (x <= xc+r_star) * (y >= yc-r_star) * (y <= yc+r_star) * (z >= zc-r_star) * (z <= zc+r_star) ) * 1 + 0" +divE_conductor.intervals = 1 + +divE_domain.type = divEReduction +divE_domain.reduction_type = integral +divE_domain.integration_type = volume +divE_domain.reduced_function(x,y,z) = " 1" +divE_domain.intervals = 1 + +my_constants.bbox = 14702 +my_constants.bbox_mindx = 14655.2 +my_constants.bbox_maxdx = 14877.5 + +Ex_xmin.type = RawEFieldReduction +Ex_xmin.reduction_type = integral +Ex_xmin.integration_type = surface +Ex_xmin.reduced_function(x,y,z) = " ( (x >= xc-bbox_maxdx) * (x <= xc-bbox_mindx) * (y >= yc-bbox) * (y <= yc+bbox) * (z >= zc-bbox) * (z <= zc+bbox) ) * 1 + 0" +Ex_xmin.intervals = 1 + +Ex_xmax.type = RawEFieldReduction +Ex_xmax.reduction_type = integral +Ex_xmax.integration_type = surface +Ex_xmax.reduced_function(x,y,z) = " ( (x >= xc+bbox_mindx) * (x <= xc+bbox_maxdx) * (y >= yc-bbox) * (y <= yc+bbox) * (z >= zc-bbox) * (z <= zc+bbox) ) * 1 + 0" +Ex_xmax.intervals = 1 + +Ey_ymin.type = RawEFieldReduction +Ey_ymin.reduction_type = integral +Ey_ymin.integration_type = surface +Ey_ymin.reduced_function(x,y,z) = " ( (x >= xc-bbox) * (x <= xc+bbox) * (y >= yc-bbox_maxdx) * (y <= yc-bbox_mindx) * (z >= zc-bbox) * (z <= zc+bbox) ) * 1 + 0" +Ey_ymin.intervals = 1 + +Ey_ymax.type = RawEFieldReduction +Ey_ymax.reduction_type = integral +Ey_ymax.integration_type = surface +Ey_ymax.reduced_function(x,y,z) = " ( (x >= xc-bbox) * (x <= xc+bbox) * (y >= yc+bbox_mindx) * (y <= yc+bbox_maxdx) * (z >= zc-bbox) * (z <= zc+bbox) ) * 1 + 0" +Ey_ymax.intervals = 1 + +Ez_zmin.type = RawEFieldReduction +Ez_zmin.reduction_type = integral +Ez_zmin.integration_type = surface +Ez_zmin.reduced_function(x,y,z) = " ( (x >= xc-bbox) * (x <= xc+bbox) * (y >= yc-bbox) * (y <= yc+bbox) * (z >= zc-bbox_maxdx) * (z <= zc-bbox_mindx) ) * 1 + 0" +Ez_zmin.intervals = 1 + +Ez_zmax.type = RawEFieldReduction +Ez_zmax.reduction_type = integral +Ez_zmax.integration_type = surface +Ez_zmax.reduced_function(x,y,z) = " ( (x >= xc-bbox) * (x <= xc+bbox) * (y >= yc-bbox) * (y <= yc+bbox) * (z >= zc+bbox_mindx) * (z <= zc+bbox_maxdx) ) * 1 + 0" +Ez_zmax.intervals = 1 + +################################# +######### PULSAR SETUP ########## +################################# +pulsar.pulsarType = "dead" # [dead/active]: sets particle injection type +pulsar.omega_star = 6245.676 # angular velocity of NS (rad/s) +pulsar.ramp_omega_time = 5e-4 # time over which to ramp up NS angular velocity (s) + # if ramp_omega_time < 0, omega = omega_star for any t + # consistency requires ramp_omega_time = my_constants.to +pulsar.center_star = 90000 90000 90000 +pulsar.R_star = 12000 # radius of NS (m) +pulsar.B_star = 8.0323e-6 # magnetic field of NS (T) +pulsar.dR = 234.375 # thickness on the surface of the pulsar + # consistency requires dR = my_constants.dR +pulsar.verbose = 0 # [0/1]: turn on verbosity for debugging print statements +pulsar.EB_external = 0 # [0/1]: to apply external E and B +pulsar.EB_corotating_maxradius = 14000 # The radius where E-field changes from corotating + # inside the star to quadrapole outside. + # Default is R_star. (r<=cor_radius) i.e. the includes + # the r specified in the input +pulsar.max_ndens = 5.54e6 # max ndens == ndens used in initializing density +pulsar.Ninj_fraction = 1.0 # fraction of sigma injected +pulsar.ModifyParticleWeight = 0 # (0/1) the fractional injection is achieved + # by modifying particle weight if 1 + # by modifying num_ppc if 0 +pulsar.particle_inj_rmin = 13672.6 # default is Rstar-dR (consistent with density function above) +pulsar.particle_inj_rmax = 14000 # default is Rstar (consistent with density function) +pulsar.max_particle_absorption_radius = 12000 # Maximum radius within which particles are + # deleted every timestep. + # Default is Rstar +pulsar.rhoGJ_scale = 1e0 # scaling down of rho_GJ +pulsar.damp_EB_internal = 0 # Damp internal electric field +pulsar.damp_EB_radius = 12000 # The radius within which EB should be damped. + # default is r_star (damping will include this radius) +pulsar.damping_scale = 1000.0 # Damping scale factor for internal electric field +pulsar.turnoffdeposition = 0 # (0/1) 0 for depositing (default) + # 1 for no deposition +pulsar.max_nodepos_radius = 0. # radius within which particle deposition for j/rho + # is turned off. (r<=max_radius) +pulsar.turnoff_plasmaEB_gather = 0 # (0/1) 0 is default. always gather + # 1 for no gather of self-consistent fields +pulsar.max_nogather_radius = 0. # radius within which self-consistent field gather +pulsar.init_dipoleBfield = 1 # default is 1 +pulsar.init_corotatingEfield = 1 # default is 1 +pulsar.init_corotatingAndExternalEField = 0 # default is 0 +pulsar.corotatingE_maxradius = 14000 # default is Rstar +pulsar.enforceDipoleB_maxradius = 12000 # default is Rstar +pulsar.enforceCorotatingE = 1 # default 1 +pulsar.enforceDipoleB = 1 # default 1 +pulsar.singleParticleTest = 0 # default 0, if 1, then a particle pair will be introduced. + # Single particle position/momentum will need to be specified +pulsar.continuous_injection = 0 +pulsar.injection_time = 0. +pulsar.DampBDipoleInRing = 0 +pulsar.Bdamping_scale = 10000 +pulsar.EnforceTheoreticalEBInGrid = 0 +pulsar.E_external_monopole = 0 # [0/1] must be 1 when using use_theoretical = 1 +pulsar.use_theoreticalEB = 0 # theoretical EB is added to particles directly +pulsar.theory_max_rstar = 0 # one cell below rstar +pulsar.AddExternalMonopoleOnly = 1 +pulsar.conductor_function(x,y,z) = " ( ((z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc))^0.5 <= r_star + skin_depth ) * 1 + 0." # sphere +pulsar.ApplyEfieldBCusingConductor = 0 # this sets Efield to 0 in the conductor +pulsar.FilterWithConductor = 0 # uses method of images for filtering. does not work at corners +pulsar.minimum_Sigma0 = 1.e-300 +pulsar.maximum_Sigma0 = 1.e+300 +pulsar.injection_tuning_interval = 10 +pulsar.plasma_injection_rate = 0.2 +pulsar.ROI_avg_size = 50 +pulsar.sigma_tune_method = relative_difference # 10percent or relative_difference +pulsar.modify_sigma_threshold = 1 +pulsar.print_injected_celldata = 0 +pulsar.print_celldata_starttime = 100000000 +# to change lower bound number density for magnetization use the input below +#pulsar.lowerBound_ndens_magnetization = 1.e-16 +# to change upperBound for relative difference method to tune sigma use input below +#pulsar.upperBound_reldiff_sigma0 = 0.1