1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
use super::*;
// ── Update system: rate-based insertion ─────────────────────────────────────
/// Update system for rate-based particle insertion during the simulation run.
///
/// Checks each registered [`RateInsertEntry`] against the current timestep, interval,
/// start/end bounds, and total limit. Uses a [`SpatialHash`] for O(1) overlap detection
/// when placing new particles. Runs in `ParticleSimScheduleSet::PreInitialIntegration`.
#[allow(clippy::too_many_arguments)]
pub fn dem_rate_insert(
comm: Res<CommResource>,
domain: Res<Domain>,
mut atom: ResMut<Atom>,
registry: Res<AtomDataRegistry>,
run_state: Res<RunState>,
mut rate_state: ResMut<RateInsertState>,
mut comm_state: ResMut<CurrentState<CommState>>,
) {
// Rate insertion runs on EVERY rank (born-in-owner): each rank generates the
// identical candidate stream from a step-derived seed and stores only the
// candidates that fall inside its own subdomain, so a new atom is born inside
// its owner and never needs a multi-hop exchange.
if rate_state.entries.is_empty() {
return;
}
let step = run_state.total_cycle;
let mut any_to_insert = false;
// Quick check if any entry needs insertion this step (before stripping ghosts)
// Quick check if any entry needs insertion this step
for entry in rate_state.entries.iter() {
let interval = entry.config.rate_interval.unwrap_or(1);
let start = entry.config.rate_start.unwrap_or(0);
if step < start {
continue;
}
if let Some(end) = entry.config.rate_end {
if step > end {
continue;
}
}
if let Some(limit) = entry.config.rate_limit {
if entry.total_inserted >= limit {
continue;
}
}
let steps_since_start = step - start;
if interval == 0 || steps_since_start % interval == 0 {
any_to_insert = true;
break;
}
}
if !any_to_insert {
return;
}
// Strip ghost atoms before inserting new local atoms.
// New atoms are appended at atom.len(), which must equal nlocal so that
// the subsequent borders() truncate_to_nlocal() doesn't discard them.
if atom.nghost > 0 {
// Keep the core ghost suffix and every DIRT extension in one
// transaction. In particular, a plugin registered after setup has a
// row here too; truncating Atom and AtomDataRegistry separately would
// leave a panic/error window between the two mutations.
ParticleStore::new(&mut atom, ®istry)
.discard_ghosts()
.expect("rate insertion requires an aligned local/ghost particle layout");
}
// Base tag must be globally consistent across ranks. (No all_reduce_max in
// the backend, so reduce -max with min and negate.) Tags advance only for
// accepted candidates, matching immediate insertion. Acceptance is globally
// replicated, so this remains rank-count invariant without another
// collective and rejected attempts cannot leave gaps in the accepted stream.
let local_max_tag = if atom.tag.is_empty() {
-1.0
} else {
atom.get_max_tag() as f64
};
let base_tag = (-comm.all_reduce_min_f64(-local_max_tag)) as i64 + 1;
let mut tag_cursor: u32 = base_tag.max(0) as u32;
for entry_idx in 0..rate_state.entries.len() {
let interval = rate_state.entries[entry_idx]
.config
.rate_interval
.unwrap_or(1);
let start = rate_state.entries[entry_idx].config.rate_start.unwrap_or(0);
let (rate, _, _) = validate_rate_insert_config(
&rate_state.entries[entry_idx].config,
"rate-based [[particles.insert]]",
)
.expect("rate insertion was validated during fallible plugin preflight");
if step < start {
continue;
}
if let Some(end) = rate_state.entries[entry_idx].config.rate_end {
if step > end {
continue;
}
}
if let Some(limit) = rate_state.entries[entry_idx].config.rate_limit {
if rate_state.entries[entry_idx].total_inserted >= limit {
continue;
}
}
let steps_since_start = step - start;
if interval > 0 && steps_since_start % interval != 0 {
continue;
}
// How many to insert this step
let mut to_insert = rate;
if let Some(limit) = rate_state.entries[entry_idx].config.rate_limit {
let remaining = limit - rate_state.entries[entry_idx].total_inserted;
to_insert = to_insert.min(remaining);
}
let prepared = &rate_state.entries[entry_idx].prepared;
// Seed the candidate stream from (config seed, step, entry) so it is
// identical on every rank yet varies between insertion events.
let stream_seed = prepared.seed
^ (step as u64).wrapping_mul(0x9E3779B97F4A7C15)
^ (entry_idx as u64).wrapping_mul(0xD1B54A32D192ED03);
// Replicated overlap scratch: ONLY the positions/radii accepted THIS step
// (the global set of new atoms). It is identical on every rank because the
// candidate stream is seeded identically and the scratch is the same on
// all ranks. This is what keeps accept/reject — and therefore the RNG
// advancement and the global accept count `inserted` — in lock-step across
// ranks, so the collective borders()/exchange() triggered below stay
// synchronized.
//
// NOTE: existing local atoms are deliberately NOT added to the scratch.
// They differ per rank, so including them would make accept/reject (and
// the RNG stream) diverge across ranks and desync the collectives. New
// atoms may therefore be born overlapping already-present particles; the
// contact model resolves that initial overlap via repulsion. Rate-insert
// regions are normally placed in free space (e.g. above a settled bed),
// so this is rare in practice.
let mut candidates =
CandidateGenerator::new(prepared, &domain, stream_seed, to_insert as usize);
let mut inserted = 0u32; // accepted globally
let mut local_inserted = 0u32; // stored on this rank
let mut attempts = 0u32;
let max_attempts = to_insert * 100;
while inserted < to_insert && attempts < max_attempts {
attempts += 1;
let Some(candidate) = candidates.next(prepared) else {
continue;
};
let tag = tag_cursor;
tag_cursor = tag_cursor.wrapping_add(1);
// Store only if this rank owns the position.
if owns_position(&domain, &candidate.pos) {
insert_single_particle(&mut atom, ®istry, candidate.particle(prepared, tag));
local_inserted += 1;
}
inserted += 1;
}
rate_state.entries[entry_idx].total_inserted += inserted;
let _ = local_inserted;
if inserted > 0 {
// Force full ghost rebuild on EVERY rank if any rank inserted, so the
// collective borders()/exchange() stay in lock-step. `inserted` is the
// global accept count and is identical on all ranks (replicated stream).
comm_state.0 = CommState::FullRebuild;
}
if inserted > 0 && attempts >= max_attempts && comm.rank() == 0 {
eprintln!(
"WARNING: Rate insertion at step {} only placed {}/{} particles (max attempts reached)",
step, inserted, to_insert
);
}
}
}