Advanced Computing Platform for Theoretical Physics

Commit 6039e3d7 authored by Chen Xia's avatar Chen Xia
Browse files

Keep binding

galpropcc: d4f0785f408af7a795e6fbade094595bcbcabd14
parent 03136cbd
......@@ -7,9 +7,35 @@ void bind_Distribution(py::module_ &m)
.def("is2", &Dimensions::is2)
.def("is3", &Dimensions::is3)
.def("ndim", &Dimensions::ndim)
.def("npos", &Dimensions::npos);
.def("npos", &Dimensions::npos)
.def_readwrite("nx", &Dimensions::nx)
.def_readwrite("nE", &Dimensions::nE)
.def_readwrite("nGamma", &Dimensions::nGamma)
.def_readwrite("nSynch", &Dimensions::nSynch)
.def_readwrite("nlon", &Dimensions::nlat)
.def_readwrite("nRbin", &Dimensions::nRbin)
.def_readwrite("dx", &Dimensions::dx)
.def_readwrite("dy", &Dimensions::dy)
.def_readwrite("dr", &Dimensions::dr)
.def_readwrite("dz", &Dimensions::dz)
.def_readwrite("Ekin_factor", &Dimensions::Ekin_factor)
.def_readwrite("EGamma_factor", &Dimensions::EGamma_factor)
.def_readwrite("ESynch_factor", &Dimensions::ESynch_factor)
.def_readwrite("dlon", &Dimensions::dlon)
.def_readwrite("dlat", &Dimensions::dlat)
.def_readwrite("x_min", &Dimensions::x_min)
.def_readwrite("y_min", &Dimensions::y_min)
.def_readwrite("r_min", &Dimensions::r_min)
.def_readwrite("z_min", &Dimensions::z_min)
.def_readwrite("Ekin_min", &Dimensions::Ekin_min)
.def_readwrite("EGamma_min", &Dimensions::EGamma_min)
.def_readwrite("ESynch_min", &Dimensions::ESynch_min)
.def_readwrite("lon_min", &Dimensions::lon_min)
.def_readwrite("lat_min", &Dimensions::lat_min);
py::bind_vector<std::vector<Distribution>>(m, "VecDistribution");
py::implicitly_convertible<py::list, std::vector<Distribution>>();
py::implicitly_convertible<py::tuple, std::vector<Distribution>>();
py::class_<Distribution>(m, "Distribution")
.def(py::init<>())
.def("init", py::overload_cast<size_t, size_t, size_t, size_t>
......@@ -19,5 +45,23 @@ void bind_Distribution(py::module_ &m)
.def("init", py::overload_cast<const Dimensions&, size_t>
(&Distribution::init), "dims"_a, "np"_a=0)
.def("max", &Distribution::max)
.def("delete_array", &Distribution::delete_array);
.def("slice", py::overload_cast<size_t, size_t, size_t>(
&Distribution::slice))
.def("slice", py::overload_cast<size_t, size_t, size_t>(
&Distribution::slice, py::const_))
.def("slice", py::overload_cast<size_t>(&Distribution::slice))
.def("slice", py::overload_cast<size_t>(&Distribution::slice, py::const_))
// TODO .at()
.def("local_spec", &Distribution::local_spec)
.def("local_spec_1d", &Distribution::local_spec_1d)
.def("local_spec_2d", &Distribution::local_spec_2d)
.def("local_spec_3d", &Distribution::local_spec_3d)
.def("delete_array", &Distribution::delete_array)
.def_readwrite("dist", &Distribution::dist)
.def_readwrite("n_spatial_dimensions", &Distribution::n_spatial_dimensions)
.def_readwrite("n_pgrid", &Distribution::n_pgrid)
.def_readwrite("n_rgrid", &Distribution::n_rgrid)
.def_readwrite("n_xgrid", &Distribution::n_xgrid)
.def_readwrite("n_ygrid", &Distribution::n_ygrid)
.def_readwrite("n_zgrid", &Distribution::n_zgrid);
}
......@@ -4,18 +4,148 @@
void bind_Galdef(py::module_ &m)
{
py::class_<Galdef>(m, "Galdef").def(py::init<>())
.def_readwrite("version", &Galdef::version)
.def_readwrite("run_no", &Galdef::run_no)
.def_readwrite("galdef_ID", &Galdef::galdef_ID)
.def_readwrite("n_spatial_dimensions", &Galdef::n_spatial_dimensions)
.def_readwrite("z_min", &Galdef::z_min)
.def_readwrite("z_max", &Galdef::z_max)
.def_readwrite("dz", &Galdef::dz)
.def_readwrite("r_min", &Galdef::r_min)
.def_readwrite("r_max", &Galdef::r_max)
.def_readwrite("dr", &Galdef::dr)
.def_readwrite("x_min", &Galdef::x_min)
.def_readwrite("x_max", &Galdef::x_max)
.def_readwrite("dx", &Galdef::dx)
.def_readwrite("y_min", &Galdef::y_min)
.def_readwrite("y_max", &Galdef::y_max)
.def_readwrite("dy", &Galdef::dy)
.def_readwrite("Ekin_min", &Galdef::Ekin_min)
.def_readwrite("Ekin_max", &Galdef::Ekin_max)
.def_readwrite("Ekin_factor", &Galdef::Ekin_factor)
.def_readwrite("E_gamma_min", &Galdef::E_gamma_min)
.def_readwrite("E_gamma_max", &Galdef::E_gamma_max)
.def_readwrite("E_gamma_factor", &Galdef::E_gamma_factor)
.def_readwrite("nu_synch_min", &Galdef::nu_synch_min)
.def_readwrite("nu_synch_max", &Galdef::nu_synch_max)
.def_readwrite("nu_synch_factor", &Galdef::nu_synch_factor)
.def_readwrite("long_min", &Galdef::long_min)
.def_readwrite("long_max", &Galdef::long_max)
.def_readwrite("d_long", &Galdef::d_long)
.def_readwrite("lat_min", &Galdef::lat_min)
.def_readwrite("lat_max", &Galdef::lat_max)
.def_readwrite("d_lat", &Galdef::d_lat)
.def_readwrite("integration_mode", &Galdef::integration_mode)
.def_readwrite("healpix_order", &Galdef::healpix_order)
.def_readwrite("lat_substep_number", &Galdef::lat_substep_number)
.def_readwrite("LoS_step", &Galdef::LoS_step)
.def_readwrite("LoS_substep_number", &Galdef::LoS_substep_number)
.def_readwrite("start_timestep", &Galdef::start_timestep)
.def_readwrite("end_timestep", &Galdef::end_timestep)
.def_readwrite("timestep_factor", &Galdef::timestep_factor)
.def_readwrite("timestep_repeat", &Galdef::timestep_repeat)
.def_readwrite("timestep_print", &Galdef::timestep_print)
.def_readwrite("timestep_diagnostics", &Galdef::timestep_diagnostics)
.def_readwrite("control_diagnostics", &Galdef::control_diagnostics)
.def_readwrite("prop_r", &Galdef::prop_r)
.def_readwrite("prop_x", &Galdef::prop_x)
.def_readwrite("prop_y", &Galdef::prop_y)
.def_readwrite("prop_z", &Galdef::prop_z)
.def_readwrite("prop_p", &Galdef::prop_p)
.def_readwrite("use_symmetry", &Galdef::use_symmetry)
.def_readwrite("D0_xx", &Galdef::D0_xx)
.def_readwrite("D_rigid_br", &Galdef::D_rigid_br)
.def_readwrite("D_g_1", &Galdef::D_g_1)
.def_readwrite("D_g_2", &Galdef::D_g_2)
.def_readwrite("D_eta", &Galdef::D_eta)
.def_readwrite("diff_reacc", &Galdef::diff_reacc)
.def_readwrite("v_Alfven", &Galdef::v_Alfven)
.def_readwrite("damping_p0", &Galdef::damping_p0)
.def_readwrite("damping_max_path_L", &Galdef::damping_max_path_L)
.def_readwrite("damping_const_G", &Galdef::damping_const_G)
.def_readwrite("convection", &Galdef::convection)
.def_readwrite("v0_conv", &Galdef::v0_conv)
.def_readwrite("dvdz_conv", &Galdef::dvdz_conv)
.def_readwrite("He_H_ratio", &Galdef::He_H_ratio)
.def_readwrite("n_X_CO", &Galdef::n_X_CO)
.def_readwrite("n_X_CO_values", &Galdef::n_X_CO_values)
.def_readwrite("X_CO", &Galdef::X_CO)
.def_readwrite("X_CO_radius", &Galdef::X_CO_radius)
.def_readwrite("X_CO_values", &Galdef::X_CO_values)
// .def_readwrite("X_CO_parameters", &Galdef::X_CO_parameters)
.def_readwrite("propagation_X_CO", &Galdef::propagation_X_CO)
.def_readwrite("nHI_model", &Galdef::nHI_model)
.def_readwrite("nH2_model", &Galdef::nH2_model)
.def_readwrite("nHII_model", &Galdef::nHII_model)
.def_readwrite("nHI_factor", &Galdef::nHI_factor)
.def_readwrite("nHII_factor", &Galdef::nHII_factor)
.def_readwrite("nH2_factor", &Galdef::nH2_factor)
.def_readwrite("COR_filename", &Galdef::COR_filename)
.def_readwrite("HIR_filename", &Galdef::HIR_filename)
.def_readwrite("fragmentation", &Galdef::fragmentation)
.def_readwrite("momentum_losses", &Galdef::momentum_losses)
.def_readwrite("radioactive_decay", &Galdef::radioactive_decay)
.def_readwrite("K_capture)", &Galdef::K_capture)
.def_readwrite("ionization_rate)", &Galdef::ionization_rate)
.def_readwrite("source_specification)", &Galdef::source_specification)
.def_readwrite("source_normalization)", &Galdef::source_normalization)
.def_readwrite("electron_source_normalization)", &Galdef::electron_source_normalization)
.def_readwrite("use_Z", &Galdef::use_Z)
.def_readwrite("max_Z", &Galdef::max_Z)
.def_readwrite("isotopic_abundance", &Galdef::isotopic_abundance)
.def_readwrite("total_cross_section", &Galdef::total_cross_section)
.def_readwrite("cross_section_option", &Galdef::cross_section_option)
.def_readwrite("t_half_limit", &Galdef::t_half_limit)
.def_readwrite("primary_electrons", &Galdef::primary_electrons)
.def_readwrite("secondary_electrons", &Galdef::secondary_electrons)
.def_readwrite("secondary_positrons", &Galdef::secondary_positrons)
.def_readwrite("knock_on_electrons", &Galdef::knock_on_electrons)
.def_readwrite("secondary_antiDeuterium", &Galdef::secondary_antiDeuterium)
.def_readwrite("secondary_antiHelium", &Galdef::secondary_antiHelium)
.def_readwrite("secondary_protons", &Galdef::secondary_protons)
.def_readwrite("secondary_antiprotons", &Galdef::secondary_antiprotons)
.def_readwrite("tertiary_antiprotons", &Galdef::tertiary_antiprotons)
.def_readwrite("primary_DM_electrons", &Galdef::primary_DM_electrons)
.def_readwrite("primary_DM_positrons", &Galdef::primary_DM_positrons)
.def_readwrite("primary_DM_antiprotons", &Galdef::primary_DM_antiprotons)
.def_readwrite("primary_DM_antiHelium", &Galdef::primary_DM_antiHelium)
.def_readwrite("primary_DM_antiDeuterium", &Galdef::primary_DM_antiDeuterium)
.def_readwrite("tertiary_DM_antiprotons", &Galdef::tertiary_DM_antiprotons)
.def_readwrite("source_model", &Galdef::source_model)
.def_readwrite("source_model_electron", &Galdef::source_model_electron)
.def_readwrite("source_parameters", &Galdef::source_parameters)
.def_readwrite("source_parameters_electron", &Galdef::source_parameters_electron)
.def_readwrite("source_radius", &Galdef::source_radius)
.def_readwrite("source_values", &Galdef::source_values)
.def_readwrite("n_source_values", &Galdef::n_source_values)
.def_readwrite("gamma_rays", &Galdef::gamma_rays)
.def_readwrite("pi0_decay", &Galdef::pi0_decay)
.def_readwrite("IC_isotropic", &Galdef::IC_isotropic)
.def_readwrite("IC_anisotropic", &Galdef::IC_anisotropic)
.def_readwrite("bremss", &Galdef::bremss)
.def_readwrite("synchrotron", &Galdef::synchrotron)
.def_readwrite("skymap_format", &Galdef::skymap_format)
.def_readwrite("B_field_model", &Galdef::B_field_model)
.def_readwrite("B_field_name", &Galdef::B_field_name)
.def_readwrite("n_B_field_parameters", &Galdef::n_B_field_parameters)
.def_readwrite("B_field_parameters", &Galdef::B_field_parameters)
.def_readwrite("ISRF_file", &Galdef::ISRF_file)
.def_readwrite("ISRF_filetype", &Galdef::ISRF_filetype)
.def_readwrite("ISRF_healpixOrder", &Galdef::ISRF_healpixOrder)
// .def_readwrite("ISRF_factors", &Galdef::ISRF_factors)
.def_readwrite("proton_norm_Ekin", &Galdef::proton_norm_Ekin, "[MeV]")
.def_readwrite("proton_norm_flux", &Galdef::proton_norm_flux)
.def_readwrite("electron_norm_Ekin", &Galdef::electron_norm_Ekin, "[MeV]")
.def_readwrite("electron_norm_flux", &Galdef::electron_norm_flux)
.def_readwrite("verbose", &Galdef::verbose)
.def_readwrite("nuc_rigid_br", &Galdef::nuc_rigid_br)
.def_readwrite("nuc_g_1", &Galdef::nuc_g_1)
.def_readwrite("nuc_g_2", &Galdef::nuc_g_2)
.def_readwrite("electron_g_0", &Galdef::electron_g_0)
.def_readwrite("electron_g_1", &Galdef::electron_g_1)
.def_readwrite("electron_g_2", &Galdef::electron_g_2)
.def_readwrite("electron_rigid_cutoff", &Galdef::electron_rigid_cutoff)
.def_readwrite("electron_rigid_br0", &Galdef::electron_rigid_br0)
.def_readwrite("electron_rigid_br", &Galdef::electron_rigid_br)
.def_readwrite("electron_rigid_br2", &Galdef::electron_rigid_br2)
.def_readwrite("max_Z", &Galdef::max_Z);
.def_readwrite("electron_rigid_br2", &Galdef::electron_rigid_br2);
}
......@@ -48,18 +48,63 @@ void bind_Particle(py::module_ &m)
.def(py::init<std::shared_ptr<const ParticleInfo>,
std::shared_ptr<const Galaxy>,
const std::string &>())
.def("info", &Particle::info)
.def("galaxy", &Particle::galaxy)
.def("dims", &Particle::dims)
.def("Z", &Particle::Z)
.def("A", &Particle::A)
.def("A_pos", &Particle::A_pos)
.def("is_lepton", &Particle::is_lepton)
.def("is_isotope", &Particle::is_isotope)
.def("depends_on", &Particle::depends_on)
.def("fill_transport_arrays", &Particle::fill_transport_arrays)
.def("fill_source_functions", &Particle::fill_source_functions)
.def("fill_diffusion_coeffs", &Particle::fill_diffusion_coeffs)
.def("fill_energy_loss_rate", &Particle::fill_energy_loss_rate)
.def("fill_destruction_rate", &Particle::fill_destruction_rate)
.def("propagate", &Particle::propagate)
.def("abundance", &Particle::abundance)
.def("local_flux", py::overload_cast<double, double>(
&Particle::local_flux, py::const_),
"2-D local flux", "r"_a=Rsun, "z"_a=0)
.def("local_flux", py::overload_cast<double, double, double>(
&Particle::local_flux, py::const_),
"3-D local flux", "x"_a, "y"_a, "z"_a)
.def("r", py::overload_cast<>(&Particle::r, py::const_))
.def("r", py::overload_cast<size_t>(&Particle::r, py::const_))
.def("x", py::overload_cast<>(&Particle::x, py::const_))
.def("x", py::overload_cast<size_t>(&Particle::x, py::const_))
.def("y", py::overload_cast<>(&Particle::y, py::const_))
.def("y", py::overload_cast<size_t>(&Particle::y, py::const_))
.def("z", py::overload_cast<>(&Particle::z, py::const_))
.def("z", py::overload_cast<size_t>(&Particle::z, py::const_))
.def("Ekin", py::overload_cast<>(&Particle::Ekin, py::const_))
.def("Ekin", py::overload_cast<size_t>(&Particle::Ekin, py::const_))
.def("Etot", py::overload_cast<>(&Particle::Etot, py::const_))
.def("Etot", py::overload_cast<size_t>(&Particle::Etot, py::const_))
.def("p", py::overload_cast<>(&Particle::p, py::const_))
.def("p", py::overload_cast<size_t>(&Particle::p, py::const_))
.def("rigidity", py::overload_cast<>(&Particle::rigidity, py::const_))
.def("rigidity", py::overload_cast<size_t>(&Particle::rigidity, py::const_))
.def("gamma", py::overload_cast<>(&Particle::gamma, py::const_))
.def("gamma", py::overload_cast<size_t>(&Particle::gamma, py::const_))
.def("beta", py::overload_cast<>(&Particle::beta, py::const_))
.def("beta", py::overload_cast<size_t>(&Particle::beta, py::const_))
.def("Egamma", py::overload_cast<>(&Particle::Egamma, py::const_))
.def("Egamma", py::overload_cast<size_t>(&Particle::Egamma, py::const_))
.def("Eradio", py::overload_cast<>(&Particle::Eradio, py::const_))
.def("Eradio", py::overload_cast<size_t>(&Particle::Eradio, py::const_))
.def("Ekin_GeV", &Particle::Ekin_GeV)
.def("t_half", &Particle::t_half)
.def_readwrite("name", &Particle::name)
.def_readwrite("K_electron", &Particle::K_electron);
.def_readwrite("K_electron", &Particle::K_electron)
.def_readwrite("cr_density", &Particle::cr_density)
.def_readwrite("source_function", &Particle::source_function)
.def_readwrite("destruct_rate", &Particle::destruct_rate)
.def_readwrite("Dxx", &Particle::Dxx)
.def_readwrite("Dpp", &Particle::Dpp)
.def_readwrite("dpdt", &Particle::dpdt)
.def_readwrite("normalization_factor", &Particle::normalization_factor)
.def_readwrite("primary_source_function", &Particle::primary_source_function)
.def_readwrite("secondary_source_function", &Particle::secondary_source_function);
}
......@@ -4,10 +4,20 @@ void bind_Structures(py::module_ &m)
{
py::class_<ParticleNumber>(m, "ParticleNumber")
.def(py::init<int, int>(), "z"_a, "a"_a=ParticleNumber::AnyA)
.def(py::init([](const py::tuple& args) {
std::size_t arg_len = py::len(args);
if (arg_len < 1 || arg_len > 2) {
throw py::value_error("tuple length incorrect");
}
int Z = args[0].cast<int>();
if (arg_len == 1) { return ParticleNumber(Z); }
return ParticleNumber(Z, args[1].cast<int>());
}))
.def_readwrite("Z", &ParticleNumber::Z)
.def_readwrite("A", &ParticleNumber::A)
.def_readonly_static("AnyA", &ParticleNumber::AnyA)
.def(py::self == py::self)
.def(py::self != py::self)
.def(py::self < py::self);
py::implicitly_convertible<py::tuple, ParticleNumber>();
}
......@@ -3,6 +3,7 @@
namespace py = pybind11;
void bind_ErrorLogger(py::module_ &);
void bind_Distribution(py::module_ &m);
void bind_SourceSpectrum(py::module_ &);
void bind_SourceDistribution(py::module_ &);
void bind_Structures(py::module_ &);
......@@ -19,6 +20,7 @@ void bind_Galprop(py::module_ &);
PYBIND11_MODULE(_bind_galpropcc, m) {
bind_ErrorLogger(m);
bind_Distribution(m);
bind_SourceSpectrum(m);
bind_SourceDistribution(m);
bind_Structures(m);
......
from ._bind_galpropcc import *
from . import tools
from . import populators
""" CR sources """
from ._bind_galpropcc import (
Galprop, GCR, SrcSpecSBPL, EleSrcDistParameterized, MapPtrParticle,
PrimaryParticle)
class PopPrimaryElecSBPL(GCR.Populator):
def __init__(self, *, segments=None, cutoff_rig=None):
GCR.Populator.__init__(self)
self._name = "primary_electrons"
self._spectral = SrcSpecSBPL()
self._spectral.abundance = 1
self._spatial = EleSrcDistParameterized()
if segments is not None:
self.set_segments(segments)
if cutoff_rig is not None:
self.set_cutoff_rig(cutoff_rig)
def set_segments(self, segments: list[tuple]) -> None:
self._spectral.segments = segments
def set_cutoff_rig(self, rig: float) -> None:
self._spectral.cutoff_rig = rig
def set_norm_ekn_GeV(self, v_GeV: float) -> None:
self._spectral.set_norm_ekn_GeV(v_GeV)
def set_norm_rig_GV(self, v_GV: float) -> None:
self._spectral.set_norm_rig_GV(v_GV)
def populate(self, g: Galprop) -> MapPtrParticle:
if not g.galdef.primary_electrons:
return MapPtrParticle()
self._spatial.pG = g
iele = g.particle_lib.get(-1, 0)
pele = PrimaryParticle(iele, g.galaxy, self._name)
pele.spec_ptr = self._spectral
pele.dist_ptr = self._spatial
particle_map = MapPtrParticle()
particle_map[self._name] = pele
return particle_map
""" Usefull tools """
import os
__all__ = ['dump_gcr']
def dump_gcr(g, outfile, particles="all", mode='w-', **kwargs):
"""Dump full cosmic ray distribution.
Parameters
----------
g : galpropy.Galprop
The Galprop instance.
outfile : str
Output file name including the path. The directory will be created if
does not exist.
particles : list[str] or set[str], or "all", optional
Specify the Particles to dump. (default "all")
mode : {'w-', 'x', 'w', 'a'}, optional
w- or x Create file, fail if exists (default)
w Create file, truncate if exists
a Read/write if exists, create otherwise
kwargs : dict, optional
Extra arguments to `h5py.File`.
"""
try:
import h5py
import numpy as np
except ImportError:
print("To store GCR please install h5py and numpy, nothing dumped")
return
outfile = os.path.abspath(outfile)
os.makedirs(os.path.dirname(outfile), exist_ok=True)
if particles == "all":
particles = {p.name for p in g.gcr}
with h5py.File(outfile, mode, **kwargs) as f:
dset = f.create_dataset('r', data=g.gcr[0].r())
dset.attrs['unit'] = 'kpc'
dset = f.create_dataset('z', data=g.gcr[0].z())
dset.attrs['unit'] = 'kpc'
dset = f.create_dataset('Ekin', data=g.GetEkin())
dset.attrs['unit'] = 'GeV/n'
for gcr in g.gcr:
if gcr.name in particles:
group = f.create_group(gcr.name)
crdist = gcr.cr_density
dist = gcr.A_pos() * 1e3 * np.array(crdist.dist).reshape(
(crdist.n_rgrid, crdist.n_zgrid, crdist.n_pgrid))
group.create_dataset('flux', data=dist)
group.attrs['Z'] = gcr.Z()
group.attrs['A'] = gcr.A()
group.attrs['unit'] = r'cm^-2 s^-1 sr^-1 GeV^-1'
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment