Skip to content

spect3d added #193

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
a-sullivan wants to merge 1 commit into
base: 1.4.0rc
Choose a base branch
from
Draft

spect3d added #193

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 205 additions & 1 deletion src/framework/domain/metadomain_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -345,8 +345,13 @@ namespace ntt {
g_writer.shouldWrite("spectra",
finished_step,
finished_time);
const auto write_spectra3D = params.template get<bool>(
"output.spectra3D.enable") and
g_writer.shouldWrite("spectra3D",
finished_step,
finished_time);
const auto extension = params.template get<std::string>("output.format");
if (not(write_fields or write_particles or write_spectra) and
if (not(write_fields or write_particles or write_spectra or write_spectra3D) and
extension != "disabled") {
return false;
}
Expand Down Expand Up @@ -788,6 +793,205 @@ namespace ntt {
g_writer.endWriting(WriteMode::Spectra);
} // end shouldWrite("spectra", step, time)

if (write_spectra3D) {
g_writer.beginWriting(WriteMode::Spectra3D, current_step, current_time);
const auto log_bins = params.template get<bool>(
"output.spectra3D.log_bins");
const auto n_bins = params.template get<std::size_t>(
"output.spectra3D.n_bins");
const auto& metric = local_domain->mesh.metric;
// extract the number of bins globally in each direction
const auto nx1_bins = params.template get<std::size_t>("output.spectra3D.nx1");
const auto nx2_bins = params.template get<std::size_t>("output.spectra3D.nx2");
const auto nx3_bins = params.template get<std::size_t>("output.spectra3D.nx3");

// select the min and max energy for the spectra
auto e_min = params.template get<real_t>("output.spectra3D.e_min");
auto e_max = params.template get<real_t>("output.spectra3D.e_max");



auto x1_min = mesh().extent(in::x1).first; // extent is in physical units, not code
decltype(x1_min) x2_min = 0;
decltype(x1_min) x3_min = 0;

auto x1_max = mesh().extent(in::x1).second; // pairs in c++ are addressed by first and secocnd
decltype(x1_max) x2_max = 0;
decltype(x1_max) x3_max = 0;



real_t dx1_bin = (x1_max - x1_min )/nx1_bins;
real_t dx2_bin = 1.;
real_t dx3_bin = 1.;



if constexpr (M::PrtlDim == Dim::_2D or M::PrtlDim == Dim::_3D){ // only pick x2 if simulation is in 2D

x2_min = mesh().extent(in::x2).first;
x2_max = mesh().extent(in::x2).second;



dx2_bin = (static_cast<real_t>(x2_max) - static_cast<real_t>(x2_min) )/static_cast<real_t>(nx2_bins);

}

if constexpr (M::PrtlDim == Dim::_3D){

x3_min = mesh().extent(in::x3).first;
x3_max = mesh().extent(in::x3).second;


dx3_bin = (static_cast<real_t>(x3_max) - static_cast<real_t>(x3_min) )/static_cast<real_t>(nx3_bins);


}



//auto dx_slice = mesh().dx1;


//auto dxslice = (x1_max - x1_min) / x_bins

if (log_bins) {
e_min = math::log10(e_min);
e_max = math::log10(e_max);
}
array_t<real_t*> energy { "energy", n_bins + 1 };
// generating the energy bins
Kokkos::parallel_for(
"GenerateEnergyBins",
n_bins + 1,
Lambda(index_t e) {
if (log_bins) {
energy(e) = math::pow(10.0, e_min + (e_max - e_min) * e / n_bins);
} else {
energy(e) = e_min + (e_max - e_min) * e / n_bins;
}
});

for (const auto& spec : g_writer.spectraWriters()) {
auto& species = local_domain->species[spec.species() - 1];
array_t<real_t****> dn3d { "dn3d", nx1_bins, nx2_bins, nx3_bins, n_bins };
auto dn3d_scatter = Kokkos::Experimental::create_scatter_view(dn3d);
auto ux1 = species.ux1;
auto ux2 = species.ux2;
auto ux3 = species.ux3;
auto i1 = species.i1;
auto dx1 = species.dx1;
//adeep_copy;
decltype(i1) i2;
decltype(i1) i3;
decltype(dx1) dx2;
decltype(dx1) dx3;
if constexpr (M::PrtlDim == Dim::_2D or M::PrtlDim == Dim::_3D){ // only pick x2 if simulation is in 2D
i2 = species.i2;
dx2 = species.dx2;
}
if constexpr (M::PrtlDim == Dim::_3D){
i3 = species.i3;
dx3 = species.dx3;
}
auto weight = species.weight;
auto tag = species.tag;
const auto is_massive = species.mass() > 0.0f;
Kokkos::parallel_for(
"ComputeSpectra",
species.rangeActiveParticles(),
Lambda(index_t p) {
if (tag(p) != ParticleTag::alive) {
return;
}

coord_t<D> x_Cd {ZERO};
if (D == Dim::_1D or D == Dim::_2D or D == Dim::_3D) {
x_Cd[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
}
if (D == Dim::_2D or D == Dim::_3D) {
x_Cd[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
}
if (D == Dim::_3D) {
x_Cd[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
}
coord_t<D> x_Ph { ZERO };
metric.template convert<Crd::Cd, Crd::Ph>(x_Cd, x_Ph);

real_t en;
if (is_massive) {
en = U2GAMMA(ux1(p), ux2(p), ux3(p)) - ONE;
} else {
en = NORM(ux1(p), ux2(p), ux3(p));
}
if (log_bins) {
en = math::log10(en);
}
std::size_t e_ind = 0;
if (en <= e_min) {
e_ind = 0;
} else if (en >= e_max) {
e_ind = n_bins-1;
} else {
e_ind = static_cast<std::size_t>(
static_cast<real_t>(n_bins) * (en - e_min) / (e_max - e_min));
}
std::cout << "[x1_min, x1_max] =[" << x1_min<< ", " << x1_max << "]" << std::endl;
std::size_t x1_ind = 0;
if (x_Ph[0]<= x1_min) {
x1_ind = 0;
} else if (x_Ph[0] >= x1_max) {
x1_ind = nx1_bins-1;
} else {
x1_ind = static_cast<std::size_t>(
static_cast<real_t>(nx1_bins) * (x_Ph[0] - x1_min) / (x1_max - x1_min));
}


std::size_t x2_ind = 0;
if (M::PrtlDim == Dim::_2D or M::PrtlDim == Dim::_3D){ // only pick x2 if simulation is in 2D
if (x_Ph[1] <= x2_min) {
x2_ind = 0;
} else if (x_Ph[1] >= x2_max) {
x2_ind = nx2_bins-1;
} else {
x2_ind = static_cast<std::size_t>(
static_cast<real_t>(nx2_bins) * (x_Ph[1] - x2_min) / (x2_max - x2_min));
}
//real_t x2_center = x2_min + (x2_ind+0.5)*dx2_bin;
//owns_x2 = (x2_center >= x2_min_local) && (x2_center <= x2_max_local);
}
std::size_t x3_ind = 0;
if (M::PrtlDim == Dim::_3D){ // only pick x3 if simulation is in 3D

if (x_Ph[2] <= x3_min) {
x3_ind = 0;
} else if (x_Ph[2] >= x3_max) {
x3_ind = nx3_bins-1;
} else {
x3_ind = static_cast<std::size_t>(
static_cast<real_t>(nx3_bins) * (x_Ph[2] - x3_min) / (x3_max - x3_min));
}


}

// now I want to save the nx_bins contained in each rank
// can save array of x1_inds which are saved?
// maybe I can ask what local x_min and x_max are for this rank, pass that to writeSpectrum3D


auto dn3d_acc = dn3d_scatter.access();
dn3d_acc(x1_ind, x2_ind, x3_ind, e_ind) += weight(p);
});
Kokkos::Experimental::contribute(dn3d, dn3d_scatter);
g_writer.writeSpectrum3D(dn3d, spec.name());
}
g_writer.writeSpectrumBins(energy, "sEbn");
g_writer.endWriting(WriteMode::Spectra3D);
} // end shouldWrite("spectra", step, time)

return true;
}

Expand Down
44 changes: 43 additions & 1 deletion src/framework/parameters/output.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ namespace ntt {
"separate_files=false is deprecated",
HERE);

for (const auto& category : { "fields", "particles", "spectra", "stats" }) {
for (const auto& category : { "fields", "particles", "spectra", "spectra3D", "stats" }) {
const auto q_int = toml::find_or<timestep_t>(toml_data,
"output",
category,
Expand Down Expand Up @@ -143,6 +143,40 @@ namespace ntt {
"spectra",
"n_bins",
defaults::output::spec_nbins);
/* Spectra3D ------------------------------------------------------------ */
spectra3D_e_min = toml::find_or(toml_data, "output", "spectra3D", "e_min", defaults::output::spec3d_emin);

spectra3D_e_max = toml::find_or(toml_data, "output", "spectra3D", "e_max", defaults::output::spec3d_emax);

spectra3D_log_bins = toml::find_or(toml_data,
"output",
"spectra3D",
"log_bins",
defaults::output::spec3d_log);

spectra3D_n_bins = toml::find_or(toml_data,
"output",
"spectra3D",
"n_bins",
defaults::output::spec3d_nbins);

spectra3D_nx1 = toml::find_or(toml_data,
"output",
"spectra3D",
"nx1",
defaults::output::spec3d_nx1);

spectra3D_nx2 = toml::find_or(toml_data,
"output",
"spectra3D",
"nx2",
defaults::output::spec3d_nx2);
spectra3D_nx3 =
toml::find_or(toml_data,
"output",
"spectra3D",
"nx3",
defaults::output::spec3d_nx3);

/* Stats ---------------------------------------------------------------- */
stats_quantities = toml::find_or(toml_data,
Expand Down Expand Up @@ -193,6 +227,14 @@ namespace ntt {
params->set("output.spectra.log_bins", spectra_log_bins);
params->set("output.spectra.n_bins", spectra_n_bins);

params->set("output.spectra3D.e_min", spectra3D_e_min);
params->set("output.spectra3D.e_max", spectra3D_e_max);
params->set("output.spectra3D.log_bins", spectra3D_log_bins);
params->set("output.spectra3D.n_bins", spectra3D_n_bins);
params->set("output.spectra3D.nx1", spectra3D_nx1);
params->set("output.spectra3D.nx2", spectra3D_nx2);
params->set("output.spectra3D.nx3", spectra3D_nx3);

params->set("output.stats.quantities", stats_quantities);
params->set("output.stats.custom", stats_custom_quantities);

Expand Down
8 changes: 8 additions & 0 deletions src/framework/parameters/output.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,14 @@ namespace ntt {
bool spectra_log_bins;
std::size_t spectra_n_bins;

real_t spectra3D_e_min;
real_t spectra3D_e_max;
bool spectra3D_log_bins;
std::size_t spectra3D_n_bins;
std::size_t spectra3D_nx1;
std::size_t spectra3D_nx2;
std::size_t spectra3D_nx3;

std::vector<std::string> stats_quantities;
std::vector<std::string> stats_custom_quantities;

Expand Down
8 changes: 8 additions & 0 deletions src/global/defaults.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,14 @@ namespace ntt::defaults {
const real_t spec_emax = 1e3;
const bool spec_log = true;
const std::size_t spec_nbins = 200;
const real_t spec3d_emin = 1e-3;
const real_t spec3d_emax = 1e3;
const bool spec3d_log = true;
const std::size_t spec3d_nbins = 200;
const std::size_t spec3d_nx1 = 1;
const std::size_t spec3d_nx2 = 1;
const std::size_t spec3d_nx3 = 1;

const std::vector<std::string> stats_quantities = { "B^2",
"E^2",
"ExB",
Expand Down
1 change: 1 addition & 0 deletions src/global/global.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,7 @@ namespace WriteMode {
Particles = 1 << 1,
Spectra = 1 << 2,
Stats = 1 << 3,
Spectra3D = 1 << 4,
};
} // namespace WriteMode

Expand Down
Loading