diff options
Diffstat (limited to 'Source')
-rw-r--r-- | Source/Diagnostics/ReducedDiags/BeamRelevant.cpp | 9 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/CMakeLists.txt | 1 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/FieldEnergy.cpp | 7 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/FieldMaximum.cpp | 7 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp | 8 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/Make.package | 1 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp | 6 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp | 7 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp | 5 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/ParticleNumber.cpp | 13 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/RhoMaximum.H | 40 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/RhoMaximum.cpp | 170 |
12 files changed, 244 insertions, 30 deletions
diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp index 57ddc64c5..fa7ff377f 100644 --- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -37,7 +37,7 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 13: rms gamma // 14,15,16: emittance x,y,z // 17: charge - m_data.resize(18,0.0); + m_data.resize(18, 0.0_rt); #elif (defined WARPX_DIM_XZ) // 0, 1: mean x,z // 2, 3, 4: mean px,py,pz @@ -47,7 +47,7 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 11: rms gamma // 12,13: emittance x,z // 14: charge - m_data.resize(15,0.0); + m_data.resize(15, 0.0_rt); #endif if (ParallelDescriptor::IOProcessor()) @@ -55,8 +55,7 @@ BeamRelevant::BeamRelevant (std::string rd_name) if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row #if (defined WARPX_DIM_3D || defined WARPX_DIM_RZ) ofs << "#"; @@ -161,7 +160,7 @@ void BeamRelevant::ComputeDiags (int step) if (w_sum < std::numeric_limits<Real>::min() ) { for (int i = 0; i < static_cast<int>(m_data.size()); ++i){ - m_data[i] = 0.0; + m_data[i] = 0.0_rt; } return; } diff --git a/Source/Diagnostics/ReducedDiags/CMakeLists.txt b/Source/Diagnostics/ReducedDiags/CMakeLists.txt index 08276f254..070b6815b 100644 --- a/Source/Diagnostics/ReducedDiags/CMakeLists.txt +++ b/Source/Diagnostics/ReducedDiags/CMakeLists.txt @@ -9,5 +9,6 @@ target_sources(WarpX ReducedDiags.cpp FieldMaximum.cpp ParticleExtrema.cpp + RhoMaximum.cpp ParticleNumber.cpp ) diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index 892ea9c7c..cf02e08fe 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -37,15 +37,14 @@ FieldEnergy::FieldEnergy (std::string rd_name) constexpr int noutputs = 3; // total energy, E-field energy and B-field energy // resize data array - m_data.resize(noutputs*nLevel,0.0); + m_data.resize(noutputs*nLevel, 0.0_rt); if (ParallelDescriptor::IOProcessor()) { if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row ofs << "#"; ofs << "[1]step()"; @@ -82,7 +81,7 @@ void FieldEnergy::ComputeDiags (int step) // Judge if the diags should be done if (!m_intervals.contains(step+1)) { return; } - // get WarpX class object + // get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); // get number of level diff --git a/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp b/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp index a7250c5aa..ad2ee1a05 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMaximum.cpp @@ -30,15 +30,14 @@ FieldMaximum::FieldMaximum (std::string rd_name) constexpr int noutputs = 8; // total energy, E-field energy and B-field energy // resize data array - m_data.resize(noutputs*nLevel,0.0); // max of Ex,Ey,Ez,|E|,Bx,By,Bz and |B| + m_data.resize(noutputs*nLevel, 0.0_rt); // max of Ex,Ey,Ez,|E|,Bx,By,Bz and |B| if (ParallelDescriptor::IOProcessor()) { if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row ofs << "#"; ofs << "[1]step()"; @@ -94,7 +93,7 @@ void FieldMaximum::ComputeDiags (int step) // Judge if the diags should be done if (!m_intervals.contains(step+1)) { return; } - // get WarpX class object + // get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); // get number of level diff --git a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp index 048aa7dab..3ec0f0d45 100644 --- a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp +++ b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp @@ -23,7 +23,7 @@ LoadBalanceCosts::LoadBalanceCosts (std::string rd_name) // function that gathers costs void LoadBalanceCosts::ComputeDiags (int step) { - // get WarpX class object + // get a reference to WarpX instance auto& warpx = WarpX::GetInstance(); const amrex::LayoutData<amrex::Real>* cost = warpx.getCosts(0); @@ -50,8 +50,8 @@ void LoadBalanceCosts::ComputeDiags (int step) const size_t dataSize = static_cast<size_t>(m_nDataFields)* static_cast<size_t>(nBoxes); - m_data.resize(dataSize, 0.0); - m_data.assign(dataSize, 0.0); + m_data.resize(dataSize, 0.0_rt); + m_data.assign(dataSize, 0.0_rt); // read in WarpX costs to local copy; compute if using `Heuristic` update amrex::Vector<std::unique_ptr<amrex::LayoutData<amrex::Real> > > costs; @@ -218,7 +218,7 @@ void LoadBalanceCosts::WriteToFile (int step) const ofs.close(); - // get WarpX class object + // get a reference to WarpX instance auto& warpx = WarpX::GetInstance(); if (!ParallelDescriptor::IOProcessor()) return; diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package index 3a8361693..22ba3dd34 100644 --- a/Source/Diagnostics/ReducedDiags/Make.package +++ b/Source/Diagnostics/ReducedDiags/Make.package @@ -7,6 +7,7 @@ CEXE_sources += LoadBalanceCosts.cpp CEXE_sources += ParticleHistogram.cpp CEXE_sources += FieldMaximum.cpp CEXE_sources += ParticleExtrema.cpp +CEXE_sources += RhoMaximum.cpp CEXE_sources += ParticleNumber.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 9c67cea9e..c7ce7b18c 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -12,6 +12,7 @@ #include "ParticleExtrema.H" #include "FieldEnergy.H" #include "FieldMaximum.H" +#include "RhoMaximum.H" #include "ParticleNumber.H" #include "MultiReducedDiags.H" @@ -62,6 +63,11 @@ MultiReducedDiags::MultiReducedDiags () m_multi_rd[i_rd] = std::make_unique<FieldMaximum>(m_rd_names[i_rd]); } + else if (rd_type.compare("RhoMaximum") == 0) + { + m_multi_rd[i_rd] = + std::make_unique<RhoMaximum>(m_rd_names[i_rd]); + } else if (rd_type.compare("BeamRelevant") == 0) { m_multi_rd[i_rd] = diff --git a/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp index 9d2b3eac2..72395cc84 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp @@ -23,7 +23,7 @@ using namespace amrex; ParticleEnergy::ParticleEnergy (std::string rd_name) : ReducedDiags{rd_name} { - // get WarpX class object + // get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); // get MultiParticleContainer class object @@ -33,7 +33,7 @@ ParticleEnergy::ParticleEnergy (std::string rd_name) const auto nSpecies = mypc.nSpecies(); // resize data array - m_data.resize(2*nSpecies+2,0.0); + m_data.resize(2*nSpecies+2, 0.0_rt); // get species names (std::vector<std::string>) const auto species_names = mypc.GetSpeciesNames(); @@ -43,8 +43,7 @@ ParticleEnergy::ParticleEnergy (std::string rd_name) if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row ofs << "#"; ofs << "[1]step()"; diff --git a/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp b/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp index f0377ca32..442f248f5 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp @@ -91,8 +91,7 @@ ParticleHistogram::ParticleHistogram (std::string rd_name) if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row ofs << "#"; ofs << "[1]step()"; @@ -122,7 +121,7 @@ void ParticleHistogram::ComputeDiags (int step) // Judge if the diags should be done if (!m_intervals.contains(step+1)) return; - // get WarpX class object + // get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); // get time at level 0 diff --git a/Source/Diagnostics/ReducedDiags/ParticleNumber.cpp b/Source/Diagnostics/ReducedDiags/ParticleNumber.cpp index 701698039..111ab07c0 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleNumber.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleNumber.cpp @@ -8,11 +8,13 @@ #include "ParticleNumber.H" #include "WarpX.H" +using namespace amrex::literals; + // constructor ParticleNumber::ParticleNumber (std::string rd_name) : ReducedDiags{rd_name} { - // get WarpX class object + // get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); // get MultiParticleContainer class object @@ -23,7 +25,7 @@ ParticleNumber::ParticleNumber (std::string rd_name) // resize data array to 2*(nSpecies+1) (each species + sum over all species // for both number of macroparticles and of physical particles) - m_data.resize(2*(nSpecies+1), amrex::Real(0.)); + m_data.resize(2*(nSpecies+1), 0.0_rt); // get species names (std::vector<std::string>) const auto species_names = mypc.GetSpeciesNames(); @@ -33,8 +35,7 @@ ParticleNumber::ParticleNumber (std::string rd_name) if ( m_IsNotRestart ) { // open file - std::ofstream ofs{m_path + m_rd_name + "." + m_extension, - std::ofstream::out | std::ofstream::app}; + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; // write header row ofs << "#"; ofs << "[1]step()"; @@ -95,8 +96,8 @@ void ParticleNumber::ComputeDiags (int step) const int idx_first_species_sum_weight = idx_total_sum_weight + 1; // Initialize total number of macroparticles and total weight (all species) to 0 - m_data[idx_total_macroparticles] = amrex::Real(0.); - m_data[idx_total_sum_weight] = amrex::Real(0.); + m_data[idx_total_macroparticles] = 0.0_rt; + m_data[idx_total_sum_weight] = 0.0_rt; // loop over species for (int i_s = 0; i_s < nSpecies; ++i_s) diff --git a/Source/Diagnostics/ReducedDiags/RhoMaximum.H b/Source/Diagnostics/ReducedDiags/RhoMaximum.H new file mode 100644 index 000000000..0638fae16 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/RhoMaximum.H @@ -0,0 +1,40 @@ +/* Copyright 2020 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_RHOMAXIMUM_H_ +#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_RHOMAXIMUM_H_ + +#include "ReducedDiags.H" +#include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" + +/** + * This class mainly contains a function that computes the extrema of the total charge density + * and of the charge density of each charged species. + */ +class RhoMaximum : public ReducedDiags +{ +public: + + /** constructor + * @param[in] rd_name reduced diags names + */ + RhoMaximum(std::string rd_name); + + /** This function computes the maximum and minimum values of rho (summed over all species) and + * the maximum absolute value of rho for each species. + */ + virtual void ComputeDiags(int step) override final; + +private: + /** Vector of (pointers to) functors to compute rho, per level, per species. We reuse here the + * same functors as those used for regular diagnostics. + * */ + amrex::Vector< amrex::Vector <std::unique_ptr<ComputeDiagFunctor > > > m_rho_functors; + +}; + +#endif // WARPX_DIAGNOSTICS_REDUCEDDIAGS_RHOMAXIMUM_H_ diff --git a/Source/Diagnostics/ReducedDiags/RhoMaximum.cpp b/Source/Diagnostics/ReducedDiags/RhoMaximum.cpp new file mode 100644 index 000000000..e804d4f3d --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/RhoMaximum.cpp @@ -0,0 +1,170 @@ +/* Copyright 2020 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "RhoMaximum.H" +#include "Diagnostics/ComputeDiagFunctors/RhoFunctor.H" +#include "WarpX.H" + +using namespace amrex::literals; + +// constructor +RhoMaximum::RhoMaximum (std::string rd_name) +: ReducedDiags{rd_name} +{ + + // RZ coordinate is not working +#if (defined WARPX_DIM_RZ) + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, + "RhoMaximum reduced diagnostics does not work for RZ coordinate."); +#endif + + // read number of levels + int nLevel = 0; + amrex::ParmParse pp("amr"); + pp.query("max_level", nLevel); + nLevel += 1; + m_rho_functors.resize(nLevel); + + // We do not use coarsening in the functors for this diag. + const amrex::IntVect crse_ratio = amrex::IntVect(1); + + // Initialize functors for the total charge density + for (int lev = 0; lev < nLevel; ++lev) + { + m_rho_functors[lev].push_back(std::make_unique<RhoFunctor>(lev, crse_ratio)); + } + + // get MultiParticleContainer class object + const auto & mypc = WarpX::GetInstance().GetPartContainer(); + + // get number of species (int) + const auto nSpecies = mypc.nSpecies(); + + // A vector to store the indices of charges species in mypc + amrex::Vector<int> indices_charged_species; + + // number of charged species + int n_charged_species = 0; + + for (int i = 0; i < nSpecies; ++i) + { + // Only charged species are relevant for this diag + if (mypc.GetParticleContainer(i).getCharge() != 0.0_rt) + { + indices_charged_species.push_back(i); + n_charged_species += 1; + for (int lev = 0; lev < nLevel; ++lev) + { + // Initialize functors for the charge density of each charged species + m_rho_functors[lev].push_back(std::make_unique<RhoFunctor>(lev, crse_ratio, i)); + } + } + } + + // get species names (std::vector<std::string>) + const auto species_names = mypc.GetSpeciesNames(); + + // Min and max of total rho + max of |rho| for each species + const int noutputs_per_level = 2+n_charged_species; + m_data.resize(static_cast<std::size_t>(nLevel*noutputs_per_level), 0.0_rt); + + if (amrex::ParallelDescriptor::IOProcessor()) + { + if ( m_IsNotRestart ) + { + // open file + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; + // write header row + ofs << "#"; + ofs << "[1]step()"; + ofs << m_sep; + ofs << "[2]time(s)"; + constexpr int shift_max_rho = 3; + constexpr int shift_min_rho = 4; + constexpr int shift_first_species = 5; + for (int lev = 0; lev < nLevel; ++lev) + { + ofs << m_sep; + ofs << "[" + std::to_string(shift_max_rho+noutputs_per_level*lev) + "]"; + ofs << "max_rho_lev"+std::to_string(lev)+" (C/m^3)"; + ofs << m_sep; + ofs << "[" + std::to_string(shift_min_rho+noutputs_per_level*lev) + "]"; + ofs << "min_rho_lev"+std::to_string(lev)+" (C/m^3)"; + for (int i = 0; i < n_charged_species; ++i) + { + ofs << m_sep; + ofs << "[" + std::to_string(shift_first_species+i+noutputs_per_level*lev) + "]"; + ofs << "max_" + species_names[indices_charged_species[i]] + + "_|rho|_lev"+std::to_string(lev)+" (C/m^3)"; + } + } + ofs << std::endl; + // close file + ofs.close(); + } + } + +} +// end constructor + +// function that computes maximum charge density values +void RhoMaximum::ComputeDiags (int step) +{ + // Judge if the diags should be done + if (!m_intervals.contains(step+1)) { return; } + + // get a reference to WarpX instance + auto & warpx = WarpX::GetInstance(); + + // get number of levels + const auto nLevel = warpx.finestLevel() + 1; + + const int n_charged_species = m_rho_functors[0].size() - 1; + // Min and max of total rho + max of |rho| for each species + const int noutputs_per_level = 2+n_charged_species; + + // loop over refinement levels + for (int lev = 0; lev < nLevel; ++lev) + { + // Declare a temporary MultiFAB to store the charge densities. + amrex::BoxArray ba = warpx.boxArray(lev); + amrex::DistributionMapping dmap = warpx.DistributionMap(lev); + constexpr int ncomp = 1; + constexpr int ngrow = 0; + amrex::MultiFab mf_temp(ba, dmap, ncomp, ngrow); + + // Fill temporary MultiFAB with total charge density + constexpr int idx_total_rho_functor = 0; + constexpr int idx_first_species_functor = 1; + constexpr int icomp = 0; + constexpr int i_buffer = 0; + m_rho_functors[lev][idx_total_rho_functor]->operator()(mf_temp, icomp, i_buffer); + + constexpr int idx_max_rho_data = 0; + constexpr int idx_min_rho_data = 1; + constexpr int idx_first_species_data = 2; + + // Fill output array with min and max of total rho + m_data[lev*noutputs_per_level + idx_max_rho_data] = mf_temp.max(icomp); + m_data[lev*noutputs_per_level + idx_min_rho_data] = mf_temp.min(icomp); + + // Loop over all charged species + for (int i = 0; i < n_charged_species; ++i) + { + // Fill temporary MultiFAB with the species charge density + m_rho_functors[lev][idx_first_species_functor+i]->operator()(mf_temp, icomp, i_buffer); + // Fill output array with max |rho| of species + m_data[lev*noutputs_per_level + idx_first_species_data + i] = mf_temp.norm0(); + } + } + // end loop over refinement levels + + /* m_data now contains up-to-date values for: + * [max(rho), min(rho), max(|rho_charged_species1|), max(|rho_charged_species2|), ...] */ + +} +// end void RhoMaximum::ComputeDiags |