diff --git a/CHANGELOG.md b/CHANGELOG.md index 6a6d708..8cc48f4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## v0.7.10 +* Documentation improvements +* Improved test coverage +* `convect_particles_and_compute_cell!` function added that computes new cell indices during the convection, a new version of `sort_particles!` +that takes advantage of these pre-computed cell indices has also been added +* `squash_pia!` now also squashes particle cell indices + ## v0.7.9 * `restore_particle_ordering!` added, this restores optimal indexing of particles and can lead to simulation speed-ups due to improved cache usage * Minor optimizations in octree merging diff --git a/PAPER_REPRODUCIBILITY.md b/PAPER_REPRODUCIBILITY.md index 1296a5f..b0afe97 100644 --- a/PAPER_REPRODUCIBILITY.md +++ b/PAPER_REPRODUCIBILITY.md @@ -32,7 +32,7 @@ For the BKW test case (see the commented-out part on "multiple runs with ensembl * `simulations/0D/BKW/bkw_varweight_nnls.jl` - for the variable-weight simulations using NNLS merging For the ionization test case (**when run over all ensembles, these produce VERY LARGE amounts of data, 100s of GBs**); -external data from [LXCat](https://us.lxcat.net/home/) is required for the cross-sections (IST-Lisbon database) (see notes above on exact XML format required): +external data from [LXCat](https://us.lxcat.net/home/) is required for the cross-sections (IST-Lisbon database) (see notes below on exact XML format required): * `simulations/0D/ionization/0D_ionization_1neutralspecies_es.jl` - for the variable-weight simulations using octree N:2 merging (uncomment lines below comment "Uncomment set-up below ..." to get the full set-up running over all parameters and ensembles; setting `paramset` to [12000, 6000] will produce results used as reference values) * `simulations/0D/ionization/0D_ionization_1neutralspecies_nnls_es.jl` - for the variable-weight simulations using different versions of NNLS merging diff --git a/Project.toml b/Project.toml index 89eefed..39eb645 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Merzbild" uuid = "de01fcb4-c117-45a4-a951-7e39b0f12516" authors = ["Georgii Oblapenko ", "Leo Basov "] -version = "0.7.9" +version = "0.7.10" [deps] ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" diff --git a/data/particles.toml b/data/particles.toml index 6876d3a..f42ef4d 100644 --- a/data/particles.toml +++ b/data/particles.toml @@ -16,4 +16,12 @@ charge = 1.0 ["e-"] mass = 9.109383632e-31 -charge = -1.0 \ No newline at end of file +charge = -1.0 + +["Xe"] +mass = 2.181e-25 +charge = 0.0 + +["Xe+"] +mass = 2.181e-25 +charge = 1.0 \ No newline at end of file diff --git a/data/test_neutral_electron_data.xml b/data/test_neutral_electron_data.xml index 9f44d8b..1e6bcc3 100644 --- a/data/test_neutral_electron_data.xml +++ b/data/test_neutral_electron_data.xml @@ -103,7 +103,7 @@ - + knstmrd.github.io @@ -194,4 +194,95 @@ + + + knstmrd.github.io + + + These data are constant (sigma=1e-19 for elastic cross-section, sigma=0.0 for ionization cross-section). + + + + email georgii.oblapenko@fastmail.com + + + + Merzbild.jl reference + + + https://github.com/merzbild/Merzbild.jl + + + + + + Const E-sigma relationship + + + + eHe + EHe + + + + E + Xe -> E + Xe + + + + + 1.37396e-5 + + + complete set + + + + Constant, sigma = 1e-19 + + + 0.000000e+0 1.0e2 5.0e2 1.0e3 + + + 1.0e-19 1.0e-19 1.0e-19 1.0e-19 + + + + + + eAr + EEAr^+ + + + + E + Xe -> E + E + Xe^+ + + + + + 1.212984e+1 + + + complete set + + + + Constant + + + 2024-10-17 + + + + 1.212984e+1 3.0e+1 1.0e+2 + + + + 0.000000e+0 0.0e-20 0.0e-20 + + + + + + + diff --git a/docs/src/modelling_ionization.md b/docs/src/modelling_ionization.md index 7c899b7..b1573ed 100644 --- a/docs/src/modelling_ionization.md +++ b/docs/src/modelling_ionization.md @@ -16,7 +16,12 @@ Examples of synthetic data used for testing can be found in `data/test_neutral_e In general, the XML file has to have the following structure ```xml + + + source. Generated on ---. All rights reserved. + - ??? database, www.lxcat.net, retrieved on --- + @@ -96,8 +101,8 @@ pair; each one is species-specific, therefore a vector for all the neutral speci [`create_computed_crosssections`](@ref) and passing an [`ElectronNeutralInteractions`](@ref) instance. ## Performing ionizing collisions -Electron-neutral collisions with elastic scattering and ionization reactions are modelled by the [`ntc_n_e!`] and -[`ntc_n_e_es!`] functions. The latter implements the Event Splitting +Electron-neutral collisions with elastic scattering and ionization reactions are modelled by the [`ntc_n_e!`](@ref) and +[`ntc_n_e_es!`](@ref) functions. The latter implements the Event Splitting of [Oblapenko et al. (2022)](https://doi.org/10.1016/j.jcp.2022.111390) and reduces the level of stochastic noise, but is suitable only for variable-weight simulations with particle merging. The function implements the collision mechanics as described in [Nanbu (2000)](https://doi.org/10.1109/27.887765); diff --git a/docs/src/overview_1d.md b/docs/src/overview_1d.md index 81f5775..39d29a4 100644 --- a/docs/src/overview_1d.md +++ b/docs/src/overview_1d.md @@ -122,7 +122,7 @@ end ## Performing convection Having set up the grid and boundary conditions, we can convect particles. -This is done by calling the `convect_particles!` function. +This is done by calling the [`convect_particles!`](@ref) function. The convection should be followed by particle sorting before any computations of physical properties are done. ```julia @@ -135,6 +135,20 @@ a `SurfProps` instance needs to be passed: convect_particles!(rng, grid, boundaries, particles[species_id], pia, species_id, species_data, surf_props, Δt) ``` +## Convection and sorting with precomputed particle/cell indices +The approach described above assumes that during the convection process, the indices of the cells the particles find themselves in +are not computed; therefore the call to `sort_particles!` requires passing in the `grid` instance, and the `sort_particles!` +calls a `get_cell` function internally. For 1-D uniform grids, this is an efficient operation, but for other grid types, +the computation of the cell index given only the particle position can be more expensive than keeping track of the cell index inside the convection routine. + +To this purpose, one can call [`convect_particles_and_compute_cell!`](@ref), which will also set the values of the `cell` array of the `ParticleVector` +instance. One can then call another version of `sort_particles!` that does **not** take the grid as a parameter and instead uses the pre-computed `cell` values +to sort the particles: + +``` +sort_particles!(gridsorter, particles[species_id], pia, species_id) +``` + ## Bringing it all together Now we can combine all the pieces to set up a simulation of a single-species Couette flow in a channel with a width of 0.5 mm, discretized with 50 cells. The y-velocity of the left wall is assumed to be -500 m/s, diff --git a/docs/src/reference_public.md b/docs/src/reference_public.md index 27eaede..de5a136 100644 --- a/docs/src/reference_public.md +++ b/docs/src/reference_public.md @@ -172,6 +172,7 @@ sort_particles! ## Particle movement ```@docs convect_particles! +convect_particles_and_compute_cell! ``` ## Particle-surface interactions diff --git a/src/Merzbild.jl b/src/Merzbild.jl index 05e6ff5..e9820f7 100644 --- a/src/Merzbild.jl +++ b/src/Merzbild.jl @@ -61,7 +61,7 @@ export estimate_sigma_g_w_max_ntc_n_e!, ntc_n_e!, ntc_n_e_es! export ParticleVector export AbstractGrid, Grid1DUniform, write_grid export GridSortInPlace, sort_particles! -export MaxwellWallBC, MaxwellWalls1D, convect_particles! +export MaxwellWallBC, MaxwellWalls1D, convect_particles!, convect_particles_and_compute_cell! export pretty_print_pia export ChunkExchanger, exchange_particles!, reset!, sort_particles_after_exchange! export count_disordered_particles, check_pia_is_correct, check_unique_index, check_unique_buffer diff --git a/src/convection/convection_1D.jl b/src/convection/convection_1D.jl index bd27f4f..8346fba 100644 --- a/src/convection/convection_1D.jl +++ b/src/convection/convection_1D.jl @@ -205,4 +205,105 @@ function convect_particles!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D surface_props_scale!(species, surf_props, species_data, Δt) end + +""" + convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, Δt) + +Convect particles on a 1-D uniform grid and write post-convection cell index to `particles.cell`. + +# Positional arguments +* `rng`: the random number generator +* `grid`: the grid on which the convection is performed +* `boundaries`: the `MaxwellWalls1D` struct describing the boundaries (it is assumed that the wall with index 1 is the left wall and + the wall with index 2 is the right wall) +* `particles`: the `ParticleVector` of particles to be convected +* `pia`: the `ParticleIndexerArray` instance +* `species`: the index of the species being convected +* `species_data`: the vector of `Species` data +* `Δt`: the convection timestep +""" +function convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, Δt) + # @inbounds @simd for i in 1:pia.n_total[species] + + @inbounds if pia.contiguous[species] + @inbounds n_tot = pia.n_total[species] + @inbounds @simd for i in 1:n_tot + convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + else + for cell in 1:grid.n_cells + @inbounds s = pia.indexer[cell, species].start1 + @inbounds e = pia.indexer[cell, species].end1 + + @inbounds @simd for i in s:e + convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + + @inbounds if pia.indexer[cell, species].n_group2 > 0 + @inbounds s = pia.indexer[cell, species].start2 + @inbounds e = pia.indexer[cell, species].end2 + + @inbounds @simd for i in s:e + convect_single_particle!(rng, grid, boundaries, particles[i], species, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + end + end + end +end + +""" + convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, surf_props::SurfProps, particles, pia, species, species_data, Δt) + +Convect particles on a 1-D uniform grid and write post-convection cell index to `particles.cell`, computing surface properties if particles hit a surface. + +# Positional arguments +* `rng`: the random number generator +* `grid`: the grid on which the convection is performed +* `boundaries`: the `MaxwellWalls1D` struct describing the boundaries (it is assumed that the wall with index 1 is the left wall and + the wall with index 2 is the right wall) +* `particles`: the `ParticleVector` of particles to be convected +* `pia`: the `ParticleIndexerArray` instance +* `species`: the index of the species being convected +* `species_data`: the vector of `Species` data +* `surf_props`: the `SurfProps` struct where the computed surface properties will be stored +* `Δt`: the convection timestep +""" +function convect_particles_and_compute_cell!(rng, grid::Grid1DUniform, boundaries::MaxwellWalls1D, particles, pia, species, species_data, surf_props::SurfProps, Δt) + # @inbounds @simd for i in 1:pia.n_total[species] + + clear_props!(surf_props) + @inbounds if pia.contiguous[species] + @inbounds n_tot = pia.n_total[species] + @inbounds @simd for i in 1:n_tot + convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + else + for cell in 1:grid.n_cells + @inbounds s = pia.indexer[cell, species].start1 + @inbounds e = pia.indexer[cell, species].end1 + + @inbounds @simd for i in s:e + convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + + @inbounds if pia.indexer[cell, species].n_group2 > 0 + @inbounds s = pia.indexer[cell, species].start2 + @inbounds e = pia.indexer[cell, species].end2 + + @inbounds @simd for i in s:e + convect_single_particle!(rng, grid, boundaries, particles[i], species, surf_props, species_data[species].mass, Δt) + particles.cell[i] = get_cell(grid, particles[i].x) + end + end + end + end + + surface_props_scale!(species, surf_props, species_data, Δt) +end + end \ No newline at end of file diff --git a/src/grids/grid_sorting.jl b/src/grids/grid_sorting.jl index 13b0a22..b46f934 100644 --- a/src/grids/grid_sorting.jl +++ b/src/grids/grid_sorting.jl @@ -44,7 +44,9 @@ GridSortInPlace(grid::G, n_particles::Integer) where {G<:AbstractGrid} = GridSor sort_particles!(gridsort::GridSortInPlace, grid, particles, pia, species) Sort particles on a grid using an in-place sorting algorithm. The `pia` instance is allowed to -have non-contiguous indices (arising for example from merging). +have non-contiguous indices (arising for example from merging). This function +assumes that at the start of the sorting, it is **not known** in which cell each particle is located, +and therefore the cell for each particle has to be determined (by calling `get_cell`). # Positional arguments * `gridsort`: the `GridSortInPlace` structure @@ -107,5 +109,74 @@ function sort_particles!(gridsort::GridSortInPlace, grid, particles, pia, specie @inbounds particles.index[i] = gridsort.sorted_indices[i] end + @inbounds pia.contiguous[species] = true +end + +""" + sort_particles!(gridsort::GridSortInPlace, particles, pia, species) + +Sort particles on a grid using an in-place sorting algorithm. The `pia` instance is allowed to +have non-contiguous indices (arising for example from merging). This function +assumes that at the start of the sorting, it is **known** in which cell each particle is located. + +# Positional arguments +* `gridsort`: the `GridSortInPlace` structure +* `particles`: the `ParticleVector` of particles to be sorted +* `pia`: the `ParticleIndexerArray` instance +* `species`: the index of the species being sorted +""" +function sort_particles!(gridsort::GridSortInPlace, particles, pia, species) + @inbounds n_cells = size(pia.indexer)[1] + @inbounds n_tot = pia.n_total[species] + @inbounds if n_tot > length(gridsort.sorted_indices) + resize!(gridsort.sorted_indices, n_tot + DELTA_PARTICLES) + end + + fill!(gridsort.cell_counts, 0) + + @inbounds if !pia.contiguous[species] + squash_pia!(particles, pia, species) + end + + @inbounds @simd for i in 1:n_tot + gridsort.cell_counts[particles.cell[i]+1] += 1 + end + + @inbounds for cell in 1:n_cells + gridsort.cell_counts[cell+1] = gridsort.cell_counts[cell+1] + gridsort.cell_counts[cell] + + cell_start = gridsort.cell_counts[cell] + 1 + cell_np = gridsort.cell_counts[cell+1] - gridsort.cell_counts[cell] + cell_end = gridsort.cell_counts[cell+1] + + indexer = pia.indexer[cell,species] + + indexer.start2 = 0 + indexer.end2 = -1 + indexer.n_group2 = 0 + + if cell_np > 0 + indexer.start1 = cell_start + indexer.end1 = cell_end + else + # this is done so that we can safely write for i in e1:s1 without worrying about accessing particles at index 0 + indexer.start1 = 0 + indexer.end1 = -1 + end + indexer.n_group1 = cell_np + indexer.n_local = cell_np + end + + @inbounds for i in n_tot:-1:1 + curr_cell = particles.cell[i] + gridsort.sorted_indices[gridsort.cell_counts[curr_cell+1]] = particles.index[i] + + gridsort.cell_counts[curr_cell+1] -= 1 + end + + @inbounds @simd for i in 1:n_tot + @inbounds particles.index[i] = gridsort.sorted_indices[i] + end + @inbounds pia.contiguous[species] = true end \ No newline at end of file diff --git a/src/particles.jl b/src/particles.jl index 6b0c41d..7f5edcc 100644 --- a/src/particles.jl +++ b/src/particles.jl @@ -637,6 +637,7 @@ function squash_pia!(pv, pia, species) @inbounds e2 = pia.indexer[1, species].end2 @inbounds for j in s2:e2 pv.index[j] = pv.index[j+offset] + pv.cell[j] = pv.cell[j+offset] end end end @@ -652,6 +653,7 @@ function squash_pia!(pv, pia, species) @inbounds e1 = pia.indexer[i+1, species].end1 @inbounds for j in s1:e1 pv.index[j] = pv.index[j+offset] + pv.cell[j] = pv.cell[j+offset] end end @inbounds last_end = pia.indexer[i+1, species].end1 > 0 ? pia.indexer[i+1, species].end1 : last_end @@ -668,6 +670,7 @@ function squash_pia!(pv, pia, species) @inbounds e2 = pia.indexer[i, species].end2 @inbounds for j in s2:e2 pv.index[j] = pv.index[j+offset] + pv.cell[j] = pv.cell[j+offset] end end @inbounds last_end = pia.indexer[i, species].end2 > 0 ? pia.indexer[i, species].end2 : last_end @@ -740,7 +743,7 @@ end """ - pretty_print_pia(pia) + pretty_print_pia(pia, species) Display a `ParticleIndexerArray` instance by showing the starting/ending indices of the groups over all cells for a specific species. diff --git a/test/data/ionization_Ar_no_es_v2_invratio.nc b/test/data/ionization_Ar_no_es_v2_invratio.nc new file mode 100644 index 0000000..c671032 Binary files /dev/null and b/test/data/ionization_Ar_no_es_v2_invratio.nc differ diff --git a/test/runtests.jl b/test/runtests.jl index 6ebc338..cfae30b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,7 @@ include("test_nnls_exact_ratepreserving_merging.jl") # test NNLS merging with e include("test_acceleration.jl") # test acceleration of charged particles include("test_collision_ionizing.jl") # test ionizing collisions include("test_ionization_sim.jl") # test 0D ionization simulation without event splitting +include("test_ionization_sim_indexing.jl") # test 0D ionization simulation with multiple neutral/ion species include("test_indexing_particlevector.jl") # test particle vector and multi-cell indexing include("test_indexing_particlebuffer.jl") # test buffer in ParticleVector include("test_pia_contiguous.jl") # test squash_pia diff --git a/test/test_convection_1D.jl b/test/test_convection_1D.jl index 59d2ba0..36e7aa9 100644 --- a/test/test_convection_1D.jl +++ b/test/test_convection_1D.jl @@ -1,4 +1,4 @@ -@testset "convection 1D" begin +@testset "convection 1D and boundaries" begin particles_data_path = joinpath(@__DIR__, "..", "data", "particles.toml") species_data = load_species_data(particles_data_path, "Ar") @@ -53,6 +53,7 @@ @test maximum(abs.(particles[1][4].x - [29.0, 6.0, -3.0])) < 2 * eps() @test particles[1][4].v == [-11.0, -3.0, 1.0] @test particles[1][4].w == 2.0 + phys_props = PhysProps(grid.n_cells, 1, [], Tref=1) compute_props!(particles, pia, species_data, phys_props) @@ -267,4 +268,113 @@ @test abs(phys_props.v[1, 4, 1] - 1.0) < eps() + # finally, test convection where particles cells are computed within the convection routine + pia = ParticleIndexerArray(grid.n_cells, 1) + pia.n_total[1] = 4 + pia.indexer[1,1].n_local = 4 + pia.indexer[1,1].n_group1 = 4 + pia.indexer[1,1].start1 = 1 + pia.indexer[1,1].end1 = 4 + + particles = [ParticleVector(4)] + + # will just move: new x_coord = 20.5 + particles[1][1] = Particle(1.0, [-1.25, -1.5, 4.0], [23.0, -8.0, 7.5]) + + # will reflect from right wall: new x_cord = 50.0 - (10.0 + 11.0) = 29.0 + particles[1][2] = Particle(2.0, [11.0, -3.0, 1.0], [49.0, 6.0, -3.0]) + + # will reflect from left wall: new x_cord = 3.0 + 20.0 = 23.0 + particles[1][3] = Particle(3.0, [-20.0, 0.0, 2.0], [17.0, 1.0, 3.0]) + + # will reflect from left wall and from right wall: new x_coord = 3.55 + particles[1][4] = Particle(4.0, [-49.0, -20.0, 13.0], [1.55, -1.0, 9.0]) + + particles[1].cell[1:4] = [1,1,1,1] + + # 1D specularly reflecting boundaries + boundaries = MaxwellWalls1D(species_data, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0) + + convect_particles_and_compute_cell!(rng, grid, boundaries, particles[1], pia, 1, species_data, 2.0) + + # cells computed correctly in routine + @test particles[1].cell[1:4] == [42, 59, 47, 8] + + sort_particles!(gridsorter, particles[1], pia, 1) + + @test particles[1].index == [4, 1, 3, 2] + + @test maximum(abs.(particles[1][1].x - [3.55, -1.0, 9.0])) < 3.65e-15 # * eps() + @test particles[1][1].v == [-49.0, -20.0, 13.0] + @test particles[1][1].w == 4.0 + + @test maximum(abs.(particles[1][2].x - [20.5, -8.0, 7.5])) < 2 * eps() + @test particles[1][2].v == [-1.25, -1.5, 4.0] + @test particles[1][2].w == 1.0 + + @test maximum(abs.(particles[1][3].x - [23.0, 1.0, 3.0])) < 2 * eps() + @test particles[1][3].v == [20.0, 0.0, 2.0] + @test particles[1][3].w == 3.0 + + @test maximum(abs.(particles[1][4].x - [29.0, 6.0, -3.0])) < 2 * eps() + @test particles[1][4].v == [-11.0, -3.0, 1.0] + @test particles[1][4].w == 2.0 + + # non-contiguous version + # we create 30 particles in cell 2, but the pia structure only points to a subset of them + # and we test that only those are convected + particles = [ParticleVector(30)] + pia = ParticleIndexerArray(grid.n_cells, 1) + + for i in 1:10 + Merzbild.add_particle!(particles[1], i, 1.0, [1.0, 0.0, 0.0], [0.75, 0.0, 0.0]) + end + for i in 11:25 + Merzbild.add_particle!(particles[1], i, 10000.0, [-1000.0, 0.0, 0.0], [0.75, 0.0, 0.0]) + end + for i in 26:30 + Merzbild.add_particle!(particles[1], i, 1.0, [1.0, 0.0, 0.0], [0.75, 0.0, 0.0]) + end + + pia.n_total[1] = 15 + pia.indexer[2,1].n_local = 15 + pia.indexer[2,1].n_group1 = 10 + pia.indexer[2,1].start1 = 1 + pia.indexer[2,1].end1 = 10 + pia.indexer[2,1].n_group2 = 5 + pia.indexer[2,1].start2 = 26 + pia.indexer[2,1].end2 = 30 + pia.contiguous[1] = false + + compute_props!(particles, pia, species_data, phys_props) + @test abs(phys_props.n[2, 1] - 1.0 * 15) < eps() + flag = true + for i in 1:100 + if i != 2 + if abs(phys_props.n[i, 1] - 0) > eps() + flag = false + end + end + end + @test flag == true + @test abs(phys_props.v[1, 2, 1] - 1.0) < eps() + + convect_particles_and_compute_cell!(rng, grid, boundaries, particles[1], pia, 1, species_data, 1.0) + sort_particles!(gridsorter, particles[1], pia, 1) + + # all particles that are actually tracked are in cell 4 + compute_props!(particles, pia, species_data, phys_props) + @test abs(phys_props.n[4, 1] - 1.0 * 15) < eps() + + flag = true + for i in 1:100 + if i != 4 + if abs(phys_props.n[i, 1] - 0) > eps() + flag = false + end + end + end + @test flag == true + + @test abs(phys_props.v[1, 4, 1] - 1.0) < eps() end \ No newline at end of file diff --git a/test/test_debugging_functions.jl b/test/test_debugging_functions.jl index 39dda70..81924d9 100644 --- a/test/test_debugging_functions.jl +++ b/test/test_debugging_functions.jl @@ -114,11 +114,13 @@ # n_group1 is not equal to e1 - s1 + 1 pia = create_pia() - pia.indexer[4,1].n_group2 = 2 + pia.indexer[4,1].n_local += 2 - pia.indexer[4,1].n_group1 + pia.indexer[4,1].n_group1 = 2 @test check_pia_is_correct(pia, 1) == (false, 4) # n_group2 is not equal to e2 - s2 + 1 pia = create_pia() + pia.indexer[3,1].n_local += 2 - pia.indexer[3,1].n_group2 pia.indexer[3,1].n_group2 = 2 @test check_pia_is_correct(pia, 1) == (false, 3) diff --git a/test/test_grid_sorting.jl b/test/test_grid_sorting.jl index c1bfdbb..0009284 100644 --- a/test/test_grid_sorting.jl +++ b/test/test_grid_sorting.jl @@ -177,7 +177,55 @@ @test pia_fine.indexer[cell,1].n_local == counts[cell] end - phys_props::PhysProps = PhysProps(grid_fine.n_cells, 1, [], Tref=1) + phys_props = PhysProps(grid_fine.n_cells, 1, [], Tref=1) + compute_props!(particles, pia_fine, species_data, phys_props) + n_per_cell = Fnum * counts + + for i in 1:grid_fine.n_cells + @test abs(phys_props.n[i, 1] - n_per_cell[i]) < 2*eps() + @test phys_props.np[i, 1] == counts[i] + end + + # repeat the test case above + # but now we use the sorting routine which assumes particle cells have already been set during convection + # so the coordinates of the particles play no role whatsoever + # now we split the particles unevenly across 3 cells + # [0.0,2.0] - 3 particles (6, 7, 8) + # [2.0, 4.0] - 1 particle (5) + # [4.0, 6.0] - 0 particles () + # [6.0, 8.0] - 4 particles (1, 2, 3, 4) + + # we reset the index + particles[1].index = [1, 2, 3, 4, 5, 6, 7, 8] + particles[1].cell = [4, 4, 4, 4, 2, 1, 1, 1] + + for i in 1:4 + particles[1][i].x = SVector{3,Float64}(0.05, 0.0, 0.0) + end + for i in 5:5 + particles[1][i].x = SVector{3,Float64}(0.1, 0.0, 0.0) + end + for i in 6:8 + particles[1][i].x = SVector{3,Float64}(0.3, 0.0, 0.0) + end + + sort_particles!(gridsorter_fine, particles[1], pia_fine, 1) + @test particles[1].index == [6, 7, 8, 5, 1, 2, 3, 4] + counts = [3, 1, 0, 4] + starts = [1, 4, 0, 5] + ends = [3, 4, -1, 8] + + for cell in 1:4 + @test pia_fine.indexer[cell,1].n_group1 == counts[cell] + @test pia_fine.indexer[cell,1].start1 == starts[cell] + @test pia_fine.indexer[cell,1].end1 == ends[cell] + @test pia_fine.indexer[cell,1].n_group2 == 0 + @test pia_fine.indexer[cell,1].start2 == 0 + @test pia_fine.indexer[cell,1].end2 == -1 + @test pia_fine.indexer[cell,1].n_local == counts[cell] + end + + phys_props = PhysProps(grid_fine.n_cells, 1, [], Tref=1) compute_props!(particles, pia_fine, species_data, phys_props) n_per_cell = Fnum * counts diff --git a/test/test_ionization_sim.jl b/test/test_ionization_sim.jl index 56bfaae..b6361b4 100644 --- a/test/test_ionization_sim.jl +++ b/test/test_ionization_sim.jl @@ -327,13 +327,22 @@ # test charge neutrality charge_neutrality = true + mass_conservation = true + for t in 1:n_t if abs(sol["ndens"][1, 2, t] - sol["ndens"][1, 3, t]) / sol["ndens"][1, 3, t] > 1.5e-14 charge_neutrality = false end + + # densitiy of ions + neutrals should stay constant + if abs(sol["ndens"][1, 1, t] + sol["ndens"][1, 2, t] + - (n_dens_ions + n_dens_neutrals)) / (n_dens_ions + n_dens_neutrals) > 1.5e-14 + mass_conservation = false + end end @test charge_neutrality == true + @test mass_conservation == true # test compared to ref solution ref_sol_path = joinpath(@__DIR__, "data", "ionization_Ar_no_es_v2.nc") @@ -344,4 +353,158 @@ close(sol) rm(sol_path) + + + + # now we test with higher T0_e and E_tn and should get ionization + # and have inverse ratio of densities to check splitting + E_Tn = 500 + T0_e = Merzbild.eV * 10.0 # T_e(t=0) = 4.5eV + + n_dens_e = 1e23 + n_dens_ions = n_dens_e + n_dens_neutrals = 1e-3 * n_dens_e + Δt = 2e-13 + + particles = [ParticleVector(np_base_heavy), + ParticleVector(np_base_heavy), + ParticleVector(np_base_electrons)] + n_sampled = [0, 0, 0] + + for (index, (n_dens, T0, nv)) in enumerate(zip([n_dens_neutrals, n_dens_ions, n_dens_e], + [T0, T0, T0_e], + [nv_heavy, nv_heavy, nv_electrons])) + n_sampled[index] = sample_maxwellian_on_grid!(rng, particles[index], nv, species_data[index].mass, T0, n_dens, + 0.0, 1.0, 0.0, 1.0, 0.0, 1.0; + v_mult=4.5, cutoff_mult=8.0, noise=0.0, v_offset=[0.0, 0.0, 0.0]) + end + + pia = ParticleIndexerArray(n_sampled) + + + sol_path = joinpath(@__DIR__, "data", "tmp_ion_no_es_with_ion_inv_ratio.nc") + # exclude v and moments + ds = NCDataHolder(sol_path, ["v", "moments"], species_data, phys_props) + + compute_props!(particles, pia, species_data, phys_props) + write_netcdf(ds, phys_props, 0) + + collision_factors = create_collision_factors_array(3) + collision_data = CollisionData() + + if pia.n_total[1] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[1], pia, 1, 1, np_target_neutrals) + end + + if pia.n_total[2] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[2], pia, 1, 2, species_data, phys_props) + end + + if pia.n_total[3] > threshold_electrons + merge_octree_N2_based!(rng, oc_electrons, particles[3], pia, 1, 3, np_target_electrons) + end + + # neutral-neutral + Fnum_neutral_mean = n_dens_neutrals / np_target_neutrals + collision_factors[1,1,1].sigma_g_w_max = estimate_sigma_g_w_max(interaction_data[1,1], + species_data[1], T0, + Fnum_neutral_mean) + + s1 = index_neutral + s2 = index_electron + s3 = index_ion + # neutral-electron + + @test collision_factors[s1,s2,1].sigma_g_w_max == 0 + + estimate_sigma_g_w_max_ntc_n_e!(rng, collision_factors[s1,s2,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[s1], particles[s2], pia, 1, s1, s2, Δt, V, min_coll=2, n_loops=1) + + cf_e_n = collision_factors[s1,s2,1].sigma_g_w_max + @test cf_e_n > 0 + + # test that more samples lead to larger estimate + estimate_sigma_g_w_max_ntc_n_e!(rng, collision_factors[s1,s2,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[s1], particles[s2], pia, 1, s1, s2, Δt, V, min_coll=15, n_loops=5) + @test collision_factors[s1,s2,1].sigma_g_w_max > cf_e_n + + for ts in 1:n_t + # collide neutrals and neutrals + ntc!(rng, collision_factors[s1,s1,1], collision_data, interaction_data, particles[s1], pia, 1, s1, Δt, V) + + ntc_n_e!(rng, collision_factors[s1,s2,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[s1], particles[s2], particles[s3], pia, 1, s1, s2, s3, Δt, V) + + if pia.n_total[1] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[1], pia, 1, 1, np_target_neutrals) + end + + if pia.n_total[2] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[2], pia, 1, 2, species_data, phys_props) + end + + if pia.n_total[3] > threshold_electrons + merge_octree_N2_based!(rng, oc_electrons, particles[3], pia, 1, 3, np_target_electrons) + end + + accelerate_constant_field_x!(particles[index_electron], + pia, 1, index_electron, species_data, + E_field, Δt) + + compute_props!(particles, pia, species_data, phys_props) + write_netcdf(ds, phys_props, ts, sync_freq=1000) + end + close_netcdf(ds) + + sol = NCDataset(sol_path, "r") + + ndens_tot = n_dens_e + n_dens_ions + n_dens_neutrals + + # test overall number density conservation, we don't have it due to ionization + ndens_conservation = true + + for t in 1:n_t + if abs(sum(sol["ndens"][1, :, t]) - ndens_tot) / ndens_tot > 1.5e-14 + ndens_conservation = false + end + end + + @test ndens_conservation == false + + # test that we actually get ionization by looking at ion density at end of simulation + @test sol["ndens"][1, 2, n_t] > n_dens_ions + + # test charge neutrality + charge_neutrality = true + + mass_conservation = true + + for t in 1:n_t + if abs(sol["ndens"][1, 2, t] - sol["ndens"][1, 3, t]) / sol["ndens"][1, 3, t] > 2.0e-14 + charge_neutrality = false + end + + # densitiy of ions + neutrals should stay constant + if abs(sol["ndens"][1, 1, t] + sol["ndens"][1, 2, t] + - (n_dens_ions + n_dens_neutrals)) / (n_dens_ions + n_dens_neutrals) > 1.5e-14 + mass_conservation = false + end + end + + @test charge_neutrality == true + + @test mass_conservation == true + + # test compared to ref solution + ref_sol_path = joinpath(@__DIR__, "data", "ionization_Ar_no_es_v2_invratio.nc") + ref_sol = NCDataset(ref_sol_path, "r") + + @test maximum(abs.(ref_sol["ndens"][1, :, 1:n_t] .- sol["ndens"][1, :, 1:n_t])) < 4 * eps() + @test maximum(abs.(ref_sol["T"][1, :, 1:n_t] .- sol["T"][1, :, 1:n_t])) < 7.5e-14 + + close(sol) + rm(sol_path) end \ No newline at end of file diff --git a/test/test_ionization_sim_indexing.jl b/test/test_ionization_sim_indexing.jl new file mode 100644 index 0000000..f6fe19a --- /dev/null +++ b/test/test_ionization_sim_indexing.jl @@ -0,0 +1,195 @@ +@testset "ionization cross-section indexing" begin + # test that internal indexing of species works correctly + + # first we test with small T0_e + E_Tn = 100 # field strength in Tn + + T0 = 300.0 + T0_e = Merzbild.eV * 5.0 + n_dens_neutrals = 1e23 + n_dens_e = 1e-5 * n_dens_neutrals + n_dens_ions = n_dens_e + + E_field = E_Tn * n_dens_neutrals * 1e-21 # convert Tn to V/m + + Δt = 5e-14 + V = 1.0 + + merging_bin_split = OctreeBinMidSplit + + n_t = 100 + + threshold_electrons = 150 + np_target_electrons = 120 + + threshold_neutrals = 100 + threshold_ion = 50 + + np_target_neutrals = 100 + Nmerging_ions = 1 # 8-16 after merging + + seed = 123 + Random.seed!(seed) + rng = StableRNG(seed) + + + particles_data_path = joinpath(@__DIR__, "..", "data", "particles.toml") + interaction_data_path = joinpath(@__DIR__, "..", "data", "vhs.toml") + cs_data_path = joinpath(@__DIR__, "..", "data", "test_neutral_electron_data.xml") + + species_data::Vector{Species} = load_species_data(particles_data_path, ["Ar", "Ar+", "Xe+", "Xe", "e-"]) + interaction_data::Array{Interaction, 2} = load_interaction_data_with_dummy(interaction_data_path, species_data) + + n_e_interactions = load_electron_neutral_interactions(species_data, cs_data_path, + Dict("Ar" => "ConstantDB", + "Xe" => "ZeroDB"), + Dict("Ar" => ScatteringIsotropic, + "Xe" => ScatteringIsotropic), + Dict("Ar" => ElectronEnergySplitEqual, + "Xe" => ElectronEnergySplitEqual)) + + n_e_cs = create_computed_crosssections(n_e_interactions) + + @test n_e_interactions.n_neutrals == 2 + @test length(n_e_interactions.neutral_indexer) == 5 + @test n_e_interactions.neutral_indexer[1] == 1 + @test n_e_interactions.neutral_indexer[4] == 2 + + @test length(n_e_interactions.mass_ratios) == 2 + @test n_e_interactions.mass_ratios[1] > n_e_interactions.mass_ratios[2] # m_Ar < m_Xe + + for i in [2,3,5] + @test n_e_interactions.neutral_indexer[i] == -1 + end + + @test length(n_e_interactions.ionization) == 2 + @test maximum(n_e_interactions.ionization[1].data.sigma) > 0.0 + @test maximum(n_e_interactions.ionization[2].data.sigma) == 0.0 + + nv_heavy = 15 # init neutrals and ions on a coarser grid + nv_electrons = 20 + np_base_heavy = nv_heavy^3 # some initial guess on # of particles in simulation + np_base_electrons = nv_electrons^3 # some initial guess on # of particles in simulation + + oc = OctreeN2Merge(OctreeBinMidSplit; init_bin_bounds=OctreeInitBinMinMaxVel, max_Nbins=6000) + oc_electrons = OctreeN2Merge(merging_bin_split; init_bin_bounds=OctreeInitBinMinMaxVel, max_Nbins=6000) + mg_ions = GridN2Merge(Nmerging_ions, Nmerging_ions, Nmerging_ions, 3.5) + + particles = [ParticleVector(np_base_heavy), + ParticleVector(np_base_heavy), + ParticleVector(np_base_heavy), + ParticleVector(np_base_heavy), + ParticleVector(np_base_electrons)] + n_sampled = [0, 0, 0, 0, 0] + + for (index, (n_dens, T0, nv)) in enumerate(zip([n_dens_neutrals, n_dens_ions, n_dens_ions, n_dens_neutrals, n_dens_e], + [T0, T0, T0, T0, T0_e], + [nv_heavy, nv_heavy, nv_heavy, nv_heavy, nv_electrons])) + n_sampled[index] = sample_maxwellian_on_grid!(rng, particles[index], nv, species_data[index].mass, T0, n_dens, + 0.0, 1.0, 0.0, 1.0, 0.0, 1.0; + v_mult=3.5, cutoff_mult=8.0, noise=0.0, v_offset=[0.0, 0.0, 0.0]) + end + + pia = ParticleIndexerArray(n_sampled) + + phys_props::PhysProps = PhysProps(1, 5, [], Tref=T0) + + collision_factors = create_collision_factors_array(5) + collision_data = CollisionData() + compute_props!(particles, pia, species_data, phys_props) + + if pia.n_total[1] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[1], pia, 1, 1, np_target_neutrals) + end + + if pia.n_total[2] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[2], pia, 1, 2, species_data, phys_props) + end + + if pia.n_total[3] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[3], pia, 1, 3, species_data, phys_props) + end + + if pia.n_total[4] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[4], pia, 1, 4, np_target_neutrals) + end + + if pia.n_total[5] > threshold_electrons + merge_octree_N2_based!(rng, oc_electrons, particles[5], pia, 1, 5, np_target_electrons) + end + + # neutral-neutral + Fnum_neutral_mean = n_dens_neutrals / np_target_neutrals + collision_factors[1,1,1].sigma_g_w_max = estimate_sigma_g_w_max(interaction_data[1,1], + species_data[1], T0, + Fnum_neutral_mean) + + collision_factors[4,4,1].sigma_g_w_max = estimate_sigma_g_w_max(interaction_data[4,4], + species_data[4], T0, + Fnum_neutral_mean) + + for s1 in [1,4] + estimate_sigma_g_w_max_ntc_n_e!(rng, collision_factors[s1,5,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[s1], particles[5], pia, 1, s1, 5, Δt, V, min_coll=2, n_loops=1) + + cf_e_n = collision_factors[s1,5,1].sigma_g_w_max + @test cf_e_n > 0 + + # test that more samples lead to larger estimate + estimate_sigma_g_w_max_ntc_n_e!(rng, collision_factors[s1,5,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[s1], particles[5], pia, 1, s1, 5, Δt, V, min_coll=15, n_loops=5) + @test collision_factors[s1,5,1].sigma_g_w_max > cf_e_n + end + + + for ts in 1:n_t + + for (neu_species, ion_species) in zip([1, 4], [2, 3]) + ntc_n_e!(rng, collision_factors[neu_species,5,1], collision_data, interaction_data, + n_e_interactions, n_e_cs, + particles[neu_species], particles[5], particles[ion_species], pia, 1, neu_species, 5, ion_species, Δt, V) + end + + + if pia.n_total[1] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[1], pia, 1, 1, np_target_neutrals) + end + + if pia.n_total[2] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[2], pia, 1, 2, species_data, phys_props) + end + + if pia.n_total[3] > threshold_ion + merge_grid_based!(rng, mg_ions, particles[3], pia, 1, 3, species_data, phys_props) + end + + if pia.n_total[4] > threshold_neutrals + merge_octree_N2_based!(rng, oc, particles[4], pia, 1, 4, np_target_neutrals) + end + + if pia.n_total[5] > threshold_electrons + merge_octree_N2_based!(rng, oc_electrons, particles[5], pia, 1, 5, np_target_electrons) + end + + accelerate_constant_field_x!(particles[5], + pia, 1, 5, species_data, + E_field, Δt) + + compute_props!(particles, pia, species_data, phys_props) + end + + + # argon gets ionized + @test phys_props.n[1, 5] > n_dens_e + @test phys_props.n[1, 1] < n_dens_neutrals + @test phys_props.n[1, 2] > n_dens_ions + + # nothing happens to Xenon due to 0 cross-section + @test abs(phys_props.n[1, 4] - n_dens_neutrals) / n_dens_neutrals < 6e-15 + @test abs(phys_props.n[1, 3] - n_dens_ions) / n_dens_ions < 1.5e-14 + + # charge neutrality + @test abs(phys_props.n[1, 2] - phys_props.n[1, 5]) / phys_props.n[1, 5] < 6e-15 +end \ No newline at end of file diff --git a/test/test_pia_contiguous.jl b/test/test_pia_contiguous.jl index 55b2ff8..01bcc8e 100644 --- a/test/test_pia_contiguous.jl +++ b/test/test_pia_contiguous.jl @@ -35,6 +35,7 @@ @test pia.contiguous[1] == true + particles[1].cell[1:10] = [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] Merzbild.delete_particle_end_group1!(particles[1], pia, 1, 1) pia.contiguous[1] = false @@ -50,6 +51,8 @@ @test pia.indexer[1,1].n_group2 == 0 + @test particles[1].cell[1:9] == [11, 12, 13, 14, 15, 16, 17, 18, 19] + n_dens = 45.0 compute_props!(particles, pia, species_data, phys_props) @test phys_props.np[1,1] == 9.0 @@ -84,6 +87,8 @@ @test pia.contiguous[1] == true + particles[1].cell[1:10] = [11, 12, 13, 14, 15, 16, 17, 18, 19, 20] + Merzbild.delete_particle_end_group1!(particles[1], pia, 1, 1) Merzbild.delete_particle_end_group2!(particles[1], pia, 1, 1) pia.contiguous[1] = false @@ -94,6 +99,9 @@ bufferindex2 = particles[1].buffer[2] squash_pia!(particles, pia) + + @test particles[1].cell[1:8] == [11, 12, 13, 14, 16, 17, 18, 19] + @test pia.contiguous[1] = true @test pia.indexer[1,1].n_group1 == 4 diff --git a/test/test_surface_props_1D_uniform.jl b/test/test_surface_props_1D_uniform.jl index ff44308..6605b22 100644 --- a/test/test_surface_props_1D_uniform.jl +++ b/test/test_surface_props_1D_uniform.jl @@ -353,9 +353,67 @@ @test abs(surf_props.flux_reflected[1,1]/f[1] - -4.0) < 3*eps() @test abs(surf_props.flux_reflected[2,1]/f[2] - -3.0) < 3*eps() # we test the computation of the quantities before, so this should be sufficient to verify that + clear_props!(surf_props) + + # test same, but with convection that computes particle cell indices as we convect them + ppc = 4 + particles = [ParticleVector(ppc * grid.n_cells)] + + # 4 particles in cell 1, all of them reach the wall + for i in 1:4 + Merzbild.update_particle_buffer_new_particle!(particles[1], i) + + w = 1.0 + particles[1][i] = Particle(w, [-1.0, 0, 0], [0.1, 0.0, 0.0]) + end + + # 4 particles in cell 8, only one of them reaches the wall + for i in 5:8 + Merzbild.update_particle_buffer_new_particle!(particles[1], i) + + w = 3.0 + vx = 0.0 + if i == 5 + vx = 1.0 + end + particles[1][i] = Particle(w, [vx, 0, 3.0], [3.9, 0.0, 0.0]) + end + + pia.n_total[1] = 8 + pia.indexer[1,1].n_local = 4 + pia.indexer[1,1].start1 = 1 + pia.indexer[1,1].end1 = 4 + pia.indexer[1,1].n_group1 = 4 + + pia.indexer[1,1].start2 = -1 + pia.indexer[1,1].end2 = -1 + pia.indexer[1,1].n_group2 = 0 + pia.indexer[8,1].n_local = 4 + pia.indexer[8,1].start1 = 5 + pia.indexer[8,1].end1 = 8 + pia.indexer[8,1].n_group1 = 4 + pia.indexer[8,1].start2 = -1 + pia.indexer[8,1].end2 = -1 + pia.indexer[8,1].n_group2 = 0 + + clear_props!(surf_props) + + convect_particles_and_compute_cell!(rng, grid, boundaries, particles[1], pia, 1, species_data, surf_props, 0.2) + + @test surf_props.np == [4.0;1.0;;] + + + f = species_data[1].mass * surf_props.inv_areas / 0.2 + + @test abs(surf_props.flux_incident[1,1]/f[1] - 4.0) < 3*eps() + @test abs(surf_props.flux_incident[2,1]/f[2] - 3.0) < 3*eps() + @test abs(surf_props.flux_reflected[1,1]/f[1] - -4.0) < 3*eps() + @test abs(surf_props.flux_reflected[2,1]/f[2] - -3.0) < 3*eps() + # we test the computation of the quantities before, so this should be sufficient to verify that clear_props!(surf_props) + # test convection, non-contiguous particles = [ParticleVector(8)] @@ -413,4 +471,63 @@ @test abs(surf_props.flux_incident[2,1]/f[2] - 3.0) < 3*eps() @test abs(surf_props.flux_reflected[1,1]/f[1] - -2.0) < 3*eps() @test abs(surf_props.flux_reflected[2,1]/f[2] - -3.0) < 3*eps() + + clear_props!(surf_props) + # same as above, but we compute particle cell indices as we convect them + particles = [ParticleVector(8)] + + # 4 particles in cell 1, all of them reach the wall + for i in 1:4 + Merzbild.update_particle_buffer_new_particle!(particles[1], i) + + w = 1.0 + particles[1][i] = Particle(w, [-1.0, 0, 0], [0.1, 0.0, 0.0]) + end + + # 4 particles in cell 8, only one of them reaches the wall + for i in 5:8 + Merzbild.update_particle_buffer_new_particle!(particles[1], i) + + w = 3.0 + vx = 0.0 + if i == 5 + vx = 1.0 + end + particles[1][i] = Particle(w, [vx, 0, 3.0], [3.9, 0.0, 0.0]) + end + + + pia.n_total[1] = 8 + pia.indexer[1,1].n_local = 4 + pia.indexer[1,1].start1 = 1 + pia.indexer[1,1].end1 = 4 + pia.indexer[1,1].n_group1 = 4 + + pia.indexer[1,1].start2 = -1 + pia.indexer[1,1].end2 = -1 + pia.indexer[1,1].n_group2 = 0 + + pia.indexer[8,1].n_local = 4 + pia.indexer[8,1].start1 = 5 + pia.indexer[8,1].end1 = 6 + pia.indexer[8,1].n_group1 = 2 + + pia.indexer[8,1].start2 = 7 + pia.indexer[8,1].end2 = 8 + pia.indexer[8,1].n_group2 = 2 + + Merzbild.delete_particle_end!(particles[1], pia, 1, 1) + Merzbild.delete_particle_end!(particles[1], pia, 1, 1) + pia.contiguous[1] = false + + convect_particles_and_compute_cell!(rng, grid, boundaries, particles[1], pia, 1, species_data, surf_props, 0.2) + + @test surf_props.np == [2.0;1.0;;] + + f = species_data[1].mass * surf_props.inv_areas / 0.2 + + @test abs(surf_props.flux_incident[1,1]/f[1] - 2.0) < 3*eps() + @test abs(surf_props.flux_incident[2,1]/f[2] - 3.0) < 3*eps() + @test abs(surf_props.flux_reflected[1,1]/f[1] - -2.0) < 3*eps() + @test abs(surf_props.flux_reflected[2,1]/f[2] - -3.0) < 3*eps() end \ No newline at end of file