diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 0da25a4..ca25877 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -1,2 +1,3 @@ c765e51344ad4a458d042ceefdf19fd61edf79ab # clang-format 59fe44fd6d08218b8aeae7612bdd5a71d0289ce5 # format /test/test_hostwrite_initialization.cpp +570b3252bcb6e080bde4113cc9d7ed2aa5796fd6 # format degas2 physics files diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml index 81cc715..d72e906 100644 --- a/.github/workflows/format-check.yml +++ b/.github/workflows/format-check.yml @@ -2,8 +2,9 @@ name: Format Check on: push: + branches: [main] pull_request: - pull_request_target: + branches: [main] workflow_dispatch: jobs: @@ -31,6 +32,11 @@ jobs: exit 0 fi + echo "Files to be checked:" + for file in "${files[@]}"; do + echo " $file" + done + clang-format-20 --version clang-format-20 -style=file --dry-run -Werror "${files[@]}" diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml index b1c6eb0..fc0d707 100644 --- a/.github/workflows/static-analysis.yml +++ b/.github/workflows/static-analysis.yml @@ -1,8 +1,10 @@ name: Static Analysis on: - push: - pull_request: + workflow_run: + workflows: ["Test"] + types: + - completed workflow_dispatch: env: @@ -11,9 +13,14 @@ env: jobs: clang-tidy: runs-on: ubuntu-latest + if: > + github.event_name == 'workflow_dispatch' || + (github.event.workflow_run.conclusion == 'success') steps: - uses: actions/checkout@v6 + with: + ref: ${{ github.event.workflow_run.head_sha || github.ref }} - name: Install dependencies uses: ./.github/actions/install-deps @@ -35,6 +42,7 @@ jobs: -DCMAKE_CXX_COMPILER=mpicxx \ -DCMAKE_C_COMPILER=mpicc \ -DBUILD_SHARED_LIBS=ON \ + -DPUMITALLY_ENABLE_DEGAS2=ON \ -DCMAKE_CXX_CLANG_TIDY=${TIDY_CMD} - name: Build with clang-tidy diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index df30e1e..72ef51a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -2,11 +2,14 @@ name: Test on: push: + branches: [main] pull_request: + branches: [main] workflow_dispatch: env: DEPS_DIR: ${{ github.workspace }}/deps + LD_LIBRARY_PATH: ${{ github.workspace }}/deps/lib jobs: ctest: @@ -26,12 +29,17 @@ jobs: -DCMAKE_PREFIX_PATH=$DEPS_DIR \ -DCMAKE_CXX_COMPILER=mpicxx \ -DCMAKE_C_COMPILER=mpicc \ - -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_BUILD_TYPE=Debug \ -DBUILD_SHARED_LIBS=ON \ - -DPUMITALLYOPENMC_ENABLE_TESTS=ON + -DPUMITALLY_ENABLE_DEGAS2=ON \ + -DPUMITALLY_ENABLE_TESTS=ON - name: Build run: cmake --build build -j$(nproc) - name: Run tests - run: ctest --test-dir build --output-on-failure + run: | + export OMP_PROC_BIND=spread + export OMP_PLACES=threads + export OMP_NUM_THREADS=$(nproc) + ctest --test-dir build --output-on-failure diff --git a/.gitignore b/.gitignore index 8cce1f8..3549186 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ build-* .idea/ config.sh Testing/ +*.vtk diff --git a/CMakeLists.txt b/CMakeLists.txt index 42f3055..6f2635b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,7 @@ option(BUILD_SHARED_LIBS "Build shared libraries" ON) option(PUMITALLY_ENABLE_TESTS "Enable tests" OFF) option(PUMITALLY_USE_KOKKOS_CUDA "Use CUDA" OFF) option(PUMITALLY_MEASURE_TIME "Measure time" OFF) +option(PUMITALLY_ENABLE_DEGAS2 "Enable Degas2" OFF) ################################################################### ################################## PUMIPIC ######################## @@ -30,6 +31,18 @@ set(PUMITALLY_SOURCES src/pumitally/PumiTallyImpl.cpp ) +if(PUMITALLY_ENABLE_DEGAS2) + set(PUMITALLY_DEGAS2_HEADERS + src/physics/degas2/DG2Physics.h + src/physics/degas2/DG2CrossSection.h + ) + set(PUMITALLY_DEGAS2_SOURCES + src/physics/degas2/DG2Physics.cpp + src/physics/degas2/DG2CrossSection.cpp + ) +endif() + + ################################## Testing ######################## if(PUMITALLY_ENABLE_TESTS) @@ -73,6 +86,21 @@ if(PUMITALLY_MEASURE_TIME) target_compile_definitions(pumitally PRIVATE -DPUMITALLY_MEASURE_TIME) endif() + +if(PUMITALLY_ENABLE_DEGAS2) + add_library(pumitally_degas2) + target_sources(pumitally_degas2 PRIVATE ${PUMITALLY_DEGAS2_SOURCES}) + target_link_libraries(pumitally_degas2 PRIVATE pumipic::pumipic) + + add_executable(degas2 src/physics/degas2/degas2-main.cpp) + target_link_libraries(degas2 PRIVATE pumitally_degas2 pumitally pumipic::pumipic) + target_include_directories(degas2 + PRIVATE + $ + $ + ${CMAKE_SOURCE_DIR}/src/pumitally + ) +endif() ################################################################### @@ -137,5 +165,25 @@ install(FILES DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pumitally ) +if(PUMITALLY_ENABLE_DEGAS2) + install(TARGETS pumitally_degas2 + EXPORT pumitallyTargets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pumitally + ) + install(EXPORT pumitallyTargets + FILE pumitally_degas2Targets.cmake + NAMESPACE pumitally:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pumitally + ) + install(FILES + ${CMAKE_CURRENT_SOURCE_DIR}/src/physics/degas2/DG2Physics.h + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pumitally/physics/degas2 + ) +endif() + ################################################################### diff --git a/src/physics/degas2/DG2CrossSection.cpp b/src/physics/degas2/DG2CrossSection.cpp new file mode 100644 index 0000000..b83745a --- /dev/null +++ b/src/physics/degas2/DG2CrossSection.cpp @@ -0,0 +1,9 @@ +// +// Created by Fuad Hasans on 6/20/25. +// + +#include "DG2CrossSection.h" + +DG2CrossSection::DG2CrossSection(std::string filename) { + // fill in the arrays +} diff --git a/src/physics/degas2/DG2CrossSection.h b/src/physics/degas2/DG2CrossSection.h new file mode 100644 index 0000000..f4ed135 --- /dev/null +++ b/src/physics/degas2/DG2CrossSection.h @@ -0,0 +1,38 @@ +// +// Created by Fuad Hasan on 6/20/25. +// + +#ifndef PUMITALLY_DG2CROSSSECTION_H +#define PUMITALLY_DG2CROSSSECTION_H + +#include + +class DG2CrossSection { +public: + // fills the cross-section arrays from the given file + explicit DG2CrossSection(std::string filename); + + // explicitly given cross-section arrays constructor + DG2CrossSection(Kokkos::View sigma_t, + Kokkos::View sigma_a, + Kokkos::View scattering_matrix, + Kokkos::View sigma_s) + : nEnergyGroups(sigma_t.extent(2)), nMaterials(sigma_t.extent(0)), + nTemperatures(sigma_t.extent(1)), sigma_t(sigma_t), sigma_a(sigma_a), + scattering_matrix(scattering_matrix), sigma_s(sigma_s) {} + + DG2CrossSection(const DG2CrossSection &other) = default; + + // no need to define readers of the cross-section arrays + + int nEnergyGroups; + int nMaterials; + int nTemperatures; + + Kokkos::View sigma_t; // mat, T, g + Kokkos::View sigma_a; // mat, T, g + Kokkos::View scattering_matrix; // mat, T, g, g + Kokkos::View sigma_s; // mat, T, g +}; + +#endif // PUMITALLY_DG2CROSSSECTION_H diff --git a/src/physics/degas2/DG2Physics.cpp b/src/physics/degas2/DG2Physics.cpp new file mode 100644 index 0000000..4b986b7 --- /dev/null +++ b/src/physics/degas2/DG2Physics.cpp @@ -0,0 +1,4 @@ +// +// Created by Fuad Hasan on 6/20/25. +// +#include "DG2Physics.h" diff --git a/src/physics/degas2/DG2Physics.h b/src/physics/degas2/DG2Physics.h new file mode 100644 index 0000000..336d2c2 --- /dev/null +++ b/src/physics/degas2/DG2Physics.h @@ -0,0 +1,345 @@ +// +// Created by Fuad Hasan on 6/20/25. +// + +#ifndef PUMITALLYOPENMC_DG2PHYSICS_H +#define PUMITALLYOPENMC_DG2PHYSICS_H + +#define SEED 12345 + +#include "DG2CrossSection.h" +#include +#include + +typedef Kokkos::Random_XorShift64_Pool + random_pool_t; + +struct ParticleInfo { + double position[3]; // Position in space (x, y, z) + double direction[3]; // Direction vector (unit vector) + double weight; + int energy_group; // Energy group *index* + int particle_index; + double alpha; // To be multiplied in when computing tally +}; + +struct FieldInfo { + double electron_temperature; + double ion_temperature; + double electron_density; + double ion_density; + double bulk_flow_velocity[3]; +}; + +class DG2Physics { +public: + DG2Physics(const std::string &cross_section_file, const int num_particles, + const int seed = SEED) + : random_pool(seed), cross_section(cross_section_file) { + // Initialize particle energy + particle_energy = Kokkos::View("particle_energy", num_particles); + // cross_sections = Kokkos::View("cross_sections", + // num_particles); + } + + DG2Physics(const DG2CrossSection cross_section, const int num_particles, + const int seed = SEED) + + : random_pool(seed), cross_section(cross_section) { + // Initialize particle energy + particle_energy = Kokkos::View("particle_energy", num_particles); + // cross_sections = Kokkos::View("cross_sections", + // num_particles); + } + + KOKKOS_FUNCTION + void set_energy(ParticleInfo &particle_info, const double energy) const { + particle_energy(particle_info.particle_index) = energy; + } + + KOKKOS_FUNCTION + double ionization_cross_section(ParticleInfo &particle_info, + const FieldInfo &field_info) const { + double energy = particle_energy(particle_info.particle_index); + double mp{938.27e6 / (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2 + double particle_velocity_squared{2 * energy / mp}; // cm^2/s^2 + + // Compute Ionization Cross Section + double coef2[9]; + + coef2[0] = -3.271396786375e1; + coef2[1] = 1.353655609057e1; + coef2[2] = -5.739328757388; + coef2[3] = 1.563154982022; + coef2[4] = -2.877056004391e-1; + coef2[5] = 3.482559773737e-2; + coef2[6] = -2.631976175590e-3; + coef2[7] = 1.119543953861e-4; + coef2[8] = -2.039149852002e-6; + + double lnrate_ion{0}; // cm^2 + + for (int i = 0; i < 9; i++) { + lnrate_ion += + coef2[i] * + Kokkos::pow(Kokkos::log(field_info.electron_temperature), i); + } + double sigma_ion = + Kokkos::exp(lnrate_ion) / Kokkos::sqrt(particle_velocity_squared); + return sigma_ion; + } + + KOKKOS_FUNCTION + double charge_exchange_cross_section(ParticleInfo &particle_info, + const FieldInfo &field_info) const { + double energy = particle_energy(particle_info.particle_index); + double mp{938.27e6 / (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2 + double particle_velocity_squared{2 * energy / mp}; // cm^2/s^2 + + // Compute Charge Exchange Cross Section + double coef[9][9]; // E index, T index + + coef[0][0] = -1.829079581680e1; + coef[0][1] = 2.169137615703e-1; + coef[0][2] = 4.307131243894e-2; + coef[0][3] = -5.754895093075e-4; + coef[0][4] = -1.552077120204e-3; + coef[0][5] = -1.876800283030e-4; + coef[0][6] = 1.125490270962e-4; + coef[0][7] = -1.238982763007e-5; + coef[0][8] = 4.163596197181e-07; + + coef[1][0] = 1.640252721210e-1; + coef[1][1] = -1.106722014459e-1; + coef[1][2] = 8.948693624917e-3; + coef[1][3] = 6.062141761233e-3; + coef[1][4] = -1.210431587568e-3; + coef[1][5] = -4.052878751584e-5; + coef[1][6] = 2.875900435895e-5; + coef[1][7] = -2.616998139678e-6; + coef[1][8] = 7.558092849125e-8; + + coef[2][0] = 3.364564509137e-2; + coef[2][1] = -1.382158680424e-3; + coef[2][2] = -1.209480567154e-02; + coef[2][3] = 1.075907881928e-3; + coef[2][4] = 8.297212635856e-4; + coef[2][5] = -1.907025662962e-04; + coef[2][6] = 1.338839628570e-5; + coef[2][7] = -1.171762874107e-7; + coef[2][8] = -1.328404104165e-8; + + coef[3][0] = 9.530225559189e-3; + coef[3][1] = 7.348786286628e-3; + coef[3][2] = -3.675019470470e-4; + coef[3][3] = -8.119301728339e-4; + coef[3][4] = 1.361661816974e-4; + coef[3][5] = 1.141663041636e-5; + coef[3][6] = -4.340802793033e-6; + coef[3][7] = 3.517971869029e-7; + coef[3][8] = -9.170850253981e-09; + + coef[4][0] = -8.519413589968e-4; + coef[4][1] = -6.343059502294e-4; + coef[4][2] = 1.039643390686e-3; + coef[4][3] = 8.911036876068e-6; + coef[4][4] = -1.008928628425e-4; + coef[4][5] = 1.775681984457e-05; + coef[4][6] = -7.003521917385e-7; + coef[4][7] = -4.928692832866e-8; + coef[4][8] = 3.208853883734e-9; + + coef[5][0] = -1.247583860943e-3; + coef[5][1] = -1.919569450380e-4; + coef[5][2] = -1.553840717902e-4; + coef[5][3] = 3.175388949811e-5; + coef[5][4] = 1.080693990468e-5; + coef[5][5] = -3.149286923815e-6; + coef[5][6] = 2.318308730487e-7; + coef[5][7] = 1.756388998863e-10; + coef[5][8] = -3.952740758950e-10; + + coef[6][0] = 3.014307545716e-4; + coef[6][1] = 4.075019351738e-5; + coef[6][2] = 2.670827249272e-6; + coef[6][3] = -4.515123641755e-6; + coef[6][4] = 5.106059413591e-7; + coef[6][5] = 3.105491554749e-8; + coef[6][6] = -6.030983538280e-9; + coef[6][7] = -1.446756795654e-10; + coef[6][8] = 2.739558475782e-11; + + coef[7][0] = -2.499323170044e-5; + coef[7][1] = -2.850044983009e-6; + coef[7][2] = 7.695300597935e-7; + coef[7][3] = 2.187439283954e-7; + coef[7][4] = -1.299275586093e-7; + coef[7][5] = 2.274394089017e-8; + coef[7][6] = -1.755944926274e-9; + coef[7][7] = 7.143183138281e-11; + coef[7][8] = -1.693040208927e-12; + + coef[8][0] = 6.932627237765e-7; + coef[8][1] = 6.966822400446e-8; + coef[8][2] = -3.783302281524e-8; + coef[8][3] = -2.911233951880e-9; + coef[8][4] = 5.117133050290e-9; + coef[8][5] = -1.130988250912e-9; + coef[8][6] = 1.005189187279e-10; + coef[8][7] = -3.989884105603e-12; + coef[8][8] = 6.388219930167e-14; + + double lnrate_cx{0}; + for (int i = 0; i < 9; i++) { + for (int j = 0; j < 9; j++) { + lnrate_cx += coef[i][j] * Kokkos::pow(Kokkos::log(energy), i) * + Kokkos::pow(Kokkos::log(field_info.ion_temperature), j); + } + } + double sigma_cx = + Kokkos::exp(lnrate_cx) / Kokkos::sqrt(particle_velocity_squared); + return sigma_cx; + } + + KOKKOS_FUNCTION // IF THIS IS CHANGED MAY NEED TO CHANGE max_sigma_cx in the + // collision code + double + analytic_cross_section_CE(double energy) const { + return (0.6937e-14 * (1 - 0.155 * Kokkos::log10(energy)) * + (1 - 0.155 * Kokkos::log10(energy))) / + (1 + 0.1112e-14 * Kokkos::pow(energy, 3.3)); + } + + KOKKOS_FUNCTION + void sample_collision_distance(ParticleInfo &particle_info, + const FieldInfo &field_info) const { + // example of sampling a random number + auto rand_gen = random_pool.get_state(); + double x = rand_gen.drand(0., 1.); + random_pool.free_state(rand_gen); + + double sigma_ion = ionization_cross_section(particle_info, field_info); + double sigma_cx = charge_exchange_cross_section(particle_info, field_info); + + // Generate distance and move particle + // printf("\n$$$ Before moving particle [%d], position: (%f, %f, %f) + //$$$\n", particle_info.particle_index, particle_info.position[0], + // particle_info.position[1], particle_info.position[2]); + double l = -Kokkos::log(x) / + (field_info.electron_density * sigma_ion + + field_info.ion_density * sigma_cx) * + 0.01; // m. n in cm^-3 + particle_info.position[0] += l * particle_info.direction[0]; + particle_info.position[1] += l * particle_info.direction[1]; + particle_info.position[2] += l * particle_info.direction[2]; + // printf("\n$$$ After moving particle [%d], position (%f, %f, %f). + // Direction was (%f, %f, %f) with l = %f $$$\n", + // particle_info.particle_index, particle_info.position[0], + // particle_info.position[1], particle_info.position[2], + // particle_info.direction[0], particle_info.direction[1], + // particle_info.direction[2], l); + double mp{938.27e6 / (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2 + particle_info.alpha = + Kokkos::sqrt(mp / + (2.0 * particle_energy(particle_info.particle_index))) * + 100; // 1/(m/s) + } + + // collision event + KOKKOS_FUNCTION + void collide_particle(ParticleInfo &particle_info, + const FieldInfo &field_info) const { + + // Generate Random Numbers + auto rand_gen = random_pool.get_state(); + double ww = rand_gen.drand(0., 1.); + + double sigma_ion = ionization_cross_section(particle_info, field_info); + double sigma_cx = charge_exchange_cross_section(particle_info, field_info); + + // Compute New Direction and Energy and set particle info + // Compute 3 Maxwellian (Gaussian) distributed velocities (cm/s) + double mp{938.27e6 / (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2 + double max_rate_cx = 1.4e-7; + + bool rejection_test = false; + auto old_mag_v = + Kokkos::sqrt(2 * particle_energy(particle_info.particle_index) / mp); + double vx{}; + double vy{}; + double vz{}; + // Loops this until it passes the rejection test + while (!rejection_test) { + // Generate random numbers for the velocity + double x1 = rand_gen.drand(0., 1.); + double x2 = rand_gen.drand(0., 1.); + double y1 = rand_gen.drand(0., 1.); + double y2 = rand_gen.drand(0., 1.); + + // Sample a new velocity from the thermal distribution (cm/s) + vx = Kokkos::sqrt(field_info.ion_temperature / mp) * + Kokkos::sqrt(-2 * Kokkos::log(x1)) * Kokkos::cos(2 * M_PI * x2); + vy = Kokkos::sqrt(field_info.ion_temperature / mp) * + Kokkos::sqrt(-2 * Kokkos::log(x1)) * Kokkos::sin(2 * M_PI * x2); + vz = Kokkos::sqrt(field_info.ion_temperature / mp) * + Kokkos::sqrt(-2 * Kokkos::log(y1)) * Kokkos::sin(2 * M_PI * y2); + + // Compute the relative velocity (note using this isn't thoroughly tested. + // Can be disabled by changing old_mag_v to 0, then it is just vx, vy, vz + // like I originally used. + auto rel_vx = vx - particle_info.direction[0] * old_mag_v; + auto rel_vy = vy - particle_info.direction[1] * old_mag_v; + auto rel_vz = vz - particle_info.direction[2] * old_mag_v; + auto mag_v2 = rel_vx * rel_vx + rel_vy * rel_vy + rel_vz * rel_vz; + + // Generate random number and compare to sigma/sigma_max + if (rand_gen.drand(0., 1.) < + Kokkos::sqrt(mag_v2) * analytic_cross_section_CE(0.5 * mp * mag_v2) / + max_rate_cx) { + rejection_test = true; + } + } + + // Now that the rejection test has passed, set the new information + double mag_v{Kokkos::sqrt(vx * vx + vy * vy + vz * vz)}; + particle_info.direction[0] = vx / mag_v; + particle_info.direction[1] = vy / mag_v; + particle_info.direction[2] = vz / mag_v; + + particle_energy(particle_info.particle_index) = 0.5 * mp * mag_v * mag_v; + // particle_info.energy_group = + // particle_energy(particle_info.particle_index); //Temporary for debugging + // Adjust Weights + double new_weight = particle_info.weight * + (1 - (field_info.electron_density * sigma_ion) / + (field_info.ion_density * sigma_cx + + field_info.electron_density * sigma_ion)); + + // Russian Roulette + double wc = 0.25; + double ws = 1.0; + + if (new_weight < wc) { + if (ww < (1 - new_weight / ws)) { + new_weight = 0; + } else { + new_weight = ws; + } + } + particle_info.weight = new_weight; + + // This definition may be redundant since it is also in the next location + // function + particle_info.alpha = + Kokkos::sqrt(mp / + (2.0 * particle_energy(particle_info.particle_index))) * + 100; // 1/(m/s) + random_pool.free_state(rand_gen); + } + // To here + random_pool_t random_pool; + DG2CrossSection cross_section; + + Kokkos::View particle_energy; +}; +#endif // PUMITALLYOPENMC_DG2PHYSICS_H diff --git a/src/physics/degas2/degas2-main.cpp b/src/physics/degas2/degas2-main.cpp new file mode 100644 index 0000000..9119861 --- /dev/null +++ b/src/physics/degas2/degas2-main.cpp @@ -0,0 +1,451 @@ +/* + * This is created as a proof of concept for the Degas2 Neutral Transport Code + * Provides parallel search capabilities and tallying to degas2 through this + * PUMI-Tally interface + */ + +#include +#include +#include +#include +#include + +#include "DG2Physics.h" +#include "PumiTallyImpl.h" + +struct InputParameters { + InputParameters(int argc, char *argv[]); + std::string mesh_name; + int num_particles = 0; + int max_iterations = 1000; + pumitally::SourceDistribution source_distribution = + pumitally::SourceDistribution::EQUAL; +}; + +void PrintInitialInfo(const std::string &mesh_name, int num_particles) noexcept; + +void SampleInitialParticleState(const Kokkos::View &energy_array, + const Kokkos::View &direction, + const Kokkos::View &alpha); + +struct Fields { + explicit Fields(const Omega_h::Reals ¢roids); + + /** + * Set the fields as mesh tag + * @param mesh Mesh to tag + */ + void SetAsTag(Omega_h::Mesh &mesh) const; + + Omega_h::Reals electron_temperature; + Omega_h::Reals ion_temperature; + Omega_h::Reals electron_density; + Omega_h::Reals ion_density; + Omega_h::Reals bulk_flow_velocity; +}; + +void Transport(pumitally::PumiTallyImpl &pumi_tally, DG2Physics &physics, + const Omega_h::Read &electron_density, + const Omega_h::Read &ion_density, + const Omega_h::Read &electron_temperature, + const Omega_h::Read &ion_temperature, + const Omega_h::Read &bulk_flow_velocity, + const Kokkos::View &initial_direction, + const Kokkos::View &initial_alpha, int max_iterations); + +// ***************************************************************************// +// *************************** Main Function ******************************** // +// ***************************************************************************// +int main(int argc, char *argv[]) { + Kokkos::initialize(argc, argv); + { + // Read input parameters + InputParameters input_params(argc, argv); + PrintInitialInfo(input_params.mesh_name, input_params.num_particles); + + // Initialize PUMI-Tally and Read Fields + auto pumi_tally = pumitally::PumiTallyImpl( + input_params.mesh_name, input_params.num_particles, argc, argv, + input_params.source_distribution); + auto &mesh = pumi_tally.full_mesh; + const auto centroids = pumitally::GetCentroids(mesh); + Fields fields(centroids); + fields.SetAsTag(mesh); + pumi_tally.is_pumipic_initialized = true; + + Kokkos::View sigma_t; // mat, T, g + Kokkos::View sigma_a; // mat, T, g + Kokkos::View scattering_matrix; // mat, T, g, g + Kokkos::View sigma_s; // mat, T, g + DG2CrossSection cross_section(sigma_t, sigma_a, scattering_matrix, sigma_s); + DG2Physics physics(cross_section, input_params.num_particles); + + Kokkos::View initial_direction("Initial Direction", + input_params.num_particles * 3); + Kokkos::View initial_alpha("Initial Alpha", + input_params.num_particles); + SampleInitialParticleState(physics.particle_energy, initial_direction, + initial_alpha); + + // Move particles + auto electron_density = + mesh.get_array(Omega_h::REGION, "electron_density"); + auto ion_density = + mesh.get_array(Omega_h::REGION, "ion_density"); + auto electron_temperature = + mesh.get_array(Omega_h::REGION, "electron_temperature"); + auto ion_temperature = + mesh.get_array(Omega_h::REGION, "ion_temperature"); + auto bulk_flow_velocity = + mesh.get_array(Omega_h::REGION, "bulk_flow_velocity"); + + Transport(pumi_tally, physics, electron_density, ion_density, + electron_temperature, ion_temperature, bulk_flow_velocity, + initial_direction, initial_alpha, input_params.max_iterations); + + // Finalize and output results + pumi_tally.WriteTallyResults(); + } + Kokkos::finalize(); + + return 0; +} + +// **************************************************************************************************************** +// // +// ********************************************* Function Definitions +// ********************************************* // +// **************************************************************************************************************** +// // + +void Transport(pumitally::PumiTallyImpl &pumi_tally, DG2Physics &physics, + const Omega_h::Read &electron_density, + const Omega_h::Read &ion_density, + const Omega_h::Read &electron_temperature, + const Omega_h::Read &ion_temperature, + const Omega_h::Read &bulk_flow_velocity, + const Kokkos::View &initial_direction, + const Kokkos::View &initial_alpha, + int max_iterations) { + + for (int iter = 0; iter < max_iterations; ++iter) { + auto particle_dest = pumi_tally.pumipic_ptcls->get<1>(); + auto particle_orig = pumi_tally.pumipic_ptcls->get<0>(); + auto particle_weight = pumi_tally.pumipic_ptcls->get<4>(); + auto particle_group = pumi_tally.pumipic_ptcls->get<5>(); + + auto last_exit = + pumi_tally.p_pumi_particle_at_elem_boundary_handler->last_exit_; + + printf("\n### iter = %d ###\n", iter); + + auto alpha = pumi_tally.p_pumi_particle_at_elem_boundary_handler->alpha_; + auto get_new_destination = + PS_LAMBDA(const int &e, const int &pid, const int &mask) { + + if (mask > 0) { // FIXME: check if the particle is at destination or at + // the boundary + // + ParticleInfo particle_info; + particle_info.position[0] = particle_orig(pid, 0); + particle_info.position[1] = particle_orig(pid, 1); + particle_info.position[2] = particle_orig(pid, 2); + + if (iter == 0) { + particle_info.direction[0] = initial_direction(pid * 3 + 0); + particle_info.direction[1] = initial_direction(pid * 3 + 1); + particle_info.direction[2] = initial_direction(pid * 3 + 2); + } else { + auto direction = Omega_h::normalize(Omega_h::Vector<3>{ + particle_info.position[0] - particle_dest(pid, 0), + particle_info.position[1] - particle_dest(pid, 1), + particle_info.position[2] - particle_dest(pid, 2)}); + particle_info.direction[0] = direction[0]; + particle_info.direction[1] = direction[1]; + particle_info.direction[2] = direction[2]; + } + particle_info.weight = particle_weight(pid); + particle_info.energy_group = particle_group(pid); + particle_info.particle_index = pid; + + FieldInfo field_info; + field_info.electron_temperature = electron_temperature[e]; + field_info.ion_temperature = ion_temperature[e]; + field_info.electron_density = electron_density[e]; + field_info.ion_density = ion_density[e]; + field_info.bulk_flow_velocity[0] = bulk_flow_velocity[e * 3 + 0]; + field_info.bulk_flow_velocity[1] = bulk_flow_velocity[e * 3 + 1]; + field_info.bulk_flow_velocity[2] = bulk_flow_velocity[e * 3 + 2]; + + // printf("\n---- Direction (pid = %d) before collision: (%f, %f, + //%f). Original pos: (%f, %f, %f), pos (%f, %f, %f) ---- \n", pid, + // direction[0], direction[1], direction[2], + // particle_orig(pid,0), particle_orig(pid,1), particle_orig(pid,2), + // particle_info.position[0], particle_info.position[1], + // particle_info.position[2]); + + if (iter != 0) { // first iteration, we do not need to collide + if (last_exit[pid] == -1) { // reached destination + physics.collide_particle(particle_info, field_info); + } + } + + physics.sample_collision_distance(particle_info, field_info); + + // Update particle position and direction + particle_dest(pid, 0) = particle_info.position[0]; + particle_dest(pid, 1) = particle_info.position[1]; + particle_dest(pid, 2) = particle_info.position[2]; + particle_weight(pid) = particle_info.weight; + particle_group(pid) = 0; // physics kernel does not handle groups yet + + alpha[pid] = (iter == 0) ? initial_alpha(pid) : particle_info.alpha; + } + }; + pumipic::parallel_for(pumi_tally.pumipic_ptcls.get(), get_new_destination, + "get new destination"); + pumi_tally.SearchAndRebuild( + false, true); // for now, always rebuild the pp structure + Kokkos::fence(); + } +} + +void initialize_equal_source(const Omega_h::Mesh &mesh, + Omega_h::Write particle_positions) { + Omega_h::LO num_particles = particle_positions.size() / 3; + Omega_h::LO particles_per_element = num_particles / mesh.nelems(); + Omega_h::LO remainder_particles = num_particles % mesh.nelems(); + const auto centroids = + mesh.get_array(Omega_h::REGION, "centroid"); + + // set the particle positions to the centroids of each element, + // and adding one particle to the first few elements if there are remainder + // particles + Omega_h::parallel_for( + "set source particles", mesh.nelems(), OMEGA_H_LAMBDA(Omega_h::LO e) { + Omega_h::Vector<3> centroid = { + centroids[e * 3 + 0], centroids[e * 3 + 1], centroids[e * 3 + 2]}; + Omega_h::LO particles_in_this_element = + particles_per_element + (e < remainder_particles ? 1 : 0); + for (Omega_h::LO p = 0; p < particles_in_this_element; ++p) { + Omega_h::LO particle_index = + e * particles_per_element + p + + (e < remainder_particles ? e : remainder_particles); + particle_positions[particle_index * 3 + 0] = centroid[0]; + particle_positions[particle_index * 3 + 1] = centroid[1]; + particle_positions[particle_index * 3 + 2] = centroid[2]; + } + }); + Kokkos::fence(); +} + +/* +void set_source_particles(Omega_h::Mesh &mesh, + const SourceDistribution source_distribution, + Omega_h::Write particle_positions) { + OMEGA_H_CHECK_PRINTF(particle_positions.size() > 0, + "Particle positions size (%ld) must be greater than 0\n", + particle_positions.size()); + Omega_h::LO num_particles = particle_positions.size() / 3; + + switch (source_distribution) { + case SourceDistribution::UNIFORM: { + auto ppe = pumiinopenmc::PPPS::kkLidView("ptcls_per_element", + mesh.nelems()); + distributeParticlesBasesOnVolume(mesh, ppe, num_particles); + initialize_uniform_source(mesh, particle_positions, ppe); + break; + } + case SourceDistribution::EQUAL: { + initialize_equal_source(mesh, particle_positions); + + break; + } + default: + throw std::runtime_error("Unsupported source distribution!\n"); + } +} +*/ + +void Fields::SetAsTag(Omega_h::Mesh &mesh) const { + OMEGA_H_CHECK_PRINTF(electron_temperature.size() == mesh.nelems(), + "Electron temperature size (%d) must be equal to number " + "of elements (%d)\n", + electron_temperature.size(), mesh.nelems()); + OMEGA_H_CHECK_PRINTF( + ion_temperature.size() == mesh.nelems(), + "Ion temperature size (%d) must be equal to number of elements (%d)\n", + ion_temperature.size(), mesh.nelems()); + OMEGA_H_CHECK_PRINTF( + electron_density.size() == mesh.nelems(), + "Electron density size (%d) must be equal to number of elements (%d)\n", + electron_density.size(), mesh.nelems()); + OMEGA_H_CHECK_PRINTF( + ion_density.size() == mesh.nelems(), + "Ion density size (%d) must be equal to number of elements (%d)\n", + ion_density.size(), mesh.nelems()); + OMEGA_H_CHECK_PRINTF(bulk_flow_velocity.size() == mesh.nelems() * 3, + "Bulk flow velocity size (%d) must be equal to number " + "of elements * 3 (%d)\n", + bulk_flow_velocity.size(), mesh.nelems() * 3); + + // Assuming the mesh has a field for each of the properties + mesh.add_tag(Omega_h::REGION, "electron_temperature", 1, + electron_temperature); + mesh.add_tag(Omega_h::REGION, "ion_temperature", 1, + ion_temperature); + mesh.add_tag(Omega_h::REGION, "electron_density", 1, + electron_density); + mesh.add_tag(Omega_h::REGION, "ion_density", 1, ion_density); + mesh.add_tag(Omega_h::REGION, "bulk_flow_velocity", 3, + bulk_flow_velocity); +} + +// TODO: Implement the function to retrieve field values based on centroids +// NOTE: This requires the python script and the plasma source mesh and data +// to be in the folder with this cpp file. The python script may need to be +// updated to read in different plasma source data. +Fields::Fields(const Omega_h::Reals ¢roids) { + auto centroids_h = Omega_h::HostRead(centroids); + + // This function should be implemented to retrieve the field values + // based on the centroids provided (read centroids_h on CPU. + // It contains centroids like e0_0, e0_1, e0_2, e1_0, e1_1, e1_2, ....) + // For now, we will just fill the fields with dummy data. + const Omega_h::LO num_elements = centroids.size() / 3; // Assuming 3D + Omega_h::HostWrite host_electron_temperature(num_elements, "host_e_temp"); + Omega_h::HostWrite host_ion_temperature(num_elements, "host_i_temp"); + Omega_h::HostWrite host_electron_density(num_elements, "host_e_density"); + Omega_h::HostWrite host_ion_density(num_elements, "host_i_density"); + Omega_h::HostWrite host_bulk_flow_velocity(3 * num_elements, "host_bulk_flow_velocity"); + + // Read in the field values + std::ifstream fieldvals{"vals.csv"}; + int i = 0; + std::string id_string; + std::string xx_string; + std::string yy_string; + std::string zz_string; + + while (std::getline(fieldvals, id_string, ',')) { + + std::getline(fieldvals, xx_string, ','); + std::getline(fieldvals, yy_string, ','); + std::getline(fieldvals, zz_string, '\n'); + + // Asign field values + host_ion_density[i] = std::stod(xx_string); + host_electron_density[i] = std::stod(xx_string); + host_electron_temperature[i] = std::stod(yy_string); + host_ion_temperature[i] = std::stod(zz_string); + + i++; + } + + ion_density = Omega_h::Reals(host_ion_density); + electron_density = Omega_h::Reals(host_electron_density); + electron_temperature = Omega_h::Reals(host_electron_temperature); + ion_temperature = Omega_h::Reals(host_ion_temperature); +} + +void PrintInitialInfo(const std::string &mesh_name, + int num_particles) noexcept { + printf("\n===================> Accelerated Degas2 <======================\n"); + std::string mesh_particle_info = + "\tMesh: " + mesh_name + + "\n\tNumber of particles: " + std::to_string(num_particles) + "\n"; + printf("%s", mesh_particle_info.c_str()); + printf("\n===============================================================\n"); +} + +InputParameters::InputParameters(const int argc, char *argv[]) { + if (argc != 5) { + throw std::runtime_error( + "Usage: " + std::string(argv[0]) + + " "); + } + mesh_name = argv[1]; + num_particles = std::stoi(argv[2]); + if (num_particles <= 0) { + throw std::runtime_error("Number of particles must be a positive integer."); + } + max_iterations = std::stoi(argv[3]); + + std::string source_dist_str = argv[4]; + // FIXME: hardcoded lower case + if (source_dist_str == "uniform") { + source_distribution = pumitally::SourceDistribution::UNIFORM; + } else if (source_dist_str == "equal") { + source_distribution = pumitally::SourceDistribution::EQUAL; + } else { + throw std::runtime_error( + "Invalid source distribution. Use 'uniform' or 'equal'."); + } +} + +// todo: split this function to sample uniform direction and sample energy +// and move to impl class +void SampleInitialParticleState(const Kokkos::View &energy_array, + const Kokkos::View &direction, + const Kokkos::View &alpha) { + const random_pool_t random_pool(SEED); + + auto sample_energy = OMEGA_H_LAMBDA(const int i) { + // Basically, feed this 4 uniformly generated random numbers, x1, x2, y1, + // y2, on the interval (0,1) and it + // will give you a velocity vector (vx,vy,vz), a unit direction vector + // (directionx, directiony, directionz), the alpha value used in the tally, + // and the energy of the particle that was sampled. All from a gas at + // temperature temp. + constexpr double kMp{938.27e6 / + (3e10 * 3e10)}; // eV/c^2 = eV*s^2/cm^2. Necessary + // constant for the distribution + + // This is what we set as the source temperature. This could in principle be + // a parameter but for now can be hard coded in as 3 eV or so and I can + // always tweak it. + constexpr double kSourceTemp = 3; // eV. + + // This uses the Box-Muller method of sampling a Gaussian, which generates + // two independent normally distributed values from two independent + // uniformly distributed numbers on the interval (0,1). (x1, x2) is the + // first pair. (y1, y2) is the second pair. Note this technically could make + // 4 independent normally distributed values, but we only need three. This + // generates a velocity vector (vx,vy,vz) sampled from an ideal gas at + // temperature temp. + auto generator = random_pool.get_state(); + const double x1 = generator.drand(0.0, 1.0); + const double x2 = generator.drand(0.0, 1.0); + const double y1 = generator.drand(0.0, 1.0); + const double y2 = generator.drand(0.0, 1.0); + random_pool.free_state(generator); + const double vx = Kokkos::sqrt(kSourceTemp / kMp) * + Kokkos::sqrt(-2 * Kokkos::log(x1)) * + Kokkos::cos(2 * M_PI * x2); + const double vy = Kokkos::sqrt(kSourceTemp / kMp) * + Kokkos::sqrt(-2 * Kokkos::log(x1)) * + Kokkos::sin(2 * M_PI * x2); + const double vz = Kokkos::sqrt(kSourceTemp / kMp) * + Kokkos::sqrt(-2 * Kokkos::log(y1)) * + Kokkos::sin(2 * M_PI * y2); + + // Compute the magnitude of the velocity vector + const double mag_v = Kokkos::sqrt(vx * vx + vy * vy + vz * vz); + + // Compute the alpha factor multiplied in when computing the tally + alpha(i) = 1 / mag_v; + + // Compute the unit normal + direction(i * 3 + 0) = vx / mag_v; + direction(i * 3 + 1) = vy / mag_v; + direction(i * 3 + 2) = vz / mag_v; + + // Compute the particle energy. This you will want to save probably to pass + // to my initialization (set_energy) code so that I have the energy for my + // computations. + const double particle_energy = 0.5 * kMp * mag_v * mag_v; + energy_array(i) = particle_energy; + }; + Kokkos::parallel_for("SampleInitialParticleState", energy_array.size(), + sample_energy); +} diff --git a/src/pumitally/PumiTallyImpl.cpp b/src/pumitally/PumiTallyImpl.cpp index d9b7057..82ca6ee 100644 --- a/src/pumitally/PumiTallyImpl.cpp +++ b/src/pumitally/PumiTallyImpl.cpp @@ -19,6 +19,33 @@ std::unique_ptr CreateParticleDS(const Omega_h::Mesh &mesh, void InitializeParticlesInElement0(Omega_h::Mesh &mesh, pumitally::PPPS *ptcls); +Omega_h::Reals GetCentroids(Omega_h::Mesh &mesh, const bool add_tag) { + const auto coords = mesh.coords(); + const auto nelems = mesh.nelems(); + const auto e2v = mesh.ask_down(Omega_h::REGION, Omega_h::VERT).ab2b; + + const Omega_h::Write centroids(nelems * 3, 0.0, "centroids"); + + // FIXME: Hardcoded for 3D tets + Omega_h::parallel_for( + "calculate centroids", nelems, OMEGA_H_LAMBDA(int e) { + const auto nodes = o::gather_verts<4>(e2v, e); + o::Few, 4> elem_coords = + o::gather_vectors<4, 3>(coords, nodes); + o::Vector<3> centroid = o::average(elem_coords); + + centroids[e * 3 + 0] = centroid[0]; + centroids[e * 3 + 1] = centroid[1]; + centroids[e * 3 + 2] = centroid[2]; + }); + + if (add_tag) { + mesh.add_tag(Omega_h::REGION, "centroid", 3, centroids); + } + + return centroids; +} + void TallyTimes::PrintTimes() const { printf("\n"); printf("[TIME] Initialization time : %f seconds\n", initialization_time); @@ -29,7 +56,8 @@ void TallyTimes::PrintTimes() const { } PumiTallyImpl::PumiTallyImpl(const std::string &mesh_filename, - const Omega_h::LO num_ptcls, int argc, char **argv) + const Omega_h::LO num_ptcls, int argc, char **argv, + const SourceDistribution source_dist) : num_particles(num_ptcls) { oh_mesh_filename = mesh_filename; @@ -43,7 +71,22 @@ PumiTallyImpl::PumiTallyImpl(const std::string &mesh_filename, // todo can track lengths be here? LoadMeshAndInitParticles(argc, argv); - InitializeParticlesInElement0(*p_picparts->mesh(), pumipic_ptcls.get()); + + switch (source_dist) { + case SourceDistribution::UNIFORM: + throw std::runtime_error( + "UNIFORM source distribution is not implemented yet"); + break; + case SourceDistribution::EQUAL: + throw std::runtime_error( + "EQUAL source distribution is not implemented yet"); + break; + case SourceDistribution::ZERO: + InitializeParticlesInElement0(*p_picparts->mesh(), pumipic_ptcls.get()); + break; + default: + throw std::runtime_error("Invalid source distribution"); + } p_particle_tracer = std::make_unique< ParticleTracer>( @@ -253,6 +296,340 @@ void UpdateCurrentElement(PPPS *ptcls, pumipic::parallel_for(ptcls, move_to_next, "move to next element"); } +void compute_boundary_normals(Omega_h::Mesh &mesh) { + const auto exposed_edges = Omega_h::mark_exposed_sides(&mesh); + const auto face2elems = mesh.ask_up(mesh.dim() - 1, mesh.dim()).ab2b; + const auto face2elemsOffset = mesh.ask_up(mesh.dim() - 1, mesh.dim()).a2ab; + const auto elem2nodes = mesh.ask_down(mesh.dim(), 0).ab2b; + const auto face2nodes = mesh.ask_down(mesh.dim() - 1, 0).ab2b; + const auto coords = mesh.coords(); + + Omega_h::Write normals(mesh.nfaces() * 3, 0.0, + "boundary_normals"); + + // calculate the normals for the exposed edges + auto calculate_normals = OMEGA_H_LAMBDA(const Omega_h::LO &face_id) { + if (exposed_edges[face_id]) { + // get this face's nodes + const auto face_nodes = Omega_h::gather_verts<3>(face2nodes, face_id); + const auto face_coords = + Omega_h::gather_vectors<3, 3>(coords, face_nodes); + const auto normal = Omega_h::cross( + face_coords[1] - face_coords[0], + face_coords[2] - face_coords[0]); // cross product to get normal + + const auto norm = Omega_h::norm(normal); + + // fix direction of the normal: get the fourth node of the element and + // check with it + const auto elem_id = + face2elems[face2elemsOffset[face_id]]; // edge sides only have one + // element + const auto elem_nodes = Omega_h::gather_verts<4>(elem2nodes, elem_id); + int fourth_node = -1; + for (int i = 0; i < 4; ++i) { + if (elem_nodes[i] != face_nodes[0] && elem_nodes[i] != face_nodes[1] && + elem_nodes[i] != face_nodes[2]) { + fourth_node = elem_nodes[i]; + break; + } + } + OMEGA_H_CHECK_PRINTF(fourth_node != -1, + "Error: fourth node not found for face %d\n", + face_id); + + const Omega_h::Vector<3> fourth_node_coord = { + coords[fourth_node * 3 + 0], coords[fourth_node * 3 + 1], + coords[fourth_node * 3 + 2]}; + + // check if the normal points towards the fourth node + Omega_h::Vector<3> fourth_2_face_vector = { + fourth_node_coord[0] - face_coords[0][0], + fourth_node_coord[1] - face_coords[0][1], + fourth_node_coord[2] - face_coords[0][2]}; + + Omega_h::Vector<3> inner_norm; + if (Omega_h::inner_product(normal, fourth_2_face_vector) < 0) { + // flip the normal if it points away from the fourth node + for (int i = 0; i < 3; ++i) { + inner_norm[i] = -normal[i] / norm; + } + } else { + for (int i = 0; i < 3; ++i) { + inner_norm[i] = normal[i] / norm; + } + } + // store the normal in the normals array + normals[face_id * 3 + 0] = inner_norm[0]; + normals[face_id * 3 + 1] = inner_norm[1]; + normals[face_id * 3 + 2] = inner_norm[2]; + } + }; + Omega_h::parallel_for(mesh.nfaces(), calculate_normals, + "compute boundary normals"); + mesh.add_tag(Omega_h::FACE, "normals", 3, Omega_h::Reals(normals)); +} + +OMEGA_H_DEVICE o::Real volume_tet(const o::Few, 4> &tet_verts) { + o::Few, 3> basis33 = {tet_verts[1] - tet_verts[0], + tet_verts[2] - tet_verts[0], + tet_verts[3] - tet_verts[0]}; + auto volume = o::tet_volume_from_basis(basis33); + return volume; +} + +o::Real volume_of_3d_mesh(o::Mesh &mesh) { + OMEGA_H_CHECK_PRINTF(mesh.dim() == 3, + "Volume calculation is only supported for 3D meshes, " + "but got %dD mesh\n", + mesh.dim()); + const auto coords = mesh.coords(); + const auto elems2nodes = mesh.ask_down(o::REGION, o::VERT).ab2b; + const auto n_elems = mesh.nelems(); + o::Real total_volume = 0.0; + + Kokkos::parallel_reduce( + n_elems, + KOKKOS_LAMBDA(const int i, o::Real &local_volume) { + auto elem_nodes = o::gather_verts<4>(elems2nodes, i); + o::Few, 4> elem_coords; + elem_coords = o::gather_vectors<4, 3>(coords, elem_nodes); + o::Real elem_volume = volume_tet(elem_coords); + local_volume += elem_volume; + }, + Kokkos::Sum(total_volume)); + + return total_volume; +} + +void apply_reflection_boundary_condition( + Omega_h::Mesh &mesh, PPPS *ptcls, Omega_h::Write &elem_ids, + Omega_h::Write &next_elems, + Omega_h::Write &ptcl_done, + Omega_h::Write &lastExit, Omega_h::Write &xFace, + Omega_h::Write &inter_points, + Omega_h::Write material_ids, bool initial) { + + // TODO: make this a member variable of the struct + auto particle_destination = ptcls->get<1>(); + auto particle_origin = ptcls->get<0>(); + + const auto class_ids = mesh.get_array(3, "class_id"); + const auto normals = mesh.get_array(Omega_h::FACE, "normals"); + + auto checkExposedEdges = + PS_LAMBDA(const int e, const int pid, const int mask) { + if (mask > 0 && !ptcl_done[pid]) { + bool reached_destination = (lastExit[pid] == -1); + bool hit_outer_boundary = + ((next_elems[pid] == -1) && (elem_ids[pid] != -1)); + + // for the initial run, we need to find the initial position of the + // particles + bool hit_material_boundary = false; + if (!initial) { // stop at geometry boundary + if (next_elems[pid] != -1) { + if (class_ids[elem_ids[pid]] != + class_ids[next_elems[pid]]) { // particle crosses geometry + // boundary + hit_material_boundary = true; + material_ids[pid] = class_ids[next_elems[pid]]; + } + } else { + material_ids[pid] = -1; // no material id if not in an element + } + } + + ptcl_done[pid] = + (reached_destination || hit_material_boundary) ? 1 : ptcl_done[pid]; + // assert that if the next element is -1, then the material id is -1 + if (!initial) { + if (next_elems[pid] == -1) { + OMEGA_H_CHECK_PRINTF( + material_ids[pid] == -1, + "Error: next_elems[%d] is -1 but material_ids[%d] " + "is %d\n", + pid, pid, material_ids[pid]); + } + // printf("Pid %d next element %d, elem id %d, material id %d\n", + // pid, next_elems[pid], elem_ids[pid], material_ids[pid]); + } + + // reflective boundary condition + if (hit_outer_boundary) { // just reached the boundary + // printf("Moving particle %4d (from -> to): element (%5d -> %5d) + // Material (%3d -> %3d)\n", + // pid, elem_ids[pid], next_elems[pid], + // class_ids[elem_ids[pid]], material_ids[pid]); + // printf("P %d in element %d hit lastExit %d xFace %d \n", pid, + // elem_ids[pid], lastExit[pid], xFace[pid]); + // xFace[pid] = lastExit[pid]; + Omega_h::LO hit_face = + (lastExit[pid] == -1) ? xFace[pid] : lastExit[pid]; + xFace[pid] = hit_face; + lastExit[pid] = hit_face; + ptcl_done[pid] = 1; // stop the particle after reflection + next_elems[pid] = elem_ids[pid]; // reflects back to the same element + OMEGA_H_CHECK_PRINTF(hit_face != -1, + "Error: xFace[%d] is -1 but " + "hit_outer_boundary is true\n", + pid); + + // change direction + auto normal = Omega_h::Vector<3>{normals[hit_face * 3 + 0], + normals[hit_face * 3 + 1], + normals[hit_face * 3 + 2]}; + Omega_h::Vector<3> incident_vector = { + particle_destination(pid, 0) - particle_origin(pid, 0), + particle_destination(pid, 1) - particle_origin(pid, 1), + particle_destination(pid, 2) - particle_origin(pid, 2)}; + // reflect the particle direction + Omega_h::Vector<3> reflected_vector = + incident_vector - + 2.0 * Omega_h::inner_product(incident_vector, normal) * normal; + + // change the particle's position + // particle reaches the boundary + particle_origin(pid, 0) = inter_points[pid * 3]; + particle_origin(pid, 1) = inter_points[pid * 3 + 1]; + particle_origin(pid, 2) = inter_points[pid * 3 + 2]; + + particle_destination(pid, 0) = + particle_origin(pid, 0) + reflected_vector[0]; + particle_destination(pid, 1) = + particle_origin(pid, 1) + reflected_vector[1]; + particle_destination(pid, 2) = + particle_origin(pid, 2) + reflected_vector[2]; + } + } + }; + pumipic::parallel_for(ptcls, checkExposedEdges, + "apply reflective boundary condition"); +} + +void distributeParticlesBasesOnVolume(Omega_h::Mesh &mesh, + pumitally::PPPS::kkLidView ppe, + const int numPtcls) { + OMEGA_H_CHECK_PRINTF(mesh.dim() == 3, + "Distributing particles based on volume is only " + "supported for 3D meshes, but got %dD mesh\n", + mesh.dim()); + o::LO ne = mesh.nelems(); + o::Real mesh_volume = volume_of_3d_mesh(mesh); + OMEGA_H_CHECK(mesh_volume > 0.0); + + auto coords = mesh.coords(); + auto element2nodes = mesh.ask_down(o::REGION, o::VERT).ab2b; + + auto distribute_based_on_volume = OMEGA_H_LAMBDA(o::LO e) { + auto verts = o::gather_verts<4>(element2nodes, e); + auto vert_coords = o::gather_vectors<4, 3>(coords, verts); + o::Real vol = volume_tet(vert_coords); +#ifdef DEBUG + OMEGA_H_CHECK(area > 0.0); +#endif + o::Real volume_fraction = vol / mesh_volume; + ppe[e] = std::round(numPtcls * volume_fraction); + }; + o::parallel_for(ne, distribute_based_on_volume); + + Omega_h::LO totPtcls = 0; + Kokkos::parallel_reduce( + ppe.size(), + OMEGA_H_LAMBDA(const int i, Omega_h::LO &lsum) { lsum += ppe[i]; }, + totPtcls); + + // remove or add particles to match the total number of particles + int extra_particles = numPtcls - totPtcls; + // go throught the first extra_particles elements and add/remove one particle + OMEGA_H_CHECK_PRINTF(extra_particles <= mesh.nelems(), + "Extra particles (%d) should be less than or equal to " + "number of elements (%d)\n", + extra_particles, mesh.nelems()); + + int add_remove = (extra_particles > 0) ? 1 : -1; + auto add_or_remove_particles = OMEGA_H_LAMBDA(o::LO e) { + ppe[e] += add_remove; + }; + o::parallel_for(std::abs(extra_particles), add_or_remove_particles); +} + +OMEGA_H_DEVICE o::Vector<3> +barycentric2real(const o::Few, 4> &tet_verts, + const o::Vector<4> &bary) { + o::Vector<3> real_coords{0, 0, 0}; + for (int i = 0; i < 4; ++i) { + real_coords += bary[i] * tet_verts[i]; + } + return real_coords; +} + +void initialize_uniform_source(Omega_h::Mesh &mesh, + Omega_h::Write particle_positions, + pumitally::PPPS::kkLidView ppe) { + int dim = mesh.dim(); + OMEGA_H_CHECK(dim == 3); + + // cumulative sum of particles per element + Omega_h::Write cumulative_particles(mesh.nelems() + 1, 0); + Omega_h::LO num_particles_cumsum = 0; + auto calculate_cumulative_number_of_particles = + KOKKOS_LAMBDA(const int &e, Omega_h::LO &cumulative, bool is_final) { + auto num_particles = ppe[e]; + cumulative += num_particles; + if (is_final) { + cumulative_particles[e + 1] = cumulative; + } + }; + Kokkos::parallel_scan("calculate_cumulative_number_of_particles", + mesh.nelems(), calculate_cumulative_number_of_particles, + num_particles_cumsum); + + OMEGA_H_CHECK_PRINTF(num_particles_cumsum == particle_positions.size() / 3, + "Total number of particles (%ld) does not match " + "cumulative particles (%ld)\n", + particle_positions.size() / 3, num_particles_cumsum); + + const auto cells2nodes = mesh.ask_down(o::REGION, o::VERT).ab2b; + const auto coords = mesh.coords(); + + Kokkos::Random_XorShift64_Pool random_pool(0); + auto set_initial_positions = OMEGA_H_LAMBDA(const int &e) { + auto pid_start = cumulative_particles[e]; + auto pid_end = cumulative_particles[e + 1]; + auto num_particles_in_element = pid_end - pid_start; + OMEGA_H_CHECK(num_particles_in_element >= 0); + + for (Omega_h::LO pid = pid_start; pid < pid_end; ++pid) { + auto gen = random_pool.get_state(); + o::Real r1 = gen.drand(0.0, 1.0); + o::Real r2 = gen.drand(0.0, 1.0); + o::Real r3 = gen.drand(0.0, 1.0); + + r1 = Kokkos::pow(r1, 1.0 / 3.0); + r2 = Kokkos::sqrt(r2); + o::Real a = 1.0 - r1; + o::Real b = r1 * (1.0 - r2); + o::Real c = r1 * r2 * (1.0 - r3); + o::Real d = r1 * r2 * r3; + + o::Vector<4> random_bcc{a, b, c, d}; + + random_pool.free_state(gen); + + auto verts = o::gather_verts<4>(cells2nodes, e); + auto vert_coords = o::gather_vectors<4, 3>(coords, verts); + + auto real_loc = barycentric2real(vert_coords, random_bcc); + + particle_positions[pid * 3 + 0] = real_loc[0]; + particle_positions[pid * 3 + 1] = real_loc[1]; + particle_positions[pid * 3 + 2] = real_loc[2]; + } + }; + o::parallel_for(mesh.nelems(), set_initial_positions); +} + void ApplyVacuumBC(const Omega_h::Mesh &mesh, PPPS *ptcls, const Omega_h::Write &elem_ids, const Omega_h::Write &next_elems, diff --git a/src/pumitally/PumiTallyImpl.h b/src/pumitally/PumiTallyImpl.h index 9740e08..59c8ba0 100644 --- a/src/pumitally/PumiTallyImpl.h +++ b/src/pumitally/PumiTallyImpl.h @@ -12,6 +12,8 @@ namespace pumitally { +Omega_h::Reals GetCentroids(Omega_h::Mesh &mesh, bool add_tag = true); + /** * Data structure to hold the timing information for different sections */ @@ -26,6 +28,12 @@ struct TallyTimes { void PrintTimes() const; }; +enum class SourceDistribution { + UNIFORM, // Source uniformly distributed across the mesh + EQUAL, // Source at centroids of each element + ZERO // in the zeroth element centroid +}; + /** * @brief PUMI-PiC Data structure Template * @details @@ -35,10 +43,11 @@ struct TallyTimes { * @n 2-ID, * @n 3-in_advance_particle_queue, * @n 4-weight + * @n 5-group */ using PPParticle = pumipic::MemberTypes; + Omega_h::I16, Omega_h::Real, Omega_h::I16>; using PPPS = pumipic::ParticleStructure; //!< PUMI-PiC Particle DS using PPExeSpace = Kokkos::DefaultExecutionSpace; //!< PUMI-PiC Default Execution Space @@ -142,6 +151,11 @@ struct ParticleAtElemBoundary { bool is_initial_track; //!< in is_initial_track run, flux is not tallied Omega_h::Write flux; //!< Flux tally array Omega_h::Write prev_xpoint; //!< Previous intersection point + + // temporary gabe merging variables + // these will be removed after the operator functinality is merged to both + Omega_h::Write last_exit_; + Omega_h::Write alpha_; }; /** @@ -186,7 +200,8 @@ struct PumiTallyImpl { TallyTimes tally_times; //!< Struct to hold times for different operations PumiTallyImpl(const std::string &mesh_filename, Omega_h::LO num_ptcls, - int argc, char **argv); + int argc, char **argv, + SourceDistribution source_dist = SourceDistribution::ZERO); ~PumiTallyImpl() = default; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 40e3261..36bfc2a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -29,4 +29,11 @@ target_include_directories(test_pumi_tally_impl_methods PRIVATE ${PROJECT_SOURCE # PumiPic needed for this since this test case uses Omega_h and PumiPic and pumitally keeps them private target_link_libraries(test_pumi_tally_impl_methods PUBLIC Catch2::Catch2WithMain pumipic::pumipic pumitally) catch_discover_tests(test_pumi_tally_impl_methods) + +if(PUMITALLY_ENABLE_DEGAS2) + add_executable(test_degas2_physics test_degas2_physics.cpp) + target_include_directories(test_degas2_physics PRIVATE ${PROJECT_SOURCE_DIR}/src/physics/degas2) + target_link_libraries(test_degas2_physics PUBLIC Catch2::Catch2WithMain pumitally_degas2 pumipic::pumipic) + catch_discover_tests(test_degas2_physics DISCOVERY_MODE PRE_TEST) +endif() #################### Test HostWrite #################### diff --git a/test/assets/tet6-222.osh/0.osh b/test/assets/tet6-222.osh/0.osh new file mode 100644 index 0000000..085bb07 Binary files /dev/null and b/test/assets/tet6-222.osh/0.osh differ diff --git a/test/assets/tet6-222.osh/nparts b/test/assets/tet6-222.osh/nparts new file mode 100644 index 0000000..d00491f --- /dev/null +++ b/test/assets/tet6-222.osh/nparts @@ -0,0 +1 @@ +1 diff --git a/test/assets/tet6-222.osh/version b/test/assets/tet6-222.osh/version new file mode 100644 index 0000000..f599e28 --- /dev/null +++ b/test/assets/tet6-222.osh/version @@ -0,0 +1 @@ +10 diff --git a/test/test_degas2_physics.cpp b/test/test_degas2_physics.cpp new file mode 100644 index 0000000..640a6fd --- /dev/null +++ b/test/test_degas2_physics.cpp @@ -0,0 +1,136 @@ +// +// Created by Fuad Hasan on 6/20/25. +// + +#include +#include + +#include "DG2Physics.h" +#include +#include +#include +#include + +TEST_CASE("Test Degas2 Physics Functions") { + Kokkos::initialize(); + + { + int numParticles = 1000; + double energy = 3.0; + + Kokkos::View sigma_t_; // mat, T, g + Kokkos::View sigma_a_; // mat, T, g + Kokkos::View scattering_matrix_; // mat, T, g, g + Kokkos::View sigma_s_; // mat, T, g + DG2CrossSection crossSection(sigma_t_, sigma_a_, scattering_matrix_, + sigma_s_); + DG2Physics physics(crossSection, numParticles); + + Kokkos::View particles("particles", numParticles); + Kokkos::View fields("fields", numParticles); + Kokkos::parallel_for( + "initialize_particles", numParticles, KOKKOS_LAMBDA(int i) { + particles(i).position[0] = 0; + particles(i).position[1] = 0; + particles(i).position[2] = 0; + particles(i).direction[0] = 1.0; + particles(i).direction[1] = 0; + particles(i).direction[2] = 0; + particles(i).energy_group = energy; + particles(i).weight = 1.0; + particles(i).particle_index = i; + particles(i).alpha = 1.0; + + fields(i).electron_density = 1.0e13; + fields(i).ion_density = 1.0e13; + fields(i).electron_temperature = 1000.0; + fields(i).ion_temperature = 1000.0; + // and others like this + }); + + // Set the particle energy view + Kokkos::parallel_for( + "Set energies", numParticles, + KOKKOS_LAMBDA(int i) { physics.set_energy(particles(i), energy); }); + + // Call the function to be tested + Kokkos::parallel_for( + "run physics", numParticles, KOKKOS_LAMBDA(int i) { + physics.sample_collision_distance(particles(i), fields(i)); + physics.collide_particle(particles(i), fields(i)); + }); + auto output = create_mirror_view(particles); + Kokkos::deep_copy(output, particles); + std::ofstream outfile("Log.txt"); + + // These 4 functions should return the commented value for SEED=12345 + + double l{0}; + for (int i = 0; i < numParticles; ++i) { + l += output(i).position[0]; + } + l /= numParticles; + outfile << "Average Distance (cm): " << l << std::endl; // 2.41 + + double varl{0}; + for (int i = 0; i < numParticles; ++i) { + varl += (output(i).position[0] - l) * (output(i).position[0] - l); + } + varl /= (numParticles - 1); + double sdl{sqrt(varl)}; + outfile << "Standard Deviation of Distance (cm): " << sdl + << std::endl; // 2.38 + + double ux{0}; + for (int i = 0; i < numParticles; ++i) { + ux += output(i).direction[0]; + } + ux /= numParticles; + outfile << "Mean x Direction: " << ux << std::endl; // 0.0159 + + double varux{0.0}; + for (int i = 0; i < numParticles; ++i) { + varux += (output(i).direction[0] - ux) * (output(i).direction[0] - ux); + } + varux /= (numParticles - 1); + double sdux{sqrt(varux)}; + + outfile << "Standard Deviation of Mean x Direction: " << sdux + << std::endl; // 0.579 + + outfile << std::left << std::setw(10) << "Particle #" << std::setw(16) + << "Velocity(m/s)" << std::setw(10) << "Weight" << std::setw(10) + << "X(cm)" << std::setw(10) << "Y(cm)" << std::setw(10) << "Z(cm)" + << std::setw(10) << "X_Dir" << std::setw(10) << "Y_Dir" + << std::setw(10) << "Z_Dir" << '\n'; + + for (int i = 0; i < numParticles; ++i) { + outfile << std::left << std::setw(10) << i << std::setw(16) + << (1.0 / output(i).alpha) << std::setw(10) << output(i).weight + << std::setw(10) << output(i).position[0] << std::setw(10) + << output(i).position[1] << std::setw(10) << output(i).position[2] + << std::setw(10) << output(i).direction[0] << std::setw(10) + << output(i).direction[1] << std::setw(10) + << output(i).direction[2] << '\n'; + } + + // print the file + std::ifstream infile("Log.txt"); + std::string line; + while (std::getline(infile, line)) { + std::cout << line << '\n'; + } + outfile.close(); + + // Fixme: Where did this tolerance come from. Looks like it's + // not enough for 1000 particles. + // Tolerance is 0.0015 for l, with 95% confidence interval for mean + // with 1000 samples. But we should use something like 2-3 standard + // deviations with more than 99% confidence. + REQUIRE_THAT(l, Catch::Matchers::WithinAbs(0.0240, 0.0015)); + REQUIRE_THAT(sdl, Catch::Matchers::WithinAbs(0.0240, 0.0015)); + REQUIRE_THAT(ux, Catch::Matchers::WithinAbs(0.0, .03)); + REQUIRE_THAT(sdux, Catch::Matchers::WithinAbs(0.577, .03)); + } + Kokkos::finalize(); +}