Skip to content

Commit aa56b07

Browse files
StaticObserverStaticObserver
committed
working version
1 parent ad88e33 commit aa56b07

4 files changed

Lines changed: 158 additions & 112 deletions

File tree

‎src/kernels/QED_process.hpp‎

Lines changed: 146 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,152 @@
1818
namespace kernel::QED{
1919
using namespace ntt;
2020

21+
22+
23+
class cdfTable {
24+
private:
25+
array_t<real_t*> cdf;
26+
array_t<real_t*> inv_cdf;
27+
28+
real_t x_min_lg;
29+
real_t x_max_lg;
30+
real_t y_min_lg;
31+
real_t y_max_lg;
32+
real_t dx;
33+
real_t dy;
34+
35+
size_t Nx;
36+
size_t Ny;
37+
38+
// 改进的文件读取方法
39+
auto read_from_file(const std::string& filename) -> std::pair<std::vector<real_t>, std::vector<real_t>> {
40+
std::fstream file;
41+
file.open(filename, std::ios::in);
42+
if (!file.is_open()) {
43+
raise::Fatal("Could not open file.", HERE);
44+
}
45+
46+
std::vector<real_t> x_data;
47+
std::vector<real_t> y_data;
48+
49+
while (!file.eof()) {
50+
real_t x, y;
51+
if (file >> x >> y) { // 确保成功读取
52+
if(x <= ZERO || y <= ZERO){
53+
raise::Fatal("Invalid data value (must be positive).", HERE);
54+
}
55+
x_data.push_back(x);
56+
y_data.push_back(y);
57+
}
58+
}
59+
file.close();
60+
61+
if (x_data.empty()) {
62+
raise::Fatal("Empty data file.", HERE);
63+
}
64+
65+
return {x_data, y_data};
66+
}
67+
68+
public:
69+
cdfTable(const std::string& cdf_filename, const std::string& inverse_cdf_filename) {
70+
// 从文件读取数据，使用std::vector暂存
71+
auto [x_cdf_data, y_cdf_data] = read_from_file(cdf_filename);
72+
auto [x_inv_cdf_data, y_inv_cdf_data] = read_from_file(inverse_cdf_filename);
73+
74+
// 设置尺寸
75+
Nx = x_cdf_data.size();
76+
Ny = x_inv_cdf_data.size();
77+
78+
// 设置范围参数
79+
x_min_lg = math::log10(x_cdf_data[0]);
80+
x_max_lg = math::log10(x_cdf_data[Nx - 1]);
81+
dx = (x_max_lg - x_min_lg) / (Nx - 1);
82+
y_min_lg = math::log10(x_inv_cdf_data[0]);
83+
y_max_lg = math::log10(x_inv_cdf_data[Ny - 1]);
84+
dy = (y_max_lg - y_min_lg) / (Ny - 1);
85+
86+
// 分配设备内存
87+
cdf = array_t<real_t*>("cdf", Nx);
88+
inv_cdf = array_t<real_t*>("inv_cdf", Ny);
89+
90+
// 创建主机镜像
91+
auto cdf_h = Kokkos::create_mirror_view(cdf);
92+
auto inv_cdf_h = Kokkos::create_mirror_view(inv_cdf);
93+
94+
// 填充数据
95+
for (size_t i = 0; i < Nx; ++i) {
96+
cdf_h(i) = y_cdf_data[i];
97+
}
98+
99+
for (size_t i = 0; i < Ny; ++i) {
100+
inv_cdf_h(i) = y_inv_cdf_data[i];
101+
}
102+
103+
// 复制到设备内存
104+
Kokkos::deep_copy(cdf, cdf_h);
105+
Kokkos::deep_copy(inv_cdf, inv_cdf_h);
106+
}
107+
108+
// 拷贝构造函数保持不变，因为是浅拷贝
109+
Inline cdfTable(const cdfTable &rhs)
110+
: cdf(rhs.cdf)
111+
, inv_cdf(rhs.inv_cdf)
112+
, x_min_lg(rhs.x_min_lg)
113+
, x_max_lg(rhs.x_max_lg)
114+
, y_min_lg(rhs.y_min_lg)
115+
, y_max_lg(rhs.y_max_lg)
116+
, dx(rhs.dx)
117+
, dy(rhs.dy)
118+
, Nx(rhs.Nx)
119+
, Ny(rhs.Ny) {}
120+
121+
~cdfTable() = default;
122+
123+
// 保持原始的CDF超出范围处理逻辑
124+
Inline auto CDF(const real_t x) const -> real_t {
125+
if (x <= ZERO) {
126+
raise::KernelError(HERE, "Invalid argument for CDF.");
127+
}
128+
const int idx = static_cast<int>((math::log10(x) - x_min_lg) / dx);
129+
if (idx < 0) {
130+
return ONE + 0.346 * x - math::pow(x, ONE / THREE) * (1.232 + 0.033 * SQR(x));
131+
}
132+
auto i = static_cast<size_t>(idx);
133+
if (i >= Nx - 1) {
134+
return cdf(Nx - 1);
135+
}
136+
137+
const real_t t = (x - math::pow(10.0, i * dx + x_min_lg))
138+
/ (math::pow(10.0, (i + 1) * dx + x_min_lg) - math::pow(10.0, i * dx + x_min_lg));
139+
140+
// 插值公式
141+
return cdf(i) * (ONE - t) + cdf(i + 1) * t;
142+
}
143+
144+
// 保持原始的Inverse_CDF超出范围处理逻辑
145+
Inline auto Inverse_CDF(const real_t y) const -> real_t {
146+
if (y <= ZERO || y >= ONE) {
147+
raise::KernelError(HERE, "Invalid argument for Inverse_CDF.");
148+
}
149+
const int idx = static_cast<int>((math::log10(y) - y_min_lg) / dy);
150+
if (idx < 0) {
151+
return math::pow(10.0, x_max_lg);
152+
}
153+
auto i = static_cast<size_t>(idx);
154+
if (i >= Ny - 1) {
155+
return CUBE((ONE - y) / 1.232);
156+
}
157+
158+
const real_t t = (y - math::pow(10.0, i * dy + y_min_lg))
159+
/ (math::pow(10.0, (i + 1) * dy + y_min_lg) - math::pow(10.0, i * dy + y_min_lg));
160+
161+
// 原始插值公式
162+
return inv_cdf(i) * (ONE - t) + inv_cdf(i + 1) * t;
163+
}
164+
}; // class cdfTable
165+
166+
21167
template <Dimension D, Coord::type C>
22168
class CurvatureEmission_kernel{
23169
static_assert(D == Dim::_1D, "Curvature emission is only implemented in 1D");
@@ -378,110 +524,6 @@ namespace kernel::QED{
378524

379525
}; // class InjectPairs_kernel
380526

381-
class cdfTable{
382-
array_t<real_t*> cdf;
383-
array_t<real_t*> inv_cdf;
384-
385-
real_t x_min;
386-
real_t x_max;
387-
real_t y_min;
388-
real_t y_max;
389-
real_t dx;
390-
real_t dy;
391-
392-
size_t Nx;
393-
size_t Ny;
394-
395-
396-
397-
auto read_from_file(const std::string& filename) -> array_t<real_t**>{
398-
std::fstream file;
399-
file.open(filename, std::ios::in);
400-
if (!file.is_open()){
401-
raise::Fatal("Could not open file.", HERE);
402-
}
403-
std::vector<real_t> x_data;
404-
std::vector<real_t> y_data;
405-
while (!file.eof()){
406-
real_t x, y;
407-
file >> x >> y;
408-
x_data.push_back(x);
409-
y_data.push_back(y);
410-
}
411-
file.close();
412-
auto N = x_data.size();
413-
auto data = array_t<real_t**>("data", N);
414-
auto data_h = Kokkos::create_mirror_view(data);
415-
for (size_t i = 0; i < N; ++i){
416-
data_h(i, 0) = x_data[i];
417-
data_h(i, 1) = y_data[i];
418-
}
419-
Kokkos::deep_copy(data, data_h);
420-
return data;
421-
}
422-
423-
public:
424-
cdfTable(const std::string& cdf_filename, const std::string& inverse_cdf_filename)
425-
{
426-
auto cdf_data = read_from_file(cdf_filename);
427-
auto inv_cdf_data = read_from_file(inverse_cdf_filename);
428-
429-
Nx = cdf_data.extent(0);
430-
Ny = inv_cdf_data.extent(0);
431-
432-
x_min = cdf_data(0, 0);
433-
x_max = cdf_data(Nx - 1, 0);
434-
dx = (math::log10(x_max) - math::log10(x_min)) / (Nx - 1);
435-
y_min = inv_cdf_data(0, 0);
436-
y_max = inv_cdf_data(Ny - 1, 0);
437-
dy = (y_max - y_min) / (Ny - 1);
438-
439-
cdf = array_t<real_t*>("cdf", Nx);
440-
inv_cdf = array_t<real_t*>("inv_cdf", Ny);
441-
442-
auto cdf_subview = Kokkos::subview(cdf_data, Kokkos::ALL, 1);
443-
auto inv_cdf_subview = Kokkos::subview(inv_cdf_data, Kokkos::ALL, 1);
444-
Kokkos::deep_copy(cdf, cdf_subview);
445-
Kokkos::deep_copy(inv_cdf, inv_cdf_subview);
446-
}
447-
448-
Inline cdfTable(const cdfTable &rhs)
449-
: cdf(rhs.cdf)
450-
, inv_cdf(rhs.inv_cdf)
451-
, x_min(rhs.x_min)
452-
, x_max(rhs.x_max)
453-
, y_min(rhs.y_min)
454-
, y_max(rhs.y_max)
455-
, dx(rhs.dx)
456-
, dy(rhs.dy)
457-
, Nx(rhs.Nx)
458-
, Ny(rhs.Ny) {}
459-
460-
~cdfTable() = default;
461-
462-
Inline auto CDF(const real_t x) const -> real_t{
463-
if (x < x_min){
464-
return ONE + 0.346 * x - math::pow(x, 1/3) * (1.232 + 0.033 * SQR(x));
465-
}
466-
if (x > x_max){
467-
return cdf(Nx - 1) * math::exp(x - x_max);
468-
}
469-
const size_t i = static_cast<size_t>((math::log10(x) - math::log10(x_min)) / dx);
470-
return cdf(i) + (x - x_min - i * dx) / dx * (cdf(i + 1) - cdf(i));
471-
}
472-
473-
Inline auto Inverse_CDF(const real_t y) const -> real_t{
474-
if (y < y_min){
475-
return x_max + math::log(y_min / y);
476-
}
477-
if (y > y_max){
478-
return CUBE((ONE - y) / 1.232);
479-
}
480-
const size_t i = static_cast<size_t>((y - y_min) / dy);
481-
return inverse_cdf(i) + (y - y_min - i * dy) / dy * (inverse_cdf(i + 1) - inverse_cdf(i));
482-
}
483-
484-
}
485527
} // namespace QED
486528

487529

‎src/kernels/tests/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,4 +33,4 @@ gen_test(gca_pusher)
3333
gen_test(prtl_bc)
3434
gen_test(flds_bc)
3535
gen_test(QED_process)
36-
gen+test(cdf)
36+
gen_test(cdf)

‎src/kernels/tests/QED_process.cpp‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -64,14 +64,17 @@ auto main(int argc, char* argv[]) -> int {
6464

6565
auto start = std::chrono::high_resolution_clock::now();
6666
array_t<size_t*> N_phs("N_phs", electron.npart());
67-
67+
68+
cdfTable cdf("cdf_table.txt", "inverse_cdf_table.txt" );
6869
// std::cout << "Begin curvature emission number." << std::endl;
69-
Curvature_Emission_Number<Dim::_1D, Coord::Cart> curvature_number(electron,
70+
Curvature_Emission_Number<Dim::_1D, Coord::Cart> curvature_number(electron,
7071
e_min,
7172
coeff * 1e9,
73+
1e5,
7274
gamma_emit,
7375
rho,
74-
N_phs);
76+
N_phs,
77+
cdf);
7578

7679

7780
Kokkos::parallel_for("CurvatureEmissionNumber", electron.rangeActiveParticles(), curvature_number);
@@ -116,7 +119,8 @@ auto main(int argc, char* argv[]) -> int {
116119
100,
117120
N_phs,
118121
offsets,
119-
random_pool);
122+
random_pool,
123+
cdf);
120124

121125
Kokkos::parallel_for("CurvatureEmission", electron.rangeActiveParticles(), curvature_emission);
122126

‎src/kernels/tests/cdf.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,9 @@ auto main(int argc, char* argv[]) -> int {
3030
array_t<real_t*> y_inv_cdf("y", N);
3131

3232
const real_t x_min { 1e-8 };
33-
const real_t x_max { 1e2 };
34-
const real_t y_min { 1e-8 };
35-
const real_t y_max { 1.0 };
33+
const real_t x_max { 11.13 };
34+
const real_t y_min { 1e-6 };
35+
const real_t y_max { 0.99 };
3636

3737
const real_t dx = (math::log10(x_max) - math::log10(x_min)) / (N - 1);
3838
const real_t dy = (math::log10(y_max) - math::log10(y_min)) / (N - 1);

0 commit comments

Comments
 (0)