Variable-weight DSMC simulations

For variable-weight DSMC simulation not much is different compared to fixed-weight DSMC simulations. Existing routines for computation of macroscopic properties, I/O, and even collisions can be re-used, since they all do not explicitly assume fixed particle weights.

Two additional things are required, however: a way to sample particles with variable weights, and an algorithm to merge particles to avoid a blow-up of the number of particles due to collisions.

Sampling variable-weight particles

In order to obtain variable-weight particles, one can either

The grid-based sampling approaches compute the values of a VDF function on a discrete velocity grid and then re-interpret the grid points and associated VDF values as particles. Thus one can get a good resolution of high-velocity tails (with the particles there having small weights). The function constructs the velocity grid via a Cartesian product of discrete grids for each velocity direction. In each direction, the interval $[-v_m\sqrt{2kT/m},v_m\sqrt{2kT/m}]$ is discretized into nv subintervals. Here the multiplier $v_m$ is specified by the v_mult parameter; so the discrete velocity grid spans several mean thermal velocities. In addition, only particles with a speed less than $c_m \sqrt{2kT/m}$ are actually created, where $c_m$ is specified by the cutoff_mult parameter. Finally, numerical noise may be added to the velocity via the noise parameter, and a velocity offset via the v_offset parameter; the spatial distribution of the particles is randomized in the cell. The function returns the total number of particles sampled (which may be less than nv^3 depending on the choice of the value of cutoff_mult).

The sample_particles_phase_box_weighted, unlike the grid-based sampling, is not deterministic, and samples particle velocities uniformly in a box in phase space (with the extent set similar to the grid-based sampling). The particle weights are computed based on the value of a prescribed Maxwellian at that velocity, and re-scaled to ensure the correct number density; the streaming velocity is added afterwards.

An example of such sampling is shown below:

using Merzbildusing Randomseed = 1Random.seed!(seed)rng = Xoshiro(seed)# our number densityndens = 1e23 # our temperatureT0 = 300.0  # load species and interaction dataspecies_data = load_species_data(joinpath(MERZBILD_DATA_PATH, "particles.toml"), "Ar")interaction_data = load_interaction_data(joinpath(MERZBILD_DATA_PATH, "pseudo_maxwell.toml"), species_data)# number of velocity grid points in each velocity directionnv = 20  # some initial guess on # of particles in simulationnp_base = 40^3  particles = [ParticleVector{0}(np_base)]# sample from a BKW distribution at t=0vdf0 = (vx, vy, vz) -> bkw(vx, vy, vz, species_data[1].mass, T0, 0.0)# returns number of particles sampledn_sampled = sample_on_grid!(rng, vdf0, particles[1], nv, species_data[1].mass, T0, ndens,                            0.0, 1.0, 0.0, 1.0, 0.0, 1.0;                            v_mult=3.5, cutoff_mult=3.5, noise=0.0, v_offset=[0.0, 0.0, 0.0])# create the ParticleIndexerArraypia = ParticleIndexerArray(n_sampled)

Colliding variable-weight particles

As mentioned above, the existing DSMC routines, unless they have equal_weight explicitly in their name, take care of the required particle splitting. For variable-weight NTC simulations, one needs to estimate $(\sigma g w)_{max}$. The simplest approach is to use estimate_sigma_g_w_max function, but it requires a (fixed) value of Fnum - so one can simply compute it as ndens/n_sampled (i.e. the number density divided by the total number of sampled particles). If we merge after sampling, an effective Fnum is then computed as ndens/n_post_merge (i.e. number density divided by the total number of post-merge particles). The ntc! function has a keyword parameter dw_tol that governs when particle collisions are treated as equal-weight so that small discrepancies in particle weights do not cause unnecessary splitting.

Merging variable-weight particles

In order to merge variable-weight particles, one needs to set up a merging algorithm and the associated data structures and parameters. A more detailed overview of merging algorithms will appear later in a Tutorials section; here, the Octree merging approach of Martin and Cambier (2016) is used. The algorithm groups particles into bins in velocity space, recursively refining the bins until the number of post-merge particles reaches the prescribed value. It then performs an $N:M$ merge in each bin, replacing all particles in a bin with M particles (currently, M=1 and M=2 are supported, N:1 merging is made conservative via re-scaling of post-merge particles' properties).

To set up this merging algorithm, one needs to create a OctreeMerge{D,M} instance, specifying

  • The extent of the first (root) bin - whether it accounts for the extent of the particles or whether it is just taken to be very large
  • How bins are split (along the middle velocity, the mean velocity, or the median velocity)
  • Whether the velocity bounds of each sub-bin are recomputed based on the particles in the bin or are based purely on the bounds of the parent bin and the splitting velocity
  • Maximum number of bins
  • Maximum refinement depth
  • M: the number of post-merge particles in a bin (defaults to 2)

It is important that the dimension D of the merging instance coincides with that used for the particle positions. For example, we can create an octree merging instance that splits velocity bins across the middle, sets the root bin bounds to the bounding box of the particle velocities, and inherits the bin bounds of the parent bin when splitting a bin. We also immediately merge our particles, setting a target particle number of 100.

# set up the merging algorithm: 0-D particles, 2 post-merge particles in a binoc = OctreeMerge{0,2}(OctreeBinMidSplit; init_bin_bounds=OctreeInitBinMinMaxVel,                   bin_bounds_compute=OctreeBinBoundsInherit, max_Nbins=6000)# set NtargetNtarget = 100# perform mergingmerge_octree!(rng, oc, particles[1], pia, 1, 1, Ntarget)

We can check the number of particles after the merging procedure by looking at pia.indexer[1,1].n_local (the number of particles of species 1 in cell 1), which should be equal to 100. Important: when estimating collision factors, one might need an estimate for the average weight of a particle; this needs to be computed after the merging is performed on the initial simulation, otherwise the estimates might be significantly off if merging reduces the number of particles by a large amount.

Example: bringing it all together

Now we can put the variable-weight collisions and merging into a timestep loop, performing merging by when pia.indexer[1,1].n_local > Nthreshold, where Nthreshold is the threshold number of particles. Alternatively, one can check phys_props.np[1,1] > Nthreshold, in case the physical properties are computed frequently enough. The following simulation computes the BKW relaxation problem.

using Merzbildusing Randomseed = 1Random.seed!(seed)rng = Xoshiro(seed)# our number densityndens = 1e23 # our temperatureT0 = 300.0  # load species and interaction dataspecies_data = load_species_data(joinpath(MERZBILD_DATA_PATH, "particles.toml"), "Ar")interaction_data = load_interaction_data(joinpath(MERZBILD_DATA_PATH, "pseudo_maxwell.toml"), species_data)# number of velocity grid points in each velocity directionnv = 20  # some initial guess on # of particles in simulationnp_base = 40^3  particles = [ParticleVector{0}(np_base)]# sample from a BKW distribution at t=0vdf0 = (vx, vy, vz) -> bkw(vx, vy, vz, species_data[1].mass, T0, 0.0)# returns number of particles sampledn_sampled = sample_on_grid!(rng, vdf0, particles[1], nv, species_data[1].mass, T0, ndens,                            0.0, 1.0, 0.0, 1.0, 0.0, 1.0;                            v_mult=3.5, cutoff_mult=3.5, noise=0.0, v_offset=[0.0, 0.0, 0.0])# create the ParticleIndexerArraypia = ParticleIndexerArray(n_sampled)# set up the merging algorithm: 0-D particles, 2 post-merge particles in a binoc = OctreeMerge{0,2}(OctreeBinMidSplit; init_bin_bounds=OctreeInitBinMinMaxVel,                   bin_bounds_compute=OctreeBinBoundsInherit, max_Nbins=6000)# set Ntarget for mergingNtarget = 100# set Nthreshold for mergingNthreshold = 120# perform initial mergemerge_octree!(rng, oc, particles[1], pia, 1, 1, Ntarget)# set some reference valuessigma_ref = π * (interaction_data[1,1].vhs_d^2)vref = sqrt(2 * k_B * T0 / species_data[1].mass)Lref = 1.0 / (ndens * sigma_ref)tref = Lref / vref# set up computation of physical propertiesphys_props = PhysProps(1, 1)compute_props!(particles, pia, species_data, phys_props)ds = NCDataHolder("output_bkw_octree.nc", species_data, phys_props)write_netcdf(ds, phys_props, 0)collision_factors = create_collision_factors_array(1)  # 1-species, 0-Dcollision_data = CollisionData()# estimate average Fnum - we need the post-merge value!Fnum = ndens/Ntarget# estimate (sigma_g_w)_maxestimate_sigma_g_w_max!(collision_factors, interaction_data, species_data, [T0], Fnum)# set number of timesteps and scaled timestepn_t = 100dt_scaled = 0.025# set timestep and cell volume (0D so V=1)Δt = dt_scaled * trefV = 1.0# start time loopfor ts in 1:n_t    # collide    ntc!(rng, collision_factors[1,1,1], collision_data, interaction_data, particles[1], pia, 1, 1, Δt, V)    # check if we need to merge    if pia.indexer[1,1].n_local > Nthreshold        merge_octree!(rng, oc, particles[1], pia, 1, 1, Ntarget)    end    # print number of particles every 10 timesteps    if ts % 10 == 0        println("t=$ts, np=$(pia.indexer[1,1].n_local)")    end        compute_props!(particles, pia, species_data, phys_props)    write_netcdf(ds, phys_props, ts)endclose_netcdf(ds)