GSI_Toolbox 1.0.0
A toolbox for Gas-Surface Interaction simulations
Loading...
Searching...
No Matches
PG_Kernel.h
1
9template <typename T>
10struct trajectory {
11 public:
12 std::vector<Matrix<T>> velocities;
13 std::vector<Matrix<T>> positions;
14};
15
16
25template <typename T>
26class PG_Kernel {
27
28 private:
29
30 Surface<T> surface; // Surface object for storing surface properties
31 Gas<T> gas; // Gas object for storing gas properties
32
33 Matrix<T> roots; // Matrix of roots for Hermite polynomials
34 Matrix<T> weights; // Matrix of weights for Hermite polynomials
35 Matrix<T> roots_2d_r; // Matrix of roots for 2D Hermite polynomials (row)
36 Matrix<T> roots_2d_c; // Matrix of roots for 2D Hermite polynomials (column)
37 Matrix<T> weights_2d; // Matrix of weights for 2D Hermite polynomials
38 Matrix<T> gamma; // Matrix of gamma values for PG_Kernel
39 Matrix<T> gamma_dot; // Matrix of gamma_dot values for PG_Kernel
40 Matrix<T> gamma_2d; // Matrix of gamma values for 2D PG_Kernel
41 Matrix<T> gamma_dot_2d; // Matrix of gamma_dot values for 2D PG_Kernel
42 Matrix<T> p_gamma; // Matrix of p_gamma values for PG_Kernel
43 Matrix<T> p_gamma_dot; // Matrix of p_gamma_dot values for PG_Kernel
44 Matrix<T> p_gamma_2d; // Matrix of p_gamma values for 2D PG_Kernel
45 Matrix<T> p_gamma_dot_2d; // Matrix of p_gamma_dot values for 2D PG_Kernel
46 Matrix<T> mu; // Matrix of mu values for PG_Kernel
47 Matrix<T> sigma; // Matrix of sigma values for PG_Kernel
48 Matrix<T> mu_2d; // Matrix of mu values for 2D PG_Kernel
49 Matrix<T> sigma_2d; // Matrix of sigma values for 2D PG_Kernel
50 Matrix<T> mu_gamma; // Matrix of mu_gamma values for PG_Kernel
51 Matrix<T> sigma_gamma; // Matrix of sigma_gamma values for PG_Kernel
52 Matrix<T> mu_gamma_2d; // Matrix of mu_gamma values for 2D PG_Kernel
53 Matrix<T> sigma_gamma_2d; // Matrix of sigma_gamma values for 2D PG_Kernel
54
55 Matrix<T> incident_velocity; // Matrix of incident velocities
56
57 unsigned int int_order; // Order of integration for Hermite polynomials
58 unsigned long num_particles; // Number of particles
59
60 std::string sim_name; // Name of the simulation
61
62 T margin, ac_length, mu_gas, T_gas; // Parameters for PG_Kernel
63
64 void (*local_kernel)(T*, T*, T*, T*, T*); // Pointer to local kernel function
65
66 Hermite<T> herm_tools; // Hermite tools object
67
68
69 public:
70
77
78 Hermite<T> herm_tools;
79
80 int_order = 25;
81 num_particles = 1000;
82 margin = 5.0;
83 ac_length = 1.0;
84 mu_gas = 4.0;
85 T_gas = 400.0;
86
87 }
88
101 PG_Kernel(Surface<T> surface, Gas<T> gas, void (*local_kernel)(T*, T*, T*, T*, T*), Matrix<T> incident_velocity, unsigned long num_particles, std::string sim_name) {
102
103 int_order = 25;
104 margin = 5.0;
105
106 this->surface = surface;
107 this->gas = gas;
108 this->ac_length = surface.autocorrelation_length();
109
110 roots = herm_tools.get_roots(int_order);
111 weights = herm_tools.get_weights(int_order);
112
113 std::tuple<Matrix<T>, Matrix<T>> roots_tuple = roots.meshgrid(roots);
114 roots_2d_r = std::get<1>(roots_tuple);
115 roots_2d_c = std::get<0>(roots_tuple);
116 weights_2d = weights.outer(weights);
117
118 gamma = roots;
119 gamma_dot = roots;
120 gamma_2d = roots_2d_r;
121 gamma_dot_2d = roots_2d_c;
122
123 mu = herm_tools.evaluate(surface.get_mu_coefficients(), gamma);
124 sigma = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma);
125 mu_2d = herm_tools.evaluate(surface.get_mu_coefficients(), gamma_2d);
126 sigma_2d = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma_2d);
127 mu_gamma = herm_tools.evaluate_deriv(surface.get_mu_coefficients(), gamma);
128 sigma_gamma = herm_tools.evaluate_deriv(surface.get_sigma_coefficients(), gamma);
129 mu_gamma_2d = herm_tools.evaluate_deriv(surface.get_mu_coefficients(), gamma_2d);
130 sigma_gamma_2d = herm_tools.evaluate_deriv(surface.get_sigma_coefficients(), gamma_2d);
131
132 p_gamma = exp(gamma * gamma * (-1.0) / 2.0) * 1.0 / std::sqrt(2.0 * M_PI);
133 p_gamma_2d = exp(gamma_2d * gamma_2d * (-1.0) / 2.0) * 1.0 / std::sqrt(2.0 * M_PI);
134 p_gamma_dot = exp(gamma_dot * gamma_dot * (-1.0) * ac_length * ac_length / 4.0) * ac_length / std::sqrt(4.0 * M_PI);
135 p_gamma_dot_2d = exp(gamma_dot_2d * gamma_dot_2d * (-1.0) * ac_length * ac_length / 4.0) * ac_length / std::sqrt(4.0 * M_PI);
136
137 mu_gas = gas.get_molar_mass();
138 T_gas = gas.get_temperature();
139
140 this->local_kernel = local_kernel;
141
142 this->incident_velocity = incident_velocity;
143 this->num_particles = num_particles;
144 this->sim_name = sim_name;
145
146 }
147
154
155 int_order = 0;
156 num_particles = 0;
157 margin = 0.0;
158 ac_length = 0.0;
159 mu_gas = 0.0;
160 T_gas = 0.0;
161 this->local_kernel = nullptr;
162
163 }
164
170 Gas<T>& get_gas() {return this->gas;}
171
177 Surface<T>& get_surface() {return this->surface;}
178
189 T random_sampler_1D(T (PG_Kernel::*pdf)(T, std::vector<T>), std::vector<T> parameters, T x0, T size, T domain_x[]);
190
204 std::tuple<T, T> random_sampler_2D(T (PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]);
205
219 std::tuple<T, T> random_sampler_angles(T (PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]);
220
229 T mixture_pdf(T gamma, T gamma_dot, std::vector<T> parameters);
230
239 T kirchhoff_pdf(T theta_r1, T theta_r2, std::vector<T> parameters);
240
249 T shadow(T theta_r1, T height, int direction);
250
260 T occlusion(T theta_r1, T height, T distance, int direction);
261
269 T height_0_pdf(T height_0, std::vector<T> parameters);
270
278 T distance_pdf(T distance, std::vector<T> parameters);
279
289 Matrix<T> global_to_local(Matrix<T> global_vector, T theta_i, T theta_r1, T theta_r2);
290
300 Matrix<T> local_to_global(Matrix<T> local_vector, T theta_i, T theta_r1, T theta_r2);
301
310 Matrix<T> angles_to_slopes(T theta_i, T theta_r1, T theta_r2);
311
318
324 std::vector<trajectory<T>> sample_batch();
325
332 void save(std::vector<trajectory<T>> trajectories, std::string filename);
333
340 void combine_files(std::string filename, int num_procs);
341
347 void set_num_particles(unsigned long num_particles) {this->num_particles = num_particles;}
348
355 std::vector<trajectory<T>> import_data(std::string filename);
356
362 std::string& get_sim_name() {return sim_name;}
363
369 unsigned long get_num_particles() {return num_particles;}
370
371
372};
373
374template <typename T>
375T PG_Kernel<T>::random_sampler_1D(T (PG_Kernel::*pdf)(T, std::vector<T>), std::vector<T> parameters, T x0, T size, T domain_x[]) {
376
377 unsigned long steps = 100;
378
379 T x = x0, x_p = 0.0, acceptance = 0.0;
380
381 std::random_device rd;
382 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
383 std::mt19937 rng(fullSeed);
384 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
385 std::normal_distribution<T> normDist(0.0f, 1.0f);
386
387 for(auto i = 0; i < steps; i ++) {
388
389 T step = normDist(rng) * size;
390 T u_sample = uniformDist(rng);
391 x_p = x + step;
392 acceptance = 0.0;
393 if(x_p >= domain_x[0] && x_p <= domain_x[1]) {
394 acceptance = std::min(1.0, (this->*pdf)(x_p, parameters) / (this->*pdf)(x, parameters));
395 }
396 if(acceptance > u_sample) {
397 x = x_p;
398 }
399 }
400
401 return x;
402}
403
404template <typename T>
405std::tuple<T, T> PG_Kernel<T>::random_sampler_2D(T (PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[]) {
406
407 unsigned long steps = 100;
408
409 T x = x0, x_p, y_p;
410 T y = y0;
411 T acceptance = 0.0;
412
413 std::random_device rd;
414 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
415 std::mt19937 rng(fullSeed);
416 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
417 std::normal_distribution<T> normDist(0.0f, 1.0f);
418
419 for(auto i = 0; i < steps; i ++) {
420
421 T step_x = normDist(rng) * size_x;
422 T step_y = normDist(rng) * size_y;
423 T u_sample = uniformDist(rng);
424 x_p = x + step_x;
425 y_p = y + step_y;
426 acceptance = 0.0;
427 if(x_p >= domain_x[0] && x_p <= domain_x[1] && y_p >= domain_y[0] && y_p <= domain_y[1]) {
428 acceptance = std::min(1.0, (this->*pdf)(x_p, y_p, parameters) / (this->*pdf)(x, y, parameters));
429 }
430 if(acceptance > u_sample) {
431 x = x_p;
432 y = y_p;
433 }
434 }
435
436 return std::make_tuple(x, y);
437}
438
439template <typename T>
440std::tuple<T, T> PG_Kernel<T>::random_sampler_angles(T (PG_Kernel::*pdf)(T, T, std::vector<T>), std::vector<T> parameters, T theta_r1, T theta_r2, T size_x, T size_y, T domain_x[], T domain_y[]) {
441
442 unsigned long steps = 100;
443
444 std::random_device rd;
445 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
446 std::mt19937 rng(fullSeed);
447 std::uniform_real_distribution<T> uniformDist(0.0f, 1.0f);
448 std::normal_distribution<T> normDist(0.0f, 1.0f);
449
450 T acceptance = 0.0;
451
452 for(auto i = 0; i < steps; i ++) {
453
454 T step_x = normDist(rng) * size_x;
455 T step_y = normDist(rng) * size_y;
456 T u_sample = uniformDist(rng);
457 T theta_r1_p = theta_r1 + step_x;
458 T theta_r2_p = theta_r2 + step_y;
459 acceptance = 0.0;
460 if(theta_r1_p < 0.0) {
461 theta_r1_p = std::abs(theta_r1_p);
462 theta_r2_p += M_PI;
463 }
464 if(theta_r2_p < 0.0) {
465 theta_r2_p += 2.0 * M_PI;
466 }
467 if(theta_r2_p > 2.0 * M_PI) {
468 theta_r2_p -= 2.0 * M_PI;
469 }
470
471 if(theta_r1_p < domain_x[1]) {
472 acceptance = std::min(1.0, (this->*pdf)(theta_r1_p, theta_r2_p, parameters) / (this->*pdf)(theta_r1, theta_r2, parameters));
473 }
474
475 if(acceptance > u_sample) {
476 theta_r1 = theta_r1_p;
477 theta_r2 = theta_r2_p;
478 }
479 }
480
481 return std::make_tuple(theta_r1, theta_r2);
482}
483
484template <typename T>
485T PG_Kernel<T>::mixture_pdf(T gamma_local, T gamma_dot_local, std::vector<T> parameters) {
486
487 bool first = (bool)parameters[0];
488 T theta_i = parameters[1];
489 T distance = parameters[2];
490
491 T mu_local = herm_tools.evaluate(surface.get_mu_coefficients(), gamma_local);
492 T mu_dot_local = herm_tools.evaluate_deriv(surface.get_mu_coefficients(), gamma_local) * gamma_dot_local;
493 T sigma_local = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma_local);
494 T sigma_gamma_local = herm_tools.evaluate_deriv(surface.get_sigma_coefficients(), gamma_local);
495
496 T w = std::sqrt(2.0 * (sigma_local * sigma_local / ac_length / ac_length) + sigma_gamma_local * sigma_gamma_local * gamma_dot_local * gamma_dot_local);
497 T eta = 1.0 / std::tan(theta_i);
498 T slope_term = 0.5 * (1.0 + std::erf((eta + mu_dot_local) / (w * std::sqrt(2.0)))) * (std::cos(theta_i) + mu_dot_local * std::sin(theta_i)) +
499 w * std::sin(theta_i) / std::sqrt(2.0 * M_PI) * std::exp(- (eta + mu_dot_local) * (eta + mu_dot_local) / 2.0 / w / w);
500 T p_gamma_local = 1.0 / std::sqrt(2.0 * M_PI) * std::exp( - gamma_local * gamma_local / 2.0);
501 T p_gamma_dot_local = ac_length / std::sqrt(4.0 * M_PI) * std::exp(- gamma_dot_local * gamma_dot_local * ac_length * ac_length / 4.0);
502
503 if(first) {
504 return p_gamma_local * p_gamma_dot_local * slope_term * shadow(theta_i, mu_local, -1);
505 }
506 else {
507 return p_gamma_local * p_gamma_dot_local * slope_term * occlusion(theta_i, mu_local, distance, -1);
508 }
509
510}
511
512template <typename T>
513T PG_Kernel<T>::kirchhoff_pdf(T theta_r1, T theta_r2, std::vector<T> parameters) {
514
515 T theta_i = parameters[0];
516 T gamma = parameters[1];
517 T gamma_dot_x = parameters[2];
518 T gamma_dot_y = parameters[3];
519 T slope_x = parameters[4];
520 T slope_y = parameters[5];
521 T tau = parameters[6];
522
523 T C = std::exp(- tau * tau / ac_length / ac_length);
524
525 T mu_gamma = herm_tools.evaluate_deriv(surface.get_mu_coefficients(), gamma);
526 T sigma = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma);
527 T sigma_gamma = herm_tools.evaluate_deriv(surface.get_sigma_coefficients(), gamma);
528
529 T mu_kirchhoff_x = mu_gamma * gamma_dot_x;
530 T mu_kirchhoff_y = mu_gamma * gamma_dot_y;
531
532 T vx = std::sin(theta_i) - std::sin(theta_r1) * std::cos(theta_r2);
533 T vy = - std::sin(theta_r1) * std::sin(theta_r2);
534 T vz = - std::cos(theta_i) - std::cos(theta_r1);
535
536 T F = (1.0 + std::cos(theta_i) * std::cos(theta_r1) - std::sin(theta_i) * std::sin(theta_r1) * std::cos(theta_r2)) /
537 (std::cos(theta_i) * std::cos(theta_i) + std::cos(theta_r1) * std::cos(theta_i));
538
539 T sigma_kirchhoff_x = std::sqrt((sigma * sigma / ac_length / ac_length) + 0.5 * sigma_gamma * sigma_gamma * gamma_dot_x * gamma_dot_x);
540 T sigma_kirchhoff_y = std::sqrt((sigma * sigma / ac_length / ac_length) + 0.5 * sigma_gamma * sigma_gamma * gamma_dot_y * gamma_dot_y);
541
542 T vx_k = (vx + mu_kirchhoff_x * vz);
543 T vy_k = (vy + mu_kirchhoff_y * vz);
544
545 T p_kirchhoff_x = std::sqrt(2.0 * M_PI) * F / vz / std::sqrt(1.0 - C * C) *
546 std::exp(- (vx_k * vx_k - 2.0 * C * vx_k * (slope_x - mu_kirchhoff_x) * vz + vz * vz * (slope_x - mu_kirchhoff_x) * (slope_x - mu_kirchhoff_x))
547 / (4.0 * vz * vz * sigma_kirchhoff_x * sigma_kirchhoff_x * (1.0 - C * C)));
548 T p_kirchhoff_y = std::sqrt(2.0 * M_PI) * F / vz / std::sqrt(1.0 - C * C) *
549 std::exp(- (vy_k * vy_k - 2.0 * C * vy_k * (slope_y - mu_kirchhoff_y) * vz + vz * vz * (slope_y - mu_kirchhoff_y) * (slope_y - mu_kirchhoff_y))
550 / (4.0 * vz * vz * sigma_kirchhoff_y * sigma_kirchhoff_y * (1.0 - C * C)));
551
552 return p_kirchhoff_x * p_kirchhoff_y * std::sin(theta_r1);
553}
554
555template <typename T>
556T PG_Kernel<T>::shadow(T theta_r1, T height, int direction) {
557
558 // Check if theta_r1 is greater than 90 degrees
559 if(theta_r1 > M_PI / 2.0) {
560 return 0.0;
561 }
562 T EPSILON = 1e-15;
563 theta_r1 = std::max(EPSILON, theta_r1);
564 T eta = 1.0 / std::tan(theta_r1);
565 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
566 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
567 (mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) * direction + eta) / (w_2d * std::sqrt(2.0)))));
568 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
569 T base = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
570 base = std::min(1.0, base);
571 return std::pow(base, exponent);
572}
573
574template <typename T>
575T PG_Kernel<T>::occlusion(T theta_r1, T height, T distance, int direction) {
576
577 theta_r1 = std::max(1e-12, theta_r1);
578 T eta = 1.0 / std::tan(theta_r1);
579 T height_new = height + distance * eta;
580 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
581 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
582 (mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) / (w_2d * std::sqrt(2.0)))));
583 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
584 T base_old = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
585 T base_new = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_new) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
586
587 return std::pow(base_old / base_new, exponent);
588}
589
590template <typename T>
591T PG_Kernel<T>::distance_pdf(T distance, std::vector<T> parameters) {
592
593 T theta_r1 = std::max(parameters[0], 1e-12);
594 T height_0 = parameters[1];
595 T gamma_0 = parameters[2];
596 T gamma_dot_0 = parameters[3];
597
598 T eta = 1.0 / std::tan(theta_r1);
599 T height_new = height_0 + eta * distance;
600 T mu_0 = herm_tools.evaluate(surface.get_mu_coefficients(), gamma_0);
601 T sigma_0 = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma_0);
602 Matrix<T> w_2d = (sigma_2d * sigma_2d / ac_length / ac_length * 2.0 + sigma_gamma_2d * sigma_gamma_2d * gamma_dot_2d * gamma_dot_2d) ^ 0.5;
603 Matrix<T> Delta = (w_2d / std::sqrt(2.0 * M_PI) * exp(((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) ^ 2.0) / ((w_2d * w_2d) * 2.0) * (-1.0)) -
604 (mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) * 0.5 * (erfc((mu_gamma_2d * gamma_dot_2d * (-1.0) + eta) / (w_2d * std::sqrt(2.0)))));
605 T exponent = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * Delta) / eta;
606 T base_old = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_0) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
607 T base_new = sum(weights_2d * p_gamma_2d * p_gamma_dot_2d * 0.5 * (erfc((mu_2d * (-1.0) + height_new) / (sigma_2d * std::sqrt(2.0))) * (-1.0) + 2.0));
608 T p_h = 1.0 / std::sqrt(2.0 * M_PI) / sigma_0 * std::exp(- (height_new - mu_0) * (height_new - mu_0) / 2.0 / sigma_0 / sigma_0);
609 T S = std::pow(base_old, exponent) / std::pow(base_new, exponent + 1.0);
610
611 return S * p_h;
612}
613
614template <typename T>
615T PG_Kernel<T>::height_0_pdf(T height_0, std::vector<T> parameters) {
616
617 T theta_i = parameters[0];
618 T gamma_0 = parameters[1];
619 T gamma_dot_0 = parameters[2];
620
621 T mu_0 = herm_tools.evaluate(surface.get_mu_coefficients(), gamma_0);
622 T sigma_0 = herm_tools.evaluate(surface.get_sigma_coefficients(), gamma_0);
623
624 T p_h = 1.0 / std::sqrt(2.0 * M_PI) / sigma_0 * std::exp(- std::pow(height_0 - mu_0, 2.0) / 2.0 / sigma_0 / sigma_0);
625
626 return p_h * shadow(theta_i, height_0, -1);
627}
628
629template <typename T>
630Matrix<T> PG_Kernel<T>::global_to_local(Matrix<T> global_vector, T theta_i, T theta_r1, T theta_r2) {
631
632 Matrix<T> unit_incident = global_vector / global_vector.norm();
633 Matrix<T> z_local(3, 1);
634 z_local(0, 0) = std::sin(theta_r1) * std::cos(theta_r2) - std::sin(theta_i);
635 z_local(1, 0) = std::sin(theta_r1) * std::sin(theta_r2);
636 z_local(2, 0) = std::cos(theta_r1) + std::cos(theta_i);
637 z_local /= z_local.norm();
638 Matrix<T> x_local(3, 1);
639 x_local = unit_incident - z_local * unit_incident.tr().dot(z_local)(0, 0);
640 x_local /= x_local.norm();
641 Matrix<T> y_local(3, 1);
642 y_local = x_local.cross(z_local);
643 y_local /= y_local.norm();
644
645 T res_x = global_vector.tr().dot(x_local)(0, 0);
646 T res_y = global_vector.tr().dot(y_local)(0, 0);
647 T res_z = global_vector.tr().dot(z_local)(0, 0);
648
649 Matrix<T> res(3, 1, {res_x, res_y, res_z});
650
651 return res;
652}
653
654template <typename T>
655Matrix<T> PG_Kernel<T>::local_to_global(Matrix<T> local_vector, T theta_i, T theta_r1, T theta_r2) {
656
657 Matrix<T> unit_incident(3, 1, {std::sin(theta_i), 0.0, -std::cos(theta_i)});
658 Matrix<T> z_local(3, 1);
659 z_local(0, 0) = std::sin(theta_r1) * std::cos(theta_r2) - std::sin(theta_i);
660 z_local(1, 0) = std::sin(theta_r1) * std::sin(theta_r2);
661 z_local(2, 0) = std::cos(theta_r1) + std::cos(theta_i);
662 z_local /= z_local.norm();
663 Matrix<T> x_local(3, 1);
664 x_local = unit_incident - z_local * unit_incident.tr().dot(z_local)(0, 0);
665 x_local /= x_local.norm();
666 Matrix<T> y_local(3, 1);
667 y_local = x_local.cross(z_local);
668 y_local /= y_local.norm();
669
670 return x_local * local_vector(0, 0) + y_local + local_vector(1, 0) + z_local * local_vector(2, 0);
671}
672
673template <typename T>
674Matrix<T> PG_Kernel<T>::angles_to_slopes(T theta_i, T theta_r1, T theta_r2) {
675
676 Matrix<T> unit_incident(3, 1, {std::sin(theta_i), 0.0, -std::cos(theta_i)});
677 Matrix<T> unit_reflected(3, 1, {std::sin(theta_r1) * std::cos(theta_r2), std::sin(theta_r1) * std::sin(theta_r2), std::cos(theta_r1)});
678 Matrix<T> normal(3, 1);
679 normal = unit_reflected - unit_incident;
680 normal /= std::abs(normal(2, 0));
681
682 return normal;
683}
684
685template <typename T>
687
688 Matrix<T> incident_velocity_copy = this->incident_velocity;
689
690 T incident_norm = incident_velocity_copy.norm();
691 T theta_i = std::acos(- incident_velocity_copy(2, 0) / incident_norm);
692 T theta_r1 = 0.0, theta_r2 = 0.0, theta_r2_total = 0.0, distance = 1e8;
693 T gamma_0 = 0.0, gamma_dot_0_x = 0.0, gamma_dot_0_y = 0.0, slope_x = 0.0, slope_y = 0.0, slope_x_new =0.0, slope_y_new = 0.0;
694
695 long num_col = 0;
696 bool collided = true;
697
698 std::vector<Matrix<T>> p_positions;
699 std::vector<Matrix<T>> p_velocities;
700
701 std::random_device rd;
702 std::seed_seq fullSeed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
703 std::mt19937 rng(fullSeed);
704 std::normal_distribution<T> normDist(0.0f, 1.0f);
705
706 T gamma, gamma_dot_x, gamma_dot_y;
707 T domain[] = {- margin, margin};
708 std::tie(gamma, gamma_dot_x) = random_sampler_2D(&PG_Kernel::mixture_pdf, {1.0, theta_i, 0.0}, normDist(rng), normDist(rng), 0.5, 0.5, domain, domain);
709 gamma_dot_y = normDist(rng) * std::sqrt(2.0) / ac_length;
710
711 T height_domain[] = {-margin, margin};
712 T height_0 = random_sampler_1D(&PG_Kernel::height_0_pdf, {theta_i, gamma, gamma_dot_x}, normDist(rng), 0.5, height_domain);
713 p_positions.push_back(Matrix<T>(3, 1, {height_0, 0.0, 0.0}));
714
715 while(collided) {
716
717 T angle_domain_x[] = {0.0, M_PI - theta_i};
718 T angle_domain_y[] = {-M_PI, M_PI};
719 std::tie(theta_r1, theta_r2) = random_sampler_angles(&PG_Kernel::kirchhoff_pdf, {theta_i, gamma, gamma_dot_x, gamma_dot_y, slope_x, slope_y, distance}, 1e-5, 1e-5,
720 (M_PI - theta_i) * 0.25, (M_PI - theta_i) * 0.25, angle_domain_x, angle_domain_y);
721 Matrix<T> slopes_matrix = angles_to_slopes(theta_i, theta_r1, theta_r2);
722 slope_x_new = - slopes_matrix(0, 0);
723 slope_y_new = - slopes_matrix(1, 0);
724
725 Matrix<T> incident_velocity_local = global_to_local(incident_velocity_copy, theta_i, theta_r1, theta_r2);
726 T vel_in_array[] = {incident_velocity_local(0, 0), incident_velocity_local(1, 0), incident_velocity_local(2, 0)};
727 T vel_refl_array[] = {0.0, 0.0, 0.0};
728 T gsi_parameters[] = {surface.get_local_parameters()[1], surface.get_local_parameters()[2]};
729 T surface_parameters[] = {surface.get_local_parameters()[0]};
730 T gas_parameters[] = {gas.get_temperature(), gas.get_density(), gas.get_molar_mass(), gas.get_speed()};
731
732 local_kernel(vel_in_array, vel_refl_array, gsi_parameters, gas_parameters, surface_parameters);
733
734 Matrix<T> reflected_velocity_local(3, 1, {vel_refl_array[0], vel_refl_array[1], vel_refl_array[2]});
735 Matrix<T> reflected_velocity = local_to_global(reflected_velocity_local, theta_i, theta_r1, theta_r2);
736
737 theta_r1 = std::acos(reflected_velocity(2, 0) / reflected_velocity.norm());
738 theta_r2 = std::atan2(reflected_velocity(1, 0), reflected_velocity(0, 0));
739
740 std::bernoulli_distribution shadow_pdf(1.0 - shadow(theta_r1, height_0, 1));
741 collided = shadow_pdf(rng);
742
743 T distance_domain[] = {0.0, margin * 5.0};
744 distance = random_sampler_1D(&PG_Kernel::distance_pdf, {theta_r1, height_0, gamma, gamma_dot_x}, 0.0, 0.5, distance_domain);
745 height_0 += distance / std::tan(theta_r1);
746 theta_i = M_PI - theta_r1;
747 theta_r2_total += theta_r2;
748
749 incident_norm = reflected_velocity.norm();
750 incident_velocity_copy(0, 0) = incident_norm * std::sin(theta_i);
751 incident_velocity_copy(1, 0) = 0.0;
752 incident_velocity_copy(2, 0) = - incident_norm * std::cos(theta_i);
753
754 Matrix<T> global_reflected_velocity(3, 1, {std::sin(theta_r1) * std::cos(theta_r2_total),
755 std::sin(theta_r1) * std::sin(theta_r2_total),
756 std::cos(theta_r1)});
757 global_reflected_velocity *= incident_norm;
758 p_velocities.push_back(global_reflected_velocity);
759
760 Matrix<T> position_change(3, 1, {std::sin(theta_r1) * std::cos(theta_r2_total),
761 std::sin(theta_r1) * std::sin(theta_r2_total), std::cos(theta_r1)});
762 position_change *= distance / std::sin(theta_r1);
763 if (collided) p_positions.push_back(p_positions.back() + position_change);
764
765 slope_x = - std::cos(theta_r2) * slope_x_new - std::sin(theta_r2) * slope_y_new;
766 slope_y = std::sin(theta_r2) * slope_x_new - std::cos(theta_r2) * slope_y_new;
767
768 std::tie(gamma, gamma_dot_x) = random_sampler_2D(&PG_Kernel::mixture_pdf, {0.0, theta_i, distance}, normDist(rng), normDist(rng), 0.5, 0.5, domain, domain);
769 gamma_dot_y = normDist(rng) * std::sqrt(2.0) / ac_length;
770 num_col ++;
771 }
772 trajectory<T> p_trajectory = {p_velocities, p_positions};
773
774 return p_trajectory;
775}
776
777template <typename T>
778std::vector<trajectory<T>> PG_Kernel<T>::sample_batch() {
779
780 std::vector<trajectory<T>> trajectories;
781
782 for(auto i = 0; i < this->num_particles; i ++) {
783 trajectory<T> p_trajectory = sample_particle();
784 trajectories.push_back(p_trajectory);
785 }
786 return trajectories;
787}
788
789template <typename T>
790void PG_Kernel<T>::save(std::vector<trajectory<T>> trajectories, std::string filename) {
791
792 std::ofstream myfile;
793 myfile.open(filename);
794
795 unsigned long num_particles = trajectories.size();
796
797 for(auto i = 0; i < num_particles; i ++) {
798
799 trajectory<T> p_trajectory = trajectories[i];
800 unsigned long num_bounces = p_trajectory.positions.size();
801
802 myfile << num_bounces << "\n";
803
804 for(auto j = 0; j < 3; j ++) {
805 for(auto k = 0; k < num_bounces; k ++) {
806 myfile << p_trajectory.positions[k](j, 0) << " ";
807 }
808 myfile << "\n";
809 }
810
811 for(auto j = 0; j < 3; j ++) {
812 for(auto k = 0; k < num_bounces; k ++) {
813 myfile << p_trajectory.velocities[k](j, 0) << " ";
814 }
815 myfile << "\n";
816 }
817
818 }
819 myfile.close();
820
821}
822
823template <typename T>
824std::vector<trajectory<T>> PG_Kernel<T>::import_data(std::string filename) {
825
826 std::ifstream linecounter(filename);
827 std::string line;
828 std::vector<trajectory<T>> trajectories;
829
830 long num_lines = 0;
831
832 while( std::getline(linecounter, line)) {
833 ++ num_lines;
834 }
835
836 long num_particles = num_lines / 7;
837 long num_bounces;
838
839 std::ifstream myfile(filename);
840
841 for(long i = 0; i < num_particles; i ++) {
842
843 myfile >> num_bounces;
844 Matrix<T> position(3, 1);
845 Matrix<T> velocity(3, 1);
846 std::vector<Matrix<T>> p_positions(num_bounces, Matrix<T>(3, 1));
847 std::vector<Matrix<T>> p_velocities(num_bounces, Matrix<T>(3, 1));
848
849 for(long k = 0; k < 3; k ++) {
850 for(long j = 0; j < num_bounces; j ++) {
851 myfile >> p_positions[j](k, 0);
852 }
853 }
854 for(long k = 0; k < 3; k ++) {
855 for(long j = 0; j < num_bounces; j ++) {
856 myfile >> p_velocities[j](k, 0);
857 }
858 }
859 trajectory<T> trajectory = {p_velocities, p_positions};
860 trajectories.push_back(trajectory);
861 }
862
863 return trajectories;
864}
865
866template <typename T>
867void PG_Kernel<T>::combine_files(std::string filename, int num_procs) {
868
869 std::vector<trajectory<T>> trajectories_total;
870
871 for( long i = 0; i < num_procs; i ++) {
872 std::vector<trajectory<T>> trajectories = import_data(filename + "_" + std::to_string(i) + ".dat");
873 trajectories_total.insert(trajectories_total.end(), trajectories.begin(), trajectories.end());
874 std::remove((filename + "_" + std::to_string(i) + ".dat").c_str());
875 }
876 save(trajectories_total, filename + ".dat");
877}
Class representing a Gas with properties such as density, temperature, molar mass,...
Definition Gas.h:10
T get_molar_mass() const
Returns the molar mass of the gas.
Definition Gas.h:141
T get_density() const
Returns the density of the gas.
Definition Gas.h:134
T get_temperature() const
Returns the temperature of the gas.
Definition Gas.h:127
T get_speed() const
Returns the speed of the gas.
Definition Gas.h:148
Class for computing Hermite polynomials and related operations.
Definition Hermite_Tools.h:11
Matrix< T > evaluate(Matrix< T > coefficients, Matrix< T > values)
Evaluates the Hermite polynomial at specified values.
Definition Hermite_Tools.h:136
Matrix< T > get_roots(unsigned int order)
Retrieves the roots of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:276
Matrix< T > get_weights(unsigned int order)
Retrieves the weights of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:299
Matrix< T > evaluate_deriv(Matrix< T > coefficients, Matrix< T > values)
Evaluates the derivative of the Hermite polynomial at specified values.
Definition Hermite_Tools.h:215
This class implements a matrix object used for linear algebra and vectorized operations.
Definition Matrix.h:25
Matrix cross(const Matrix &other)
Calculates the cross product of two matrices.
Definition Matrix.h:624
std::tuple< Matrix< T >, Matrix< T > > meshgrid(const Matrix &other)
Calculates the meshgrid of 2 1D matrices.
Definition Matrix.h:656
Matrix dot(const Matrix &other)
Calculates the dot product of two matrices.
Definition Matrix.h:604
T norm()
Calculates the Frobenius norm of the matrix.
Definition Matrix.h:674
Matrix outer(const Matrix &other)
Calculates the outer product of two matrices.
Definition Matrix.h:639
Matrix tr()
Transposes the matrix.
Definition Matrix.h:705
Represents a class for the proposed rough surface kernel applied to Poly-Gaussian surfaces.
Definition PG_Kernel.h:26
trajectory< T > sample_particle()
Sample a single particle trajectory.
Definition PG_Kernel.h:686
void set_num_particles(unsigned long num_particles)
Set the number of particles.
Definition PG_Kernel.h:347
T distance_pdf(T distance, std::vector< T > parameters)
Calculate the PDF for the distance parameter.
Definition PG_Kernel.h:591
unsigned long get_num_particles()
Get the number of particles.
Definition PG_Kernel.h:369
void save(std::vector< trajectory< T > > trajectories, std::string filename)
Save the trajectories to a file.
Definition PG_Kernel.h:790
PG_Kernel(Surface< T > surface, Gas< T > gas, void(*local_kernel)(T *, T *, T *, T *, T *), Matrix< T > incident_velocity, unsigned long num_particles, std::string sim_name)
Constructor for PG_Kernel.
Definition PG_Kernel.h:101
std::vector< trajectory< T > > import_data(std::string filename)
Import trajectory data from a file.
Definition PG_Kernel.h:824
T mixture_pdf(T gamma, T gamma_dot, std::vector< T > parameters)
Calculate the PDF of the control process gamma and its derivative gamma_dot, for Poly-Gaussian surfac...
Definition PG_Kernel.h:485
Matrix< T > local_to_global(Matrix< T > local_vector, T theta_i, T theta_r1, T theta_r2)
Convert local coordinates to global coordinates.
Definition PG_Kernel.h:655
Matrix< T > angles_to_slopes(T theta_i, T theta_r1, T theta_r2)
Convert angles to slopes.
Definition PG_Kernel.h:674
std::string & get_sim_name()
Get the simulation name.
Definition PG_Kernel.h:362
std::tuple< T, T > random_sampler_angles(T(PG_Kernel::*pdf)(T, T, std::vector< T >), std::vector< T > parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[])
Perform random sampling of angles in 2D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:440
T shadow(T theta_r1, T height, int direction)
Calculate the probability of a reflected particle to escape the surface.
Definition PG_Kernel.h:556
T random_sampler_1D(T(PG_Kernel::*pdf)(T, std::vector< T >), std::vector< T > parameters, T x0, T size, T domain_x[])
Perform random sampling in 1D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:375
T occlusion(T theta_r1, T height, T distance, int direction)
Calculate the probability of a reflected particle to escape the surface, given a certain travelled ho...
Definition PG_Kernel.h:575
~PG_Kernel()
Destructor for PG_Kernel.
Definition PG_Kernel.h:153
T height_0_pdf(T height_0, std::vector< T > parameters)
Calculate the PDF for the height_0 parameter.
Definition PG_Kernel.h:615
Matrix< T > global_to_local(Matrix< T > global_vector, T theta_i, T theta_r1, T theta_r2)
Convert global coordinates to local coordinates.
Definition PG_Kernel.h:630
std::vector< trajectory< T > > sample_batch()
Sample a batch of particle trajectories.
Definition PG_Kernel.h:778
std::tuple< T, T > random_sampler_2D(T(PG_Kernel::*pdf)(T, T, std::vector< T >), std::vector< T > parameters, T x0, T y0, T size_x, T size_y, T domain_x[], T domain_y[])
Perform random sampling in 2D using the Metropolis-Hastings algorithm.
Definition PG_Kernel.h:405
Surface< T > & get_surface()
Get the Surface object.
Definition PG_Kernel.h:177
T kirchhoff_pdf(T theta_r1, T theta_r2, std::vector< T > parameters)
Calculate the reflected particle angular PDF according to the Kirchhoff model.
Definition PG_Kernel.h:513
Gas< T > & get_gas()
Get the Gas object.
Definition PG_Kernel.h:170
PG_Kernel()
Default constructor for PG_Kernel.
Definition PG_Kernel.h:76
void combine_files(std::string filename, int num_procs)
Combine multiple trajectory files into one.
Definition PG_Kernel.h:867
A class representing a surface with geometry and statistical properties.
Definition Surface.h:11
T autocorrelation_length()
Retrieves the autocorrelation length of the surface.
Definition Surface.h:100
T * get_local_parameters()
Retrieves the local parameters of the surface.
Definition Surface.h:107
Matrix< T > get_sigma_coefficients()
Retrieves the sigma coefficients for the surface.
Definition Surface.h:86
Matrix< T > get_mu_coefficients()
Retrieves the mu coefficients for the surface.
Definition Surface.h:79
Represents a trajectory of a certain type.
Definition PG_Kernel.h:10