PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
test_runtime_kernels.c
Go to the documentation of this file.
1/**
2 * @file test_runtime_kernels.c
3 * @brief C unit tests for runtime, particle, wall, and walltime-guard helpers.
4 */
5
6#include "test_support.h"
7
8#include "BC_Handlers.h"
9#include "ParticleMotion.h"
10#include "ParticlePhysics.h"
11#include "ParticleSwarm.h"
12#include "interpolation.h"
13#include "initialcondition.h"
14#include "les.h"
15#include "runloop.h"
16#include "setup.h"
17#include "wallfunction.h"
18#include "walkingsearch.h"
19
20/**
21 * @brief Synchronizes the minimal runtime fixture's global fields into their persistent local ghosts.
22 */
23static PetscErrorCode SyncRuntimeFieldGhosts(UserCtx *user)
24{
25 PetscFunctionBeginUser;
26 PetscCall(DMGlobalToLocalBegin(user->da, user->P, INSERT_VALUES, user->lP));
27 PetscCall(DMGlobalToLocalEnd(user->da, user->P, INSERT_VALUES, user->lP));
28 PetscCall(DMGlobalToLocalBegin(user->da, user->Psi, INSERT_VALUES, user->lPsi));
29 PetscCall(DMGlobalToLocalEnd(user->da, user->Psi, INSERT_VALUES, user->lPsi));
30 PetscCall(DMGlobalToLocalBegin(user->da, user->Diffusivity, INSERT_VALUES, user->lDiffusivity));
31 PetscCall(DMGlobalToLocalEnd(user->da, user->Diffusivity, INSERT_VALUES, user->lDiffusivity));
32 PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
33 PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
34 PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
35 PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont));
36 PetscCall(DMGlobalToLocalBegin(user->fda, user->DiffusivityGradient, INSERT_VALUES, user->lDiffusivityGradient));
37 PetscCall(DMGlobalToLocalEnd(user->fda, user->DiffusivityGradient, INSERT_VALUES, user->lDiffusivityGradient));
38 PetscFunctionReturn(0);
39}
40
41/**
42 * @brief Seeds one localized swarm particle with the cell, position, weight, and status data used by runtime tests.
43 */
44static PetscErrorCode SeedSingleParticle(UserCtx *user,
45 PetscInt ci,
46 PetscInt cj,
47 PetscInt ck,
48 PetscReal x,
49 PetscReal y,
50 PetscReal z,
51 PetscReal wx,
52 PetscReal wy,
53 PetscReal wz,
54 PetscInt status_value)
55{
56 PetscReal *positions = NULL;
57 PetscReal *weights = NULL;
58 PetscInt *cell_ids = NULL;
59 PetscInt *status = NULL;
60
61 PetscFunctionBeginUser;
62 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
63 PetscCall(DMSwarmGetField(user->swarm, "weight", NULL, NULL, (void **)&weights));
64 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
65 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
66
67 positions[0] = x;
68 positions[1] = y;
69 positions[2] = z;
70 weights[0] = wx;
71 weights[1] = wy;
72 weights[2] = wz;
73 cell_ids[0] = ci;
74 cell_ids[1] = cj;
75 cell_ids[2] = ck;
76 status[0] = status_value;
77
78 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
79 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
80 PetscCall(DMSwarmRestoreField(user->swarm, "weight", NULL, NULL, (void **)&weights));
81 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
82 PetscFunctionReturn(0);
83}
84/**
85 * @brief Tests particle distribution remainder handling across ranks.
86 */
87
89{
90 PetscInt particles_per_rank = -1;
91 PetscInt remainder = -1;
92
93 PetscFunctionBeginUser;
94 PetscCall(DistributeParticles(10, 0, 3, &particles_per_rank, &remainder));
95 PetscCall(PicurvAssertIntEqual(4, particles_per_rank, "rank 0 should receive one remainder particle"));
96 PetscCall(PicurvAssertIntEqual(1, remainder, "remainder should be reported correctly"));
97
98 PetscCall(DistributeParticles(10, 2, 3, &particles_per_rank, &remainder));
99 PetscCall(PicurvAssertIntEqual(3, particles_per_rank, "last rank should receive base particle count"));
100 PetscCall(PicurvAssertIntEqual(1, remainder, "remainder should remain unchanged across ranks"));
101 PetscFunctionReturn(0);
102}
103/**
104 * @brief Tests basic particle-inside-bounding-box classification cases.
105 */
106
108{
109 BoundingBox bbox;
110 Particle particle;
111
112 PetscFunctionBeginUser;
113 PetscCall(PetscMemzero(&bbox, sizeof(bbox)));
114 PetscCall(PetscMemzero(&particle, sizeof(particle)));
115
116 bbox.min_coords.x = 0.0;
117 bbox.min_coords.y = 0.0;
118 bbox.min_coords.z = 0.0;
119 bbox.max_coords.x = 1.0;
120 bbox.max_coords.y = 2.0;
121 bbox.max_coords.z = 3.0;
122
123 particle.loc.x = 0.25;
124 particle.loc.y = 1.0;
125 particle.loc.z = 2.5;
126 PetscCall(PicurvAssertBool(IsParticleInsideBoundingBox(&bbox, &particle), "particle should be inside bounding box"));
127
128 particle.loc.x = 1.5;
129 PetscCall(PicurvAssertBool((PetscBool)!IsParticleInsideBoundingBox(&bbox, &particle), "particle should be outside bounding box"));
130 PetscFunctionReturn(0);
131}
132/**
133 * @brief Tests particle weight updates against expected ratios.
134 */
135
137{
138 Particle particle;
139 PetscReal distances[NUM_FACES] = {1.0, 3.0, 2.0, 2.0, 4.0, 1.0};
140 PetscReal clamped[NUM_FACES] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
141
142 PetscFunctionBeginUser;
143 PetscCall(PetscMemzero(&particle, sizeof(particle)));
144 PetscCall(UpdateParticleWeights(distances, &particle));
145 PetscCall(PicurvAssertRealNear(0.25, particle.weights.x, 1.0e-12, "x interpolation weight"));
146 PetscCall(PicurvAssertRealNear(0.5, particle.weights.y, 1.0e-12, "y interpolation weight"));
147 PetscCall(PicurvAssertRealNear(0.2, particle.weights.z, 1.0e-12, "z interpolation weight"));
148
149 PetscCall(UpdateParticleWeights(clamped, &particle));
150 PetscCall(PicurvAssertRealNear(0.5, particle.weights.x, 1.0e-12, "clamped x weight remains centered"));
151 PetscCall(PicurvAssertRealNear(0.5, particle.weights.y, 1.0e-12, "clamped y weight remains centered"));
152 PetscCall(PicurvAssertRealNear(0.5, particle.weights.z, 1.0e-12, "clamped z weight remains centered"));
153 PetscFunctionReturn(0);
154}
155/**
156 * @brief Tests particle position updates without Brownian forcing.
157 */
158
160{
161 SimCtx *simCtx = NULL;
162 UserCtx *user = NULL;
163 Particle particle;
164
165 PetscFunctionBeginUser;
166 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
167 simCtx->dt = 0.25;
168
169 PetscCall(PetscMemzero(&particle, sizeof(particle)));
170 particle.loc.x = 1.0;
171 particle.loc.y = -2.0;
172 particle.loc.z = 3.0;
173 particle.vel.x = 0.5;
174 particle.vel.y = -1.0;
175 particle.vel.z = 2.0;
176 particle.diffusivitygradient.x = 0.1;
177 particle.diffusivitygradient.y = 0.2;
178 particle.diffusivitygradient.z = -0.3;
179 particle.diffusivity = 0.0;
180
181 PetscCall(UpdateParticlePosition(user, &particle));
182 PetscCall(PicurvAssertRealNear(1.15, particle.loc.x, 1.0e-12, "x position update"));
183 PetscCall(PicurvAssertRealNear(-2.2, particle.loc.y, 1.0e-12, "y position update"));
184 PetscCall(PicurvAssertRealNear(3.425, particle.loc.z, 1.0e-12, "z position update"));
185
186 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
187 PetscFunctionReturn(0);
188}
189/**
190 * @brief Tests particle position updates driven only by diffusivity-gradient drift.
191 */
192
194{
195 SimCtx *simCtx = NULL;
196 UserCtx *user = NULL;
197 Particle particle;
198
199 PetscFunctionBeginUser;
200 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
201 simCtx->dt = 0.5;
202
203 PetscCall(PetscMemzero(&particle, sizeof(particle)));
204 particle.loc.x = 0.25;
205 particle.loc.y = 0.5;
206 particle.loc.z = 0.75;
207 particle.vel.x = 0.0;
208 particle.vel.y = 0.0;
209 particle.vel.z = 0.0;
210 particle.diffusivitygradient.x = 0.2;
211 particle.diffusivitygradient.y = -0.1;
212 particle.diffusivitygradient.z = 0.3;
213 particle.diffusivity = 0.0;
214
215 PetscCall(UpdateParticlePosition(user, &particle));
216 PetscCall(PicurvAssertRealNear(0.35, particle.loc.x, 1.0e-12, "drift-only x position update"));
217 PetscCall(PicurvAssertRealNear(0.45, particle.loc.y, 1.0e-12, "drift-only y position update"));
218 PetscCall(PicurvAssertRealNear(0.90, particle.loc.z, 1.0e-12, "drift-only z position update"));
219
220 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
221 PetscFunctionReturn(0);
222}
223/**
224 * @brief Tests IEM relaxation updates for particle-carried fields.
225 */
226
227static PetscErrorCode TestUpdateParticleFieldIEMRelaxation(void)
228{
229 PetscReal psi = 1.0;
230 const PetscReal dt = 0.5;
231 const PetscReal diffusivity = 0.2;
232 const PetscReal mean_val = 3.0;
233 const PetscReal cell_vol = 8.0;
234 const PetscReal c_model = 2.0;
235 PetscReal unchanged = 7.0;
236
237 PetscFunctionBeginUser;
238 PetscCall(UpdateParticleField("Psi", dt, &psi, diffusivity, mean_val, cell_vol, c_model));
239 PetscCall(PicurvAssertRealNear(
240 mean_val + (1.0 - mean_val) * PetscExpReal(-(c_model * diffusivity / PetscPowReal(cell_vol, 0.6666667)) * dt),
241 psi,
242 1.0e-12,
243 "IEM update should match analytical relaxation"));
244
245 PetscCall(UpdateParticleField("UnrelatedField", dt, &unchanged, diffusivity, mean_val, cell_vol, c_model));
246 PetscCall(PicurvAssertRealNear(7.0, unchanged, 1.0e-12, "unknown field should remain unchanged"));
247 PetscFunctionReturn(0);
248}
249/**
250 * @brief Tests that non-Ucont requests do not modify interior field initialization.
251 */
252
254{
255 SimCtx *simCtx = NULL;
256 UserCtx *user = NULL;
257
258 PetscFunctionBeginUser;
259 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
260 PetscCall(VecSet(user->Ucont, 7.0));
261
262 PetscCall(SetInitialInteriorField(user, "P"));
263 PetscCall(PicurvAssertVecConstant(user->Ucont, 7.0, 1.0e-12, "non-Ucont request should not modify Ucont"));
264
265 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
266 PetscFunctionReturn(0);
267}
268/**
269 * @brief Tests cartesian Constant IC: Cart2Contra sets contravariant flux via metric dot product.
270 *
271 * The minimal runtime fixture uses identity face metrics. UniformCart2Contra
272 * therefore maps (u,v,w)=(0,0,2) directly to a Zeta flux of 2.
273 */
275{
276 SimCtx *simCtx = NULL;
277 UserCtx *user = NULL;
278 Cmpnts ***ucont = NULL;
279
280 PetscFunctionBeginUser;
281 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
283 simCtx->InitialConstantContra.x = 0.0;
284 simCtx->InitialConstantContra.y = 0.0;
285 simCtx->InitialConstantContra.z = 2.0;
286 PetscCall(VecSet(user->Ucont, 0.0));
287
288 PetscCall(SetInitialInteriorField(user, "Ucont"));
289 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
290 PetscCall(PicurvAssertRealNear(0.0, ucont[1][1][1].x, 1.0e-10, "Xi flux is zero for pure-z velocity on Cartesian grid"));
291 PetscCall(PicurvAssertRealNear(0.0, ucont[1][1][1].y, 1.0e-10, "Eta flux is zero for pure-z velocity on Cartesian grid"));
292 PetscCall(PicurvAssertRealNear(2.0, ucont[1][1][1].z, 1.0e-10, "Zeta flux follows the identity fixture metric"));
293 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
294
295 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
296 PetscFunctionReturn(0);
297}
298
299/**
300 * @brief Tests curvilinear Constant IC: flow_direction selects the streamwise axis.
301 *
302 * icVelocityPhysical=2.0, FLOW_DIR_POS_ZETA. The minimal runtime fixture uses
303 * identity face metrics, so ucont.z = 2.0.
304 * Xi and Eta components are left at zero regardless of InitialConstantContra.
305 */
307{
308 SimCtx *simCtx = NULL;
309 UserCtx *user = NULL;
310 Cmpnts ***ucont = NULL;
311
312 PetscFunctionBeginUser;
313 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
315 simCtx->icVelocityPhysical = 2.0;
317 user->GridOrientation = 1;
318 /* InitialConstantContra intentionally non-zero to confirm only z-component is set */
319 simCtx->InitialConstantContra.x = 99.0;
320 simCtx->InitialConstantContra.y = 99.0;
321 simCtx->InitialConstantContra.z = 99.0;
322 PetscCall(VecSet(user->Ucont, 0.0));
323
324 PetscCall(SetInitialInteriorField(user, "Ucont"));
325 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
326 PetscCall(PicurvAssertRealNear(0.0, ucont[1][1][1].x, 1.0e-10, "Xi flux is zero in curvilinear Zeta mode"));
327 PetscCall(PicurvAssertRealNear(0.0, ucont[1][1][1].y, 1.0e-10, "Eta flux is zero in curvilinear Zeta mode"));
328 PetscCall(PicurvAssertRealNear(2.0, ucont[1][1][1].z, 1.0e-10, "Zeta flux follows the identity fixture metric"));
329 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
330
331 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
332 PetscFunctionReturn(0);
333}
334
335/**
336 * @brief Tests zero IC clears physical-cell contravariant velocity.
337 */
339{
340 SimCtx *simCtx = NULL;
341 UserCtx *user = NULL;
342 Cmpnts ***ucont = NULL;
343
344 PetscFunctionBeginUser;
345 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
347 PetscCall(VecSet(user->Ucont, 7.0));
348
349 PetscCall(SetInitialInteriorField(user, "Ucont"));
350 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
351 PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].x, 1.0e-12, "zero IC clears interior Xi flux"));
352 PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].y, 1.0e-12, "zero IC clears interior Eta flux"));
353 PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].z, 1.0e-12, "zero IC clears interior Zeta flux"));
354 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
355
356 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
357 PetscFunctionReturn(0);
358}
359
360/**
361 * @brief Tests Poiseuille IC follows the discrete cross-stream profile and reaches zero at edges.
362 */
364{
365 SimCtx *simCtx = NULL;
366 UserCtx *user = NULL;
367 Cmpnts ***ucont = NULL;
368
369 PetscFunctionBeginUser;
370 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
372 simCtx->icVelocityPhysical = 3.0;
374 user->GridOrientation = 1;
375 PetscCall(VecZeroEntries(user->Ucont));
376
377 PetscCall(SetInitialInteriorField(user, "Ucont"));
378 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
379 PetscCall(PicurvAssertRealNear(3.0 * 64.0 / 81.0, ucont[2][2][2].z, 1.0e-12,
380 "Poiseuille IC follows the discrete center-adjacent profile value"));
381 PetscCall(PicurvAssertRealNear(0.0, ucont[2][1][1].z, 1.0e-12, "Poiseuille IC is zero at cross-stream edge"));
382 PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].x, 1.0e-12, "Poiseuille IC leaves Xi flux zero"));
383 PetscCall(PicurvAssertRealNear(0.0, ucont[2][2][2].y, 1.0e-12, "Poiseuille IC leaves Eta flux zero"));
384 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
385
386 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
387 PetscFunctionReturn(0);
388}
389
390/**
391 * @brief Tests spatially varying Cartesian-field conversion to contravariant fluxes.
392 */
393static PetscErrorCode TestCart2ContraConvertsCartesianField(void)
394{
395 SimCtx *simCtx = NULL;
396 UserCtx *user = NULL;
397 Cmpnts ***ucat = NULL;
398 Cmpnts ***ucont = NULL;
399
400 PetscFunctionBeginUser;
401 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
402 PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
403 for (PetscInt k = 0; k < 6; k++)
404 for (PetscInt j = 0; j < 6; j++)
405 for (PetscInt i = 0; i < 6; i++)
406 ucat[k][j][i] = (Cmpnts){(PetscReal)i, 2.0, 3.0};
407 PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
408 PetscCall(UpdateLocalGhosts(user, "Ucat"));
409 PetscCall(Cart2Contra(user));
410
411 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
412 PetscCall(PicurvAssertRealNear(1.5, ucont[1][1][1].x, 1.0e-10,
413 "Xi flux uses face-interpolated Ucat"));
414 PetscCall(PicurvAssertRealNear(2.0, ucont[1][1][1].y, 1.0e-10,
415 "Eta flux uses Cartesian y velocity"));
416 PetscCall(PicurvAssertRealNear(3.0, ucont[1][1][1].z, 1.0e-10,
417 "Zeta flux uses Cartesian z velocity"));
418 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
419
420 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
421 PetscFunctionReturn(0);
422}
423
424/** @brief Verifies periodic Ucat endpoints are repaired before Cart2Contra terminal faces. */
426{
427 SimCtx *simCtx=NULL; UserCtx *user=NULL; Cmpnts ***ucat=NULL,***ucont=NULL; PetscInt mx;
428 const char *cell_fields[]={"Ucat"}, *staggered_fields[]={"Ucont"};
429 PetscFunctionBeginUser;
430 PetscCall(PicurvCreateMinimalContextsWithPeriodicity(&simCtx,&user,4,4,4,PETSC_TRUE,PETSC_FALSE,PETSC_FALSE));
433 mx=user->info.mx;
434 PetscCall(VecSet(user->Ucat,0.0));
435 PetscCall(DMDAVecGetArray(user->fda,user->Ucat,&ucat));
436 for(PetscInt k=1;k<user->info.mz-1;k++) for(PetscInt j=1;j<user->info.my-1;j++)
437 for(PetscInt i=1;i<mx-1;i++) ucat[k][j][i]=(Cmpnts){10.0+(PetscReal)i,2.0,3.0};
438 ucat[2][2][0]=(Cmpnts){-1000.0,-1000.0,-1000.0};
439 ucat[2][2][mx-1]=(Cmpnts){-2000.0,-2000.0,-2000.0};
440 PetscCall(DMDAVecRestoreArray(user->fda,user->Ucat,&ucat));
441 PetscCall(SynchronizePeriodicCellFields(user,1,cell_fields));
442 PetscCall(UpdateLocalGhosts(user,"Ucat"));
443 PetscCall(Cart2Contra(user));
444 PetscCall(SynchronizePeriodicStaggeredFields(user,1,staggered_fields));
445 PetscCall(DMDAVecGetArrayRead(user->fda,user->Ucont,&ucont));
446 PetscCall(PicurvAssertRealNear(0.5*((10.0+mx-2)+(10.0+1)),ucont[2][2][mx-2].x,1e-12,"terminal Xi face uses corrected periodic Ucat endpoint"));
447 PetscCall(DMDAVecRestoreArrayRead(user->fda,user->Ucont,&ucont));
448 PetscCall(PicurvDestroyMinimalContexts(&simCtx,&user));
449 PetscFunctionReturn(0);
450}
451
452/**
453 * @brief Tests loading a staged Ucat file IC and converting it to Ucont.
454 */
456{
457 SimCtx *simCtx = NULL;
458 UserCtx *user = NULL;
459 Cmpnts ***ucont = NULL;
460 char tmpdir[PETSC_MAX_PATH_LEN];
461
462 PetscFunctionBeginUser;
463 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
464 PetscCall(PicurvMakeTempDir(tmpdir, sizeof(tmpdir)));
465 PetscCall(VecSet(user->Ucat, 2.5));
466 PetscCall(PetscStrncpy(simCtx->_io_context_buffer, tmpdir, sizeof(simCtx->_io_context_buffer)));
468 PetscCall(WriteFieldData(user, "ufield", user->Ucat, 0, "dat"));
469 simCtx->current_io_directory = NULL;
470
471 PetscCall(VecZeroEntries(user->Ucat));
472 PetscCall(VecZeroEntries(user->Ucont));
475 PetscCall(PetscStrncpy(simCtx->initialConditionDirectory, tmpdir, sizeof(simCtx->initialConditionDirectory)));
476 PetscCall(PopulateInitialUcont(user));
477
478 PetscCall(PicurvAssertVecConstant(user->Ucat, 2.5, 1.0e-12, "file IC should restore Ucat"));
479 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
480 PetscCall(PicurvAssertRealNear(2.5, ucont[1][1][1].x, 1.0e-12, "file Ucat IC should populate interior Xi flux"));
481 PetscCall(PicurvAssertRealNear(2.5, ucont[1][1][1].y, 1.0e-12, "file Ucat IC should populate interior Eta flux"));
482 PetscCall(PicurvAssertRealNear(2.5, ucont[1][1][1].z, 1.0e-12, "file Ucat IC should populate interior Zeta flux"));
483 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
484 PetscCall(PicurvRemoveTempDir(tmpdir));
485 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
486 PetscFunctionReturn(0);
487}
488
489/**
490 * @brief Tests loading a staged Ucont file IC without Cartesian conversion.
491 */
493{
494 SimCtx *simCtx = NULL;
495 UserCtx *user = NULL;
496 char tmpdir[PETSC_MAX_PATH_LEN];
497
498 PetscFunctionBeginUser;
499 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
500 PetscCall(PicurvMakeTempDir(tmpdir, sizeof(tmpdir)));
501 PetscCall(VecSet(user->Ucont, 4.5));
502 PetscCall(PetscStrncpy(simCtx->_io_context_buffer, tmpdir, sizeof(simCtx->_io_context_buffer)));
504 PetscCall(WriteFieldData(user, "vfield", user->Ucont, 0, "dat"));
505 simCtx->current_io_directory = NULL;
506
507 PetscCall(VecZeroEntries(user->Ucont));
510 PetscCall(PetscStrncpy(simCtx->initialConditionDirectory, tmpdir, sizeof(simCtx->initialConditionDirectory)));
511 PetscCall(PopulateInitialUcont(user));
512
513 PetscCall(PicurvAssertVecConstant(user->Ucont, 4.5, 1.0e-12, "file Ucont IC should restore Ucont directly"));
514 PetscCall(PicurvRemoveTempDir(tmpdir));
515 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
516 PetscFunctionReturn(0);
517}
518/**
519 * @brief Tests direct interpolation from Eulerian fields to one localized swarm particle.
520 */
521
523{
524 SimCtx *simCtx = NULL;
525 UserCtx *user = NULL;
526 Cmpnts ***grad = NULL;
527 PetscReal *velocity = NULL;
528 PetscReal *diffusivity = NULL;
529 PetscReal *diffusivity_gradient = NULL;
530
531 PetscFunctionBeginUser;
532 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
533 PetscCall(PicurvCreateSwarmPair(user, 1, "ske"));
534 PetscCall(VecSet(user->Ucat, 2.0));
535 PetscCall(VecSet(user->Diffusivity, 0.25));
536
537 PetscCall(DMDAVecGetArray(user->fda, user->DiffusivityGradient, &grad));
538 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
539 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
540 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
541 grad[k][j][i].x = 0.1;
542 grad[k][j][i].y = 0.2;
543 grad[k][j][i].z = 0.3;
544 }
545 }
546 }
547 PetscCall(DMDAVecRestoreArray(user->fda, user->DiffusivityGradient, &grad));
548 PetscCall(SyncRuntimeFieldGhosts(user));
549 PetscCall(SeedSingleParticle(user, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, ACTIVE_AND_LOCATED));
550
551 PetscCall(InterpolateAllFieldsToSwarm(user));
552
553 PetscCall(DMSwarmGetField(user->swarm, "velocity", NULL, NULL, (void **)&velocity));
554 PetscCall(DMSwarmGetField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
555 PetscCall(DMSwarmGetField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
556 PetscCall(PicurvAssertRealNear(2.0, velocity[0], 1.0e-12, "Interpolated velocity x should match constant Eulerian field"));
557 PetscCall(PicurvAssertRealNear(2.0, velocity[1], 1.0e-12, "Interpolated velocity y should match constant Eulerian field"));
558 PetscCall(PicurvAssertRealNear(2.0, velocity[2], 1.0e-12, "Interpolated velocity z should match constant Eulerian field"));
559 PetscCall(PicurvAssertRealNear(0.25, diffusivity[0], 1.0e-12, "Interpolated scalar diffusivity should match constant Eulerian field"));
560 PetscCall(PicurvAssertRealNear(0.1, diffusivity_gradient[0], 1.0e-12, "Interpolated diffusivity-gradient x component"));
561 PetscCall(PicurvAssertRealNear(0.2, diffusivity_gradient[1], 1.0e-12, "Interpolated diffusivity-gradient y component"));
562 PetscCall(PicurvAssertRealNear(0.3, diffusivity_gradient[2], 1.0e-12, "Interpolated diffusivity-gradient z component"));
563 PetscCall(DMSwarmRestoreField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
564 PetscCall(DMSwarmRestoreField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
565 PetscCall(DMSwarmRestoreField(user->swarm, "velocity", NULL, NULL, (void **)&velocity));
566
567 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
568 PetscFunctionReturn(0);
569}
570/**
571 * @brief Tests the corner-averaged (legacy) interpolation path on constant fields.
572 */
573
575{
576 SimCtx *simCtx = NULL;
577 UserCtx *user = NULL;
578 Cmpnts ***grad = NULL;
579 PetscReal *velocity = NULL;
580 PetscReal *diffusivity = NULL;
581 PetscReal *diffusivity_gradient = NULL;
582
583 PetscFunctionBeginUser;
584 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
586 PetscCall(PicurvCreateSwarmPair(user, 1, "ske"));
587 PetscCall(VecSet(user->Ucat, 2.0));
588 PetscCall(VecSet(user->Diffusivity, 0.25));
589
590 PetscCall(DMDAVecGetArray(user->fda, user->DiffusivityGradient, &grad));
591 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
592 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
593 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
594 grad[k][j][i].x = 0.1;
595 grad[k][j][i].y = 0.2;
596 grad[k][j][i].z = 0.3;
597 }
598 }
599 }
600 PetscCall(DMDAVecRestoreArray(user->fda, user->DiffusivityGradient, &grad));
601 PetscCall(SyncRuntimeFieldGhosts(user));
602 PetscCall(SeedSingleParticle(user, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, ACTIVE_AND_LOCATED));
603
604 PetscCall(InterpolateAllFieldsToSwarm(user));
605
606 PetscCall(DMSwarmGetField(user->swarm, "velocity", NULL, NULL, (void **)&velocity));
607 PetscCall(DMSwarmGetField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
608 PetscCall(DMSwarmGetField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
609 PetscCall(PicurvAssertRealNear(2.0, velocity[0], 1.0e-12, "CornerAveraged: interpolated velocity x should match constant Eulerian field"));
610 PetscCall(PicurvAssertRealNear(2.0, velocity[1], 1.0e-12, "CornerAveraged: interpolated velocity y should match constant Eulerian field"));
611 PetscCall(PicurvAssertRealNear(2.0, velocity[2], 1.0e-12, "CornerAveraged: interpolated velocity z should match constant Eulerian field"));
612 PetscCall(PicurvAssertRealNear(0.25, diffusivity[0], 1.0e-12, "CornerAveraged: interpolated scalar diffusivity should match constant Eulerian field"));
613 PetscCall(PicurvAssertRealNear(0.1, diffusivity_gradient[0], 1.0e-12, "CornerAveraged: interpolated diffusivity-gradient x component"));
614 PetscCall(PicurvAssertRealNear(0.2, diffusivity_gradient[1], 1.0e-12, "CornerAveraged: interpolated diffusivity-gradient y component"));
615 PetscCall(PicurvAssertRealNear(0.3, diffusivity_gradient[2], 1.0e-12, "CornerAveraged: interpolated diffusivity-gradient z component"));
616 PetscCall(DMSwarmRestoreField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
617 PetscCall(DMSwarmRestoreField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
618 PetscCall(DMSwarmRestoreField(user->swarm, "velocity", NULL, NULL, (void **)&velocity));
619
620 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
621 PetscFunctionReturn(0);
622}
623/**
624 * @brief Tests particle-to-grid scattering using known cell occupancy and scalar values.
625 */
626
628{
629 SimCtx *simCtx = NULL;
630 UserCtx *user = NULL;
631 PetscInt *cell_ids = NULL;
632 PetscReal *psi = NULL;
633 PetscReal ***psi_grid = NULL;
634
635 PetscFunctionBeginUser;
636 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
637 PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
638
639 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
640 cell_ids[0] = 0; cell_ids[1] = 0; cell_ids[2] = 0;
641 cell_ids[3] = 0; cell_ids[4] = 0; cell_ids[5] = 0;
642 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
643
644 PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
645 psi[0] = 1.0;
646 psi[1] = 3.0;
647 PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
648
650
651 PetscCall(DMDAVecGetArrayRead(user->da, user->Psi, &psi_grid));
652 PetscCall(PicurvAssertRealNear(2.0, psi_grid[1][1][1], 1.0e-12, "Scatter should average particle Psi values into the owning cell"));
653 PetscCall(DMDAVecRestoreArrayRead(user->da, user->Psi, &psi_grid));
654
655 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
656 PetscFunctionReturn(0);
657}
658/**
659 * @brief Tests particle counting by geometric cell IDs using the production +1 storage shift.
660 */
661
663{
664 SimCtx *simCtx = NULL;
665 UserCtx *user = NULL;
666 PetscInt *cell_ids = NULL;
667 PetscReal ***counts = NULL;
668
669 PetscFunctionBeginUser;
670 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
671 PetscCall(PicurvCreateSwarmPair(user, 3, "ske"));
672
673 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
674 cell_ids[0] = 0; cell_ids[1] = 0; cell_ids[2] = 0;
675 cell_ids[3] = 0; cell_ids[4] = 0; cell_ids[5] = 0;
676 cell_ids[6] = 1; cell_ids[7] = 0; cell_ids[8] = 0;
677 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
678
679 PetscCall(CalculateParticleCountPerCell(user));
680
681 PetscCall(DMDAVecGetArrayRead(user->da, user->ParticleCount, &counts));
682 PetscCall(PicurvAssertRealNear(2.0, counts[1][1][1], 1.0e-12, "Two particles in cell (0,0,0) should accumulate at shifted index (1,1,1)"));
683 PetscCall(PicurvAssertRealNear(1.0, counts[1][1][2], 1.0e-12, "One particle in cell (1,0,0) should accumulate at shifted index (2,1,1)"));
684 PetscCall(DMDAVecRestoreArrayRead(user->da, user->ParticleCount, &counts));
685
686 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
687 PetscFunctionReturn(0);
688}
689/**
690 * @brief Tests localized particle-status reset behavior for restart of the location workflow.
691 */
692
694{
695 SimCtx *simCtx = NULL;
696 UserCtx *user = NULL;
697 PetscInt *status = NULL;
698
699 PetscFunctionBeginUser;
700 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
701 PetscCall(PicurvCreateSwarmPair(user, 3, "ske"));
702
703 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
704 status[0] = ACTIVE_AND_LOCATED;
705 status[1] = LOST;
706 status[2] = NEEDS_LOCATION;
707 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
708
709 PetscCall(ResetAllParticleStatuses(user));
710
711 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
712 PetscCall(PicurvAssertIntEqual(NEEDS_LOCATION, status[0], "ACTIVE_AND_LOCATED particles should be reset to NEEDS_LOCATION"));
713 PetscCall(PicurvAssertIntEqual(LOST, status[1], "LOST particles should remain LOST"));
714 PetscCall(PicurvAssertIntEqual(NEEDS_LOCATION, status[2], "NEEDS_LOCATION particles should remain unchanged"));
715 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
716
717 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
718 PetscFunctionReturn(0);
719}
720/**
721 * @brief Tests direct removal of particles that leave every rank bounding box.
722 */
723
725{
726 SimCtx *simCtx = NULL;
727 UserCtx *user = NULL;
728 PetscReal *positions = NULL;
729 PetscInt removed_local = 0;
730 PetscInt removed_global = 0;
731 PetscInt nlocal = 0;
732
733 PetscFunctionBeginUser;
734 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
735 PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
736
737 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
738 positions[0] = 0.5; positions[1] = 0.5; positions[2] = 0.5;
739 positions[3] = 9.0; positions[4] = 9.0; positions[5] = 9.0;
740 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
741
742 PetscCall(CheckAndRemoveOutOfBoundsParticles(user, &removed_local, &removed_global, simCtx->bboxlist));
743 PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
744 PetscCall(PicurvAssertIntEqual(1, removed_local, "Exactly one particle should be removed as out-of-bounds on a single rank"));
745 PetscCall(PicurvAssertIntEqual(1, removed_global, "Global out-of-bounds removal count should match the local single-rank result"));
746 PetscCall(PicurvAssertIntEqual(1, nlocal, "One in-bounds particle should remain after out-of-bounds removal"));
747
748 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
749 PetscFunctionReturn(0);
750}
751/**
752 * @brief Tests direct removal of particles already marked LOST by the location workflow.
753 */
754
756{
757 SimCtx *simCtx = NULL;
758 UserCtx *user = NULL;
759 PetscInt *status = NULL;
760 PetscInt removed_local = 0;
761 PetscInt removed_global = 0;
762 PetscInt nlocal = 0;
763
764 PetscFunctionBeginUser;
765 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
766 PetscCall(PicurvCreateSwarmPair(user, 3, "ske"));
767
768 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
769 status[0] = ACTIVE_AND_LOCATED;
770 status[1] = LOST;
771 status[2] = LOST;
772 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
773
774 PetscCall(CheckAndRemoveLostParticles(user, &removed_local, &removed_global));
775 PetscCall(DMSwarmGetLocalSize(user->swarm, &nlocal));
776 PetscCall(PicurvAssertIntEqual(2, removed_local, "Two LOST particles should be removed locally"));
777 PetscCall(PicurvAssertIntEqual(2, removed_global, "Global LOST-particle removal count should match the local single-rank result"));
778 PetscCall(PicurvAssertIntEqual(1, nlocal, "One non-LOST particle should remain"));
779
780 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
781 PetscFunctionReturn(0);
782}
783/**
784 * @brief Tests Brownian displacement generation against a duplicated seeded RNG stream.
785 */
786
788{
789 SimCtx *simCtx = NULL;
790 UserCtx *user = NULL;
791 char tmpdir[PETSC_MAX_PATH_LEN];
792 Cmpnts first;
793 Cmpnts second;
794
795 PetscFunctionBeginUser;
796 PetscCall(PicurvBuildTinyRuntimeContext(NULL, PETSC_FALSE, &simCtx, &user, tmpdir, sizeof(tmpdir)));
797 simCtx->dt = 0.25;
798 PetscCall(PicurvAssertBool((PetscBool)(simCtx->BrownianMotionRNG != NULL),
799 "runtime setup path should initialize the Brownian RNG"));
800
801 PetscCall(PetscRandomSetSeed(simCtx->BrownianMotionRNG, 12345));
802 PetscCall(PetscRandomSeed(simCtx->BrownianMotionRNG));
803
804 PetscCall(CalculateBrownianDisplacement(user, 0.5, &first));
805 PetscCall(PetscRandomSetSeed(simCtx->BrownianMotionRNG, 12345));
806 PetscCall(PetscRandomSeed(simCtx->BrownianMotionRNG));
807 PetscCall(CalculateBrownianDisplacement(user, 0.5, &second));
808
809 PetscCall(PicurvAssertRealNear(first.x, second.x, 1.0e-12, "Resetting the Brownian RNG seed should reproduce the x displacement"));
810 PetscCall(PicurvAssertRealNear(first.y, second.y, 1.0e-12, "Resetting the Brownian RNG seed should reproduce the y displacement"));
811 PetscCall(PicurvAssertRealNear(first.z, second.z, 1.0e-12, "Resetting the Brownian RNG seed should reproduce the z displacement"));
812
813 PetscCall(PicurvDestroyRuntimeContext(&simCtx));
814 PetscCall(PicurvRemoveTempDir(tmpdir));
815 PetscFunctionReturn(0);
816}
817/**
818 * @brief Tests swarm-wide particle position updates using the same transport path as the runtime loop.
819 */
821{
822 SimCtx *simCtx = NULL;
823 UserCtx *user = NULL;
824 PetscReal *positions = NULL;
825 PetscReal *velocities = NULL;
826 PetscReal *diffusivity = NULL;
827 Cmpnts *diffusivity_gradient = NULL;
828 PetscReal *psi = NULL;
829 PetscReal *weights = NULL;
830 PetscInt *cell_ids = NULL;
831 PetscInt *status = NULL;
832 PetscInt64 *pid = NULL;
833
834 PetscFunctionBeginUser;
835 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
836 PetscCall(PicurvCreateSwarmPair(user, 1, "ske"));
837 simCtx->dt = 0.25;
838
839 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
840 PetscCall(DMSwarmGetField(user->swarm, "velocity", NULL, NULL, (void **)&velocities));
841 PetscCall(DMSwarmGetField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
842 PetscCall(DMSwarmGetField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
843 PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
844 PetscCall(DMSwarmGetField(user->swarm, "weight", NULL, NULL, (void **)&weights));
845 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
846 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
847 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
848
849 positions[0] = 0.20; positions[1] = 0.30; positions[2] = 0.40;
850 velocities[0] = 0.40; velocities[1] = -0.20; velocities[2] = 0.10;
851 diffusivity[0] = 0.0;
852 diffusivity_gradient[0].x = 0.10;
853 diffusivity_gradient[0].y = 0.20;
854 diffusivity_gradient[0].z = -0.10;
855 psi[0] = 0.5;
856 weights[0] = 0.5; weights[1] = 0.5; weights[2] = 0.5;
857 cell_ids[0] = 0; cell_ids[1] = 0; cell_ids[2] = 0;
858 status[0] = ACTIVE_AND_LOCATED;
859 pid[0] = 7;
860
861 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
862 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
863 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
864 PetscCall(DMSwarmRestoreField(user->swarm, "weight", NULL, NULL, (void **)&weights));
865 PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
866 PetscCall(DMSwarmRestoreField(user->swarm, "DiffusivityGradient", NULL, NULL, (void **)&diffusivity_gradient));
867 PetscCall(DMSwarmRestoreField(user->swarm, "Diffusivity", NULL, NULL, (void **)&diffusivity));
868 PetscCall(DMSwarmRestoreField(user->swarm, "velocity", NULL, NULL, (void **)&velocities));
869 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
870
871 PetscCall(UpdateAllParticlePositions(user));
872
873 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
874 PetscCall(PicurvAssertRealNear(0.325, positions[0], 1.0e-12, "UpdateAllParticlePositions should advect x"));
875 PetscCall(PicurvAssertRealNear(0.300, positions[1], 1.0e-12, "UpdateAllParticlePositions should advect y"));
876 PetscCall(PicurvAssertRealNear(0.400, positions[2], 1.0e-12, "UpdateAllParticlePositions should advect z"));
877 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
878
879 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
880 PetscFunctionReturn(0);
881}
882/**
883 * @brief Tests the location orchestrator fast path when a particle already carries a valid prior cell.
884 */
886{
887 SimCtx *simCtx = NULL;
888 UserCtx *user = NULL;
889 PetscReal *positions = NULL;
890 PetscReal *weights = NULL;
891 PetscInt *cell_ids = NULL;
892 PetscInt *status = NULL;
893 PetscInt64 *pid = NULL;
894
895 PetscFunctionBeginUser;
896 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
897 PetscCall(PetscMalloc1(simCtx->size, &user->RankCellInfoMap));
898 PetscCall(GetOwnedCellRange(&user->info, 0, &user->RankCellInfoMap[0].xs_cell, &user->RankCellInfoMap[0].xm_cell));
899 PetscCall(GetOwnedCellRange(&user->info, 1, &user->RankCellInfoMap[0].ys_cell, &user->RankCellInfoMap[0].ym_cell));
900 PetscCall(GetOwnedCellRange(&user->info, 2, &user->RankCellInfoMap[0].zs_cell, &user->RankCellInfoMap[0].zm_cell));
901 PetscCall(PicurvCreateSwarmPair(user, 1, "ske"));
902
903 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
904 PetscCall(DMSwarmGetField(user->swarm, "weight", NULL, NULL, (void **)&weights));
905 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
906 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
907 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
908 positions[0] = 0.375; positions[1] = 0.375; positions[2] = 0.375;
909 weights[0] = 0.5; weights[1] = 0.5; weights[2] = 0.5;
910 cell_ids[0] = 1; cell_ids[1] = 1; cell_ids[2] = 1;
911 status[0] = NEEDS_LOCATION;
912 pid[0] = 11;
913 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
914 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
915 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
916 PetscCall(DMSwarmRestoreField(user->swarm, "weight", NULL, NULL, (void **)&weights));
917 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
918
919 PetscCall(LocateAllParticlesInGrid(user, simCtx->bboxlist));
920
921 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
922 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
923 PetscCall(PicurvAssertIntEqual(1, cell_ids[0], "prior-cell fast path should preserve the i cell id"));
924 PetscCall(PicurvAssertIntEqual(1, cell_ids[1], "prior-cell fast path should preserve the j cell id"));
925 PetscCall(PicurvAssertIntEqual(1, cell_ids[2], "prior-cell fast path should preserve the k cell id"));
926 PetscCall(PicurvAssertIntEqual(ACTIVE_AND_LOCATED, status[0], "prior-cell fast path should mark the particle ACTIVE_AND_LOCATED"));
927 PetscCall(PicurvAssertIntEqual(1, simCtx->searchMetrics.searchAttempts, "prior-cell fast path should record one search attempt"));
928 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.searchPopulation, "prior-cell fast path should record one input particle"));
929 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.searchLocatedCount, "prior-cell fast path should count one located particle"));
930 PetscCall(PicurvAssertIntEqual(0, (PetscInt)simCtx->searchMetrics.searchLostCount, "prior-cell fast path should not lose the particle"));
931 PetscCall(PicurvAssertIntEqual(0, (PetscInt)simCtx->searchMetrics.reSearchCount, "prior-cell fast path should not re-search on later passes"));
932 PetscCall(PicurvAssertBool((PetscBool)(simCtx->searchMetrics.traversalStepsSum > 0), "prior-cell fast path should accumulate traversal steps"));
933 PetscCall(PicurvAssertIntEqual(1, simCtx->searchMetrics.maxParticlePassDepth, "prior-cell fast path should report one settlement pass"));
934 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
935 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
936
937 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
938 PetscFunctionReturn(0);
939}
940/**
941 * @brief Tests the guess-then-verify orchestrator path for a local particle with an unknown prior cell.
942 */
944{
945 SimCtx *simCtx = NULL;
946 UserCtx *user = NULL;
947 PetscReal *positions = NULL;
948 PetscReal *weights = NULL;
949 PetscInt *cell_ids = NULL;
950 PetscInt *status = NULL;
951 PetscInt64 *pid = NULL;
952
953 PetscFunctionBeginUser;
954 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
955 PetscCall(PetscMalloc1(simCtx->size, &user->RankCellInfoMap));
956 PetscCall(GetOwnedCellRange(&user->info, 0, &user->RankCellInfoMap[0].xs_cell, &user->RankCellInfoMap[0].xm_cell));
957 PetscCall(GetOwnedCellRange(&user->info, 1, &user->RankCellInfoMap[0].ys_cell, &user->RankCellInfoMap[0].ym_cell));
958 PetscCall(GetOwnedCellRange(&user->info, 2, &user->RankCellInfoMap[0].zs_cell, &user->RankCellInfoMap[0].zm_cell));
959 PetscCall(PicurvCreateSwarmPair(user, 1, "ske"));
960
961 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
962 PetscCall(DMSwarmGetField(user->swarm, "weight", NULL, NULL, (void **)&weights));
963 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
964 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
965 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
966 positions[0] = 0.625; positions[1] = 0.625; positions[2] = 0.625;
967 weights[0] = 0.5; weights[1] = 0.5; weights[2] = 0.5;
968 cell_ids[0] = -1; cell_ids[1] = -1; cell_ids[2] = -1;
969 status[0] = NEEDS_LOCATION;
970 pid[0] = 22;
971 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_pid", NULL, NULL, (void **)&pid));
972 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
973 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
974 PetscCall(DMSwarmRestoreField(user->swarm, "weight", NULL, NULL, (void **)&weights));
975 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
976
977 PetscCall(LocateAllParticlesInGrid(user, simCtx->bboxlist));
978
979 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
980 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
981 PetscCall(PicurvAssertIntEqual(2, cell_ids[0], "guess-path location should resolve the i cell id"));
982 PetscCall(PicurvAssertIntEqual(2, cell_ids[1], "guess-path location should resolve the j cell id"));
983 PetscCall(PicurvAssertIntEqual(2, cell_ids[2], "guess-path location should resolve the k cell id"));
984 PetscCall(PicurvAssertIntEqual(ACTIVE_AND_LOCATED, status[0], "guess-path location should mark the particle ACTIVE_AND_LOCATED"));
985 PetscCall(PicurvAssertIntEqual(1, simCtx->searchMetrics.searchAttempts, "guess-path location should perform one robust search"));
986 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.searchPopulation, "guess-path location should count one input particle"));
987 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.searchLocatedCount, "guess-path location should count one located particle"));
988 PetscCall(PicurvAssertIntEqual(0, (PetscInt)simCtx->searchMetrics.searchLostCount, "guess-path location should not lose the particle"));
989 PetscCall(PicurvAssertIntEqual(0, (PetscInt)simCtx->searchMetrics.reSearchCount, "guess-path location should not count later-pass re-searches"));
990 PetscCall(PicurvAssertIntEqual(1, simCtx->searchMetrics.bboxGuessFallbackCount, "guess-path location should record one bbox fallback"));
991 PetscCall(PicurvAssertIntEqual(0, simCtx->searchMetrics.bboxGuessSuccessCount, "guess-path local resolution should not count as remote bbox success"));
992 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", NULL, NULL, (void **)&status));
993 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", NULL, NULL, (void **)&cell_ids));
994
995 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
996 PetscFunctionReturn(0);
997}
998/**
999 * @brief Verifies that later settlement passes increment re-search metrics.
1000 */
1002{
1003 SimCtx *simCtx = NULL;
1004 UserCtx *user = NULL;
1005 Particle particle;
1007
1008 PetscFunctionBeginUser;
1009 PetscCall(PetscMemzero(&particle, sizeof(particle)));
1010 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
1011 PetscCall(PetscMalloc1(simCtx->size, &user->RankCellInfoMap));
1012 PetscCall(GetOwnedCellRange(&user->info, 0, &user->RankCellInfoMap[0].xs_cell, &user->RankCellInfoMap[0].xm_cell));
1013 PetscCall(GetOwnedCellRange(&user->info, 1, &user->RankCellInfoMap[0].ys_cell, &user->RankCellInfoMap[0].ym_cell));
1014 PetscCall(GetOwnedCellRange(&user->info, 2, &user->RankCellInfoMap[0].zs_cell, &user->RankCellInfoMap[0].zm_cell));
1015
1017 particle.PID = 33;
1018 particle.cell[0] = 1;
1019 particle.cell[1] = 1;
1020 particle.cell[2] = 1;
1021 particle.loc.x = 0.375;
1022 particle.loc.y = 0.375;
1023 particle.loc.z = 0.375;
1024 particle.weights.x = 0.5;
1025 particle.weights.y = 0.5;
1026 particle.weights.z = 0.5;
1027
1028 PetscCall(LocateParticleOrFindMigrationTarget(user, &particle, &status));
1029
1030 PetscCall(PicurvAssertIntEqual(ACTIVE_AND_LOCATED, status, "direct re-search test should locate the particle"));
1031 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.searchAttempts, "direct re-search test should record one robust walk"));
1032 PetscCall(PicurvAssertIntEqual(1, (PetscInt)simCtx->searchMetrics.reSearchCount, "direct re-search test should increment re_search_count on later passes"));
1033 PetscCall(PicurvAssertIntEqual(0, (PetscInt)simCtx->searchMetrics.maxTraversalFailCount, "direct re-search test should not hit MAX_TRAVERSAL"));
1034
1035 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1036 PetscFunctionReturn(0);
1037}
1038/**
1039 * @brief Tests no-slip and free-slip wall helper kernels.
1040 */
1041
1042static PetscErrorCode TestWallNoSlipAndFreeSlipHelpers(void)
1043{
1044 Cmpnts wall_velocity = {0.0, 0.0, 0.0};
1045 Cmpnts reference_velocity = {2.0, 4.0, 6.0};
1046 Cmpnts boundary_velocity = {0.0, 0.0, 0.0};
1047 Cmpnts free_slip_reference = {2.0, 3.0, 4.0};
1048
1049 PetscFunctionBeginUser;
1050 noslip(NULL, 2.0, 1.0, wall_velocity, reference_velocity, &boundary_velocity, 1.0, 0.0, 0.0);
1051 PetscCall(PicurvAssertRealNear(1.0, boundary_velocity.x, 1.0e-12, "no-slip interpolated x"));
1052 PetscCall(PicurvAssertRealNear(2.0, boundary_velocity.y, 1.0e-12, "no-slip interpolated y"));
1053 PetscCall(PicurvAssertRealNear(3.0, boundary_velocity.z, 1.0e-12, "no-slip interpolated z"));
1054
1055 freeslip(NULL, 2.0, 1.0, wall_velocity, free_slip_reference, &boundary_velocity, 1.0, 0.0, 0.0);
1056 PetscCall(PicurvAssertRealNear(1.0, boundary_velocity.x, 1.0e-12, "free-slip interpolated normal component"));
1057 PetscCall(PicurvAssertRealNear(3.0, boundary_velocity.y, 1.0e-12, "free-slip tangential y preserved"));
1058 PetscCall(PicurvAssertRealNear(4.0, boundary_velocity.z, 1.0e-12, "free-slip tangential z preserved"));
1059 PetscFunctionReturn(0);
1060}
1061/**
1062 * @brief Tests wall-model scalar helper kernels.
1063 */
1064
1065static PetscErrorCode TestWallModelScalarHelpers(void)
1066{
1067 const PetscReal expected_smooth_e = PetscExpReal(0.41 * 5.5);
1068 PetscReal e_coeff = 0.0;
1069 PetscReal utau = 0.0;
1070 PetscReal residual = 0.0;
1071
1072 PetscFunctionBeginUser;
1073 e_coeff = E_coeff(0.1, 0.0, 1.0e-3);
1074 PetscCall(PicurvAssertRealNear(expected_smooth_e, e_coeff, 1.0e-10, "smooth-wall E coefficient"));
1075
1076 utau = find_utau_hydset(1.0e-3, 1.0, 1.0e-2, 0.1, 0.0);
1077 PetscCall(PicurvAssertBool((PetscBool)(utau > 0.0), "friction velocity should remain positive"));
1078 residual = f_hydset(1.0e-3, 1.0, 1.0e-2, utau, 0.0);
1079 PetscCall(PicurvAssertBool((PetscBool)(PetscAbsReal(residual) < 1.0e-5), "Newton solve residual should be small"));
1080
1081 PetscCall(PicurvAssertRealNear(0.0, nu_t(0.0), 1.0e-12, "eddy viscosity ratio at wall"));
1082 PetscCall(PicurvAssertBool((PetscBool)(integrate_1(1.0e-3, 1.0e-2, 0.1, 0) > 0.0), "integral helper should be positive"));
1083 PetscFunctionReturn(0);
1084}
1085/**
1086 * @brief Tests closed-form and iterative wall-model velocity helpers against inverse reconstructions.
1087 */
1088
1089static PetscErrorCode TestWallModelVelocityHelpers(void)
1090{
1091 const PetscReal kinematic_viscosity = 1.0e-3;
1092 const PetscReal wall_distance = 2.0e-2;
1093 const PetscReal target_velocity = 1.0;
1094 const PetscReal roughness_length = 1.0e-4;
1095 PetscReal utau_loglaw = 0.0;
1096 PetscReal utau_werner = 0.0;
1097 PetscReal utau_cabot = 0.0;
1098 PetscReal wall_shear_velocity = 0.0;
1099 PetscReal wall_shear_normal = 0.0;
1100
1101 PetscFunctionBeginUser;
1102 utau_loglaw = find_utau_loglaw(target_velocity, wall_distance, roughness_length);
1103 PetscCall(PicurvAssertRealNear(target_velocity, u_loglaw(wall_distance, utau_loglaw, roughness_length), 1.0e-12,
1104 "simple log-law inversion should reconstruct the target velocity"));
1105
1106 utau_werner = find_utau_Werner(kinematic_viscosity, target_velocity, wall_distance, 0.1);
1107 PetscCall(PicurvAssertBool((PetscBool)(utau_werner > 0.0), "Werner-Wengle friction velocity should remain positive"));
1108 PetscCall(PicurvAssertRealNear(target_velocity, u_Werner(kinematic_viscosity, wall_distance, utau_werner), 1.0e-6,
1109 "Werner-Wengle inversion should reconstruct the target velocity"));
1110
1111 find_utau_Cabot(kinematic_viscosity, target_velocity, wall_distance, 0.1, 0.0, 0.0,
1112 &utau_cabot, &wall_shear_velocity, &wall_shear_normal);
1113 PetscCall(PicurvAssertBool((PetscBool)(utau_cabot > 0.0), "Cabot friction velocity should remain positive"));
1114 PetscCall(PicurvAssertRealNear(target_velocity, u_Cabot(kinematic_viscosity, wall_distance, utau_cabot, 0.0, wall_shear_velocity), 1.0e-6,
1115 "Cabot inversion should reconstruct the target velocity when pressure gradient is zero"));
1116 PetscCall(PicurvAssertRealNear(0.0, wall_shear_normal, 1.0e-10,
1117 "zero normal pressure gradient should keep Cabot normal wall shear at zero"));
1118 PetscFunctionReturn(0);
1119}
1120/**
1121 * @brief Tests the vector wall-function wrappers on a tangential reference flow.
1122 */
1123static PetscErrorCode TestWallFunctionVectorWrappers(void)
1124{
1125 SimCtx *simCtx = NULL;
1126 UserCtx *user = NULL;
1127 Cmpnts wall_velocity = {0.0, 0.0, 0.0};
1128 Cmpnts reference_velocity = {0.0, 1.0, 0.0};
1129 Cmpnts boundary_velocity = {0.0, 0.0, 0.0};
1130 PetscReal friction_velocity = 0.0;
1131
1132 PetscFunctionBeginUser;
1133 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
1134 simCtx->ren = 1000.0;
1135
1136 wall_function(user, 2.0e-2, 1.0e-2, wall_velocity, reference_velocity, &boundary_velocity, &friction_velocity, 1.0, 0.0, 0.0);
1137 PetscCall(PicurvAssertRealNear(0.0, boundary_velocity.x, 1.0e-12, "Werner wall function should preserve zero normal velocity"));
1138 PetscCall(PicurvAssertBool((PetscBool)(boundary_velocity.y > 0.0 && boundary_velocity.y < 1.0), "Werner wall function should damp tangential velocity"));
1139 PetscCall(PicurvAssertBool((PetscBool)(friction_velocity > 0.0), "Werner wall function should compute positive friction velocity"));
1140
1141 wall_function_loglaw(user, 1.0e-4, 2.0e-2, 1.0e-2, wall_velocity, reference_velocity, &boundary_velocity, &friction_velocity, 1.0, 0.0, 0.0);
1142 PetscCall(PicurvAssertRealNear(0.0, boundary_velocity.x, 1.0e-12, "log-law wall function should preserve zero normal velocity"));
1143 PetscCall(PicurvAssertBool((PetscBool)(boundary_velocity.y > 0.0 && boundary_velocity.y <= 1.0), "log-law wall function should keep tangential velocity bounded"));
1144 PetscCall(PicurvAssertBool((PetscBool)(friction_velocity > 0.0), "log-law wall function should compute positive friction velocity"));
1145
1146 wall_function_Cabot(user, 1.0e-4, 2.0e-2, 1.0e-2, wall_velocity, reference_velocity, &boundary_velocity, &friction_velocity,
1147 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10);
1148 PetscCall(PicurvAssertRealNear(0.0, boundary_velocity.x, 1.0e-12, "Cabot wall function should preserve zero normal velocity"));
1149 PetscCall(PicurvAssertBool((PetscBool)(boundary_velocity.y > 0.0 && boundary_velocity.y <= 1.0), "Cabot wall function should keep tangential velocity bounded"));
1150 PetscCall(PicurvAssertBool((PetscBool)(friction_velocity > 0.0), "Cabot wall function should compute positive friction velocity"));
1151
1152 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1153 PetscFunctionReturn(0);
1154}
1155/**
1156 * @brief Tests driven-flow validation when no driven handlers are present.
1157 */
1158
1160{
1161 SimCtx *simCtx = NULL;
1162 UserCtx *user = NULL;
1163
1164 PetscFunctionBeginUser;
1165 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
1166 PetscCall(Validate_DrivenFlowConfiguration(user));
1167 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1168 PetscFunctionReturn(0);
1169}
1170/**
1171 * @brief Tests the constant Smagorinsky model helper path.
1172 */
1173
1175{
1176 SimCtx *simCtx = NULL;
1177 UserCtx *user = NULL;
1178 PetscReal ***lcs = NULL;
1179
1180 PetscFunctionBeginUser;
1181 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
1182 PetscCall(DMCreateGlobalVector(user->da, &user->CS));
1183 PetscCall(DMCreateLocalVector(user->da, &user->lCs));
1184 simCtx->step = 2;
1185 simCtx->StartStep = 0;
1186 simCtx->les = CONSTANT_SMAGORINSKY;
1187 simCtx->Const_CS = 0.17;
1188
1189 PetscCall(ComputeSmagorinskyConstant(user));
1190 PetscCall(PicurvAssertVecConstant(user->CS, 0.17, 1.0e-12, "constant Smagorinsky branch should fill CS"));
1191 PetscCall(DMDAVecGetArrayRead(user->da, user->lCs, &lcs));
1192 PetscCall(PicurvAssertRealNear(0.17, lcs[2][2][2], 1.0e-12,
1193 "constant Smagorinsky branch should refresh local CS"));
1194 PetscCall(DMDAVecRestoreArrayRead(user->da, user->lCs, &lcs));
1195
1196 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1197 PetscFunctionReturn(0);
1198}
1199/**
1200 * @brief Tests that the shared minimal fixture mirrors the production DA contract.
1201 */
1202
1204{
1205 SimCtx *simCtx = NULL;
1206 UserCtx *user = NULL;
1207 DM coord_dm = NULL;
1208 PetscInt mx = 0, my = 0, mz = 0;
1209
1210 PetscFunctionBeginUser;
1211 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 6, 4));
1212
1213 PetscCall(DMDAGetInfo(user->da, NULL, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1214 PetscCall(PicurvAssertIntEqual(user->IM + 1, mx, "minimal fixture should size da with IM+1 nodes"));
1215 PetscCall(PicurvAssertIntEqual(user->JM + 1, my, "minimal fixture should size da with JM+1 nodes"));
1216 PetscCall(PicurvAssertIntEqual(user->KM + 1, mz, "minimal fixture should size da with KM+1 nodes"));
1217 PetscCall(PicurvAssertIntEqual(mx, user->info.mx, "user->info should be sourced from the production da"));
1218 PetscCall(PicurvAssertIntEqual(my, user->info.my, "user->info my should match the da dimensions"));
1219 PetscCall(PicurvAssertIntEqual(mz, user->info.mz, "user->info mz should match the da dimensions"));
1220
1221 PetscCall(DMGetCoordinateDM(user->da, &coord_dm));
1222 PetscCall(PicurvAssertBool((PetscBool)(coord_dm == user->fda),
1223 "minimal fixture should derive fda from the coordinate-DM path"));
1224
1225 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1226 PetscFunctionReturn(0);
1227}
1228/**
1229 * @brief Tests that the shared swarm fixture registers the production field set.
1230 */
1231
1233{
1234 SimCtx *simCtx = NULL;
1235 UserCtx *user = NULL;
1236 PetscInt bs = 0;
1237 void *field_ptr = NULL;
1238
1239 PetscFunctionBeginUser;
1240 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
1241 PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
1242
1243 PetscCall(DMSwarmGetField(user->swarm, "position", &bs, NULL, &field_ptr));
1244 PetscCall(PicurvAssertIntEqual(3, bs, "solver swarm should register particle position"));
1245 PetscCall(PicurvAssertBool((PetscBool)(field_ptr != NULL), "position field should be retrievable"));
1246 PetscCall(DMSwarmRestoreField(user->swarm, "position", &bs, NULL, &field_ptr));
1247
1248 PetscCall(DMSwarmGetField(user->swarm, "velocity", &bs, NULL, &field_ptr));
1249 PetscCall(PicurvAssertIntEqual(3, bs, "solver swarm should register particle velocity"));
1250 PetscCall(DMSwarmRestoreField(user->swarm, "velocity", &bs, NULL, &field_ptr));
1251
1252 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_CellID", &bs, NULL, &field_ptr));
1253 PetscCall(PicurvAssertIntEqual(3, bs, "solver swarm should register DMSwarm_CellID"));
1254 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_CellID", &bs, NULL, &field_ptr));
1255
1256 PetscCall(DMSwarmGetField(user->swarm, "weight", &bs, NULL, &field_ptr));
1257 PetscCall(PicurvAssertIntEqual(3, bs, "solver swarm should register particle weight"));
1258 PetscCall(DMSwarmRestoreField(user->swarm, "weight", &bs, NULL, &field_ptr));
1259
1260 PetscCall(DMSwarmGetField(user->swarm, "Diffusivity", &bs, NULL, &field_ptr));
1261 PetscCall(PicurvAssertIntEqual(1, bs, "solver swarm should register particle diffusivity"));
1262 PetscCall(DMSwarmRestoreField(user->swarm, "Diffusivity", &bs, NULL, &field_ptr));
1263
1264 PetscCall(DMSwarmGetField(user->swarm, "DiffusivityGradient", &bs, NULL, &field_ptr));
1265 PetscCall(PicurvAssertIntEqual(3, bs, "solver swarm should register particle diffusivity gradients"));
1266 PetscCall(DMSwarmRestoreField(user->swarm, "DiffusivityGradient", &bs, NULL, &field_ptr));
1267
1268 PetscCall(DMSwarmGetField(user->swarm, "Psi", &bs, NULL, &field_ptr));
1269 PetscCall(PicurvAssertIntEqual(1, bs, "solver swarm should register particle scalar Psi"));
1270 PetscCall(DMSwarmRestoreField(user->swarm, "Psi", &bs, NULL, &field_ptr));
1271
1272 PetscCall(DMSwarmGetField(user->swarm, "DMSwarm_location_status", &bs, NULL, &field_ptr));
1273 PetscCall(PicurvAssertIntEqual(1, bs, "solver swarm should register particle location status"));
1274 PetscCall(DMSwarmRestoreField(user->swarm, "DMSwarm_location_status", &bs, NULL, &field_ptr));
1275
1276 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1277 PetscFunctionReturn(0);
1278}
1279/**
1280 * @brief Tests solver history-vector shifting between time levels.
1281 */
1282
1284{
1285 SimCtx *simCtx = NULL;
1286 UserCtx *user = NULL;
1287
1288 PetscFunctionBeginUser;
1289 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
1290
1291 PetscCall(VecSet(user->Ucont, 11.0));
1292 PetscCall(VecSet(user->Ucont_o, 7.0));
1293 PetscCall(VecSet(user->Ucont_rm1, 3.0));
1294 PetscCall(VecSet(user->Ucat, 5.0));
1295 PetscCall(VecSet(user->Ucat_o, -1.0));
1296 PetscCall(VecSet(user->P, 9.0));
1297 PetscCall(VecSet(user->P_o, -2.0));
1298
1299 PetscCall(UpdateSolverHistoryVectors(user));
1300
1301 PetscCall(PicurvAssertVecConstant(user->Ucont_o, 11.0, 1.0e-12, "Ucont_o should receive current Ucont"));
1302 PetscCall(PicurvAssertVecConstant(user->Ucont_rm1, 7.0, 1.0e-12, "Ucont_rm1 should receive prior Ucont_o"));
1303 PetscCall(PicurvAssertVecConstant(user->Ucat_o, 5.0, 1.0e-12, "Ucat_o should receive current Ucat"));
1304 PetscCall(PicurvAssertVecConstant(user->P_o, 9.0, 1.0e-12, "P_o should receive current P"));
1305 PetscCall(PicurvAssertVecConstant(user->lUcont_o, 11.0, 1.0e-12, "lUcont_o ghost sync should match Ucont_o"));
1306 PetscCall(PicurvAssertVecConstant(user->lUcont_rm1, 7.0, 1.0e-12, "lUcont_rm1 ghost sync should match Ucont_rm1"));
1307
1308 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1309 PetscFunctionReturn(0);
1310}
1311/**
1312 * @brief Tests owned-cell range accounting on a single MPI rank.
1313 */
1314
1316{
1317 SimCtx *simCtx = NULL;
1318 UserCtx *user = NULL;
1319 PetscInt xs = -1, ys = -1, zs = -1;
1320 PetscInt xm = -1, ym = -1, zm = -1;
1321
1322 PetscFunctionBeginUser;
1323 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 6, 4));
1324
1325 PetscCall(GetOwnedCellRange(&user->info, 0, &xs, &xm));
1326 PetscCall(GetOwnedCellRange(&user->info, 1, &ys, &ym));
1327 PetscCall(GetOwnedCellRange(&user->info, 2, &zs, &zm));
1328
1329 PetscCall(PicurvAssertIntEqual(0, xs, "single-rank x cell-start index"));
1330 PetscCall(PicurvAssertIntEqual(0, ys, "single-rank y cell-start index"));
1331 PetscCall(PicurvAssertIntEqual(0, zs, "single-rank z cell-start index"));
1332 PetscCall(PicurvAssertIntEqual(user->info.mx - 2, xm, "single-rank x owned cell count"));
1333 PetscCall(PicurvAssertIntEqual(user->info.my - 2, ym, "single-rank y owned cell count"));
1334 PetscCall(PicurvAssertIntEqual(user->info.mz - 2, zm, "single-rank z owned cell count"));
1335
1336 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1337 PetscFunctionReturn(0);
1338}
1339/**
1340 * @brief Tests neighbor-rank discovery on a single MPI rank.
1341 */
1342
1344{
1345 SimCtx *simCtx = NULL;
1346 UserCtx *user = NULL;
1347
1348 PetscFunctionBeginUser;
1349 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 6, 6, 6));
1350 PetscCall(ComputeAndStoreNeighborRanks(user));
1351
1352 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_xm, "single-rank xm neighbor should be null"));
1353 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_xp, "single-rank xp neighbor should be null"));
1354 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_ym, "single-rank ym neighbor should be null"));
1355 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_yp, "single-rank yp neighbor should be null"));
1356 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_zm, "single-rank zm neighbor should be null"));
1357 PetscCall(PicurvAssertIntEqual(MPI_PROC_NULL, user->neighbors.rank_zp, "single-rank zp neighbor should be null"));
1358
1359 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
1360 PetscFunctionReturn(0);
1361}
1362/**
1363 * @brief Tests parsing of positive runtime walltime metadata values.
1364 */
1365
1367{
1368 PetscReal seconds = 0.0;
1369
1370 PetscFunctionBeginUser;
1371 PetscCall(PicurvAssertBool(RuntimeWalltimeGuardParsePositiveSeconds("300", &seconds), "plain positive seconds should parse"));
1372 PetscCall(PicurvAssertRealNear(300.0, seconds, 1.0e-12, "parsed walltime seconds"));
1373 PetscCall(PicurvAssertBool(RuntimeWalltimeGuardParsePositiveSeconds(" 42.5 ", &seconds), "whitespace-wrapped decimal seconds should parse"));
1374 PetscCall(PicurvAssertRealNear(42.5, seconds, 1.0e-12, "parsed decimal walltime seconds"));
1375 PetscCall(PicurvAssertBool((PetscBool)!RuntimeWalltimeGuardParsePositiveSeconds("nope", &seconds), "non-numeric metadata should fail parsing"));
1376 PetscCall(PicurvAssertBool((PetscBool)!RuntimeWalltimeGuardParsePositiveSeconds("-10", &seconds), "negative metadata should fail parsing"));
1377 PetscFunctionReturn(0);
1378}
1379/**
1380 * @brief Tests walltime-guard estimator helper calculations.
1381 */
1382
1384{
1385 PetscReal ewma_fast = 0.0;
1386 PetscReal ewma_slow = 0.0;
1387 PetscReal conservative_fast = 0.0;
1388 PetscReal conservative_slow = 0.0;
1389 PetscReal required_headroom = 0.0;
1390
1391 PetscFunctionBeginUser;
1392 ewma_fast = RuntimeWalltimeGuardUpdateEWMA(PETSC_TRUE, 4.0, 6.0, 0.5);
1393 ewma_slow = RuntimeWalltimeGuardUpdateEWMA(PETSC_TRUE, ewma_fast, 12.0, 0.5);
1394 conservative_fast = RuntimeWalltimeGuardConservativeEstimate(5.0, ewma_fast, 6.0);
1395 conservative_slow = RuntimeWalltimeGuardConservativeEstimate(5.0, ewma_slow, 12.0);
1396 required_headroom = RuntimeWalltimeGuardRequiredHeadroom(8.0, 2.0, conservative_slow);
1397
1398 PetscCall(PicurvAssertRealNear(5.0, ewma_fast, 1.0e-12, "EWMA after moderate step"));
1399 PetscCall(PicurvAssertRealNear(8.5, ewma_slow, 1.0e-12, "EWMA after newer slow step"));
1400 PetscCall(PicurvAssertRealNear(6.0, conservative_fast, 1.0e-12, "conservative estimate tracks latest moderate step"));
1401 PetscCall(PicurvAssertRealNear(12.0, conservative_slow, 1.0e-12, "conservative estimate tracks newest slow step"));
1402 PetscCall(PicurvAssertRealNear(24.0, required_headroom, 1.0e-12, "required headroom scales with conservative estimate"));
1403 PetscFunctionReturn(0);
1404}
1405/**
1406 * @brief Tests runtime walltime-guard shutdown trigger decisions.
1407 */
1408
1410{
1411 PetscBool should_trigger = PETSC_FALSE;
1412 PetscReal required_headroom = 0.0;
1413
1414 PetscFunctionBeginUser;
1415 should_trigger = RuntimeWalltimeGuardShouldTrigger(9, 10, 15.0, 5.0, 2.0, 6.0, 6.0, 6.0, &required_headroom);
1416 PetscCall(PicurvAssertBool((PetscBool)!should_trigger, "guard should not trigger before warmup completes"));
1417
1418 should_trigger = RuntimeWalltimeGuardShouldTrigger(10, 10, 40.0, 5.0, 2.0, 10.0, 12.0, 14.0, &required_headroom);
1419 PetscCall(PicurvAssertBool((PetscBool)!should_trigger, "guard should not trigger when remaining walltime exceeds required headroom"));
1420 PetscCall(PicurvAssertRealNear(28.0, required_headroom, 1.0e-12, "required headroom after warmup"));
1421
1422 should_trigger = RuntimeWalltimeGuardShouldTrigger(10, 10, 28.0, 5.0, 2.0, 10.0, 12.0, 14.0, &required_headroom);
1423 PetscCall(PicurvAssertBool(should_trigger, "guard should trigger when remaining walltime reaches required headroom"));
1424 PetscCall(PicurvAssertRealNear(28.0, required_headroom, 1.0e-12, "required headroom remains unchanged at trigger threshold"));
1425 PetscFunctionReturn(0);
1426}
1427/**
1428 * @brief Runs the unit-runtime PETSc test binary.
1429 */
1430
1431int main(int argc, char **argv)
1432{
1433 PetscErrorCode ierr;
1434 const PicurvTestCase cases[] = {
1435 {"distribute-particles-remainder-handling", TestDistributeParticlesRemainderHandling},
1436 {"is-particle-inside-bbox-basic-cases", TestIsParticleInsideBoundingBoxBasicCases},
1437 {"update-particle-weights-computes-expected-ratios", TestUpdateParticleWeightsComputesExpectedRatios},
1438 {"update-particle-position-without-brownian-contribution", TestUpdateParticlePositionWithoutBrownianContribution},
1439 {"update-particle-position-diffusivity-gradient-only", TestUpdateParticlePositionDiffusivityGradientOnly},
1440 {"update-particle-field-iem-relaxation", TestUpdateParticleFieldIEMRelaxation},
1441 {"set-initial-interior-field-ignores-non-ucont-request", TestSetInitialInteriorFieldIgnoresNonUcontRequest},
1442 {"set-initial-interior-field-cartesian-constant-sets-contravariant-flux", TestSetInitialInteriorFieldCartesianConstantSetsContravariantFlux},
1443 {"set-initial-interior-field-curvilinear-constant-via-flow-direction", TestSetInitialInteriorFieldCurvilinearConstantViaFlowDirection},
1444 {"set-initial-interior-field-zero-clears-interior", TestSetInitialInteriorFieldZeroClearsInterior},
1445 {"set-initial-interior-field-poiseuille-profile", TestSetInitialInteriorFieldPoiseuilleProfile},
1446 {"cart2contra-converts-cartesian-field", TestCart2ContraConvertsCartesianField},
1447 {"cart2contra-uses-finalized-periodic-ucat", TestCart2ContraUsesFinalizedPeriodicUcat},
1448 {"populate-initial-ucont-loads-staged-ucat", TestPopulateInitialUcontLoadsStagedUcat},
1449 {"populate-initial-ucont-loads-staged-ucont", TestPopulateInitialUcontLoadsStagedUcont},
1450 {"interpolate-all-fields-to-swarm-constant-fields", TestInterpolateAllFieldsToSwarmConstantFields},
1451 {"interpolate-all-fields-to-swarm-corner-averaged-constant-fields", TestInterpolateAllFieldsToSwarmCornerAveragedConstantFields},
1452 {"scatter-all-particle-fields-to-euler-fields-averages-psi", TestScatterAllParticleFieldsToEulerFieldsAveragesPsi},
1453 {"calculate-particle-count-per-cell-counts-global-cell-ids", TestCalculateParticleCountPerCellCountsGlobalCellIDs},
1454 {"reset-all-particle-statuses-leaves-lost-particles-untouched", TestResetAllParticleStatusesLeavesLostParticlesUntouched},
1455 {"check-and-remove-out-of-bounds-particles-removes-escaped-particle", TestCheckAndRemoveOutOfBoundsParticlesRemovesEscapedParticle},
1456 {"check-and-remove-lost-particles-removes-lost-entries", TestCheckAndRemoveLostParticlesRemovesLostEntries},
1457 {"calculate-brownian-displacement-deterministic-seed", TestCalculateBrownianDisplacementDeterministicSeed},
1458 {"update-all-particle-positions-moves-swarm-entries", TestUpdateAllParticlePositionsMovesSwarmEntries},
1459 {"locate-all-particles-in-grid-prior-cell-fast-path", TestLocateAllParticlesInGridPriorCellFastPath},
1460 {"locate-all-particles-in-grid-guess-path-resolves-local-particle", TestLocateAllParticlesInGridGuessPathResolvesLocalParticle},
1461 {"locate-particle-or-find-migration-target-counts-research", TestLocateParticleOrFindMigrationTargetCountsReSearch},
1462 {"wall-noslip-and-freeslip-helpers", TestWallNoSlipAndFreeSlipHelpers},
1463 {"wall-model-scalar-helpers", TestWallModelScalarHelpers},
1464 {"wall-model-velocity-helpers", TestWallModelVelocityHelpers},
1465 {"wall-function-vector-wrappers", TestWallFunctionVectorWrappers},
1466 {"validate-driven-flow-configuration-no-driven-handlers", TestValidateDrivenFlowConfigurationNoDrivenHandlers},
1467 {"compute-smagorinsky-constant-constant-model", TestComputeSmagorinskyConstantConstantModel},
1468 {"minimal-fixture-mirrors-production-dm-layout", TestMinimalFixtureMirrorsProductionDMLayout},
1469 {"minimal-fixture-registers-production-swarm-fields", TestMinimalFixtureRegistersProductionSwarmFields},
1470 {"update-solver-history-vectors-shifts-states", TestUpdateSolverHistoryVectorsShiftsStates},
1471 {"get-owned-cell-range-single-rank-accounting", TestGetOwnedCellRangeSingleRankAccounting},
1472 {"compute-and-store-neighbor-ranks-single-rank", TestComputeAndStoreNeighborRanksSingleRank},
1473 {"runtime-walltime-guard-parses-positive-seconds", TestRuntimeWalltimeGuardParsesPositiveSeconds},
1474 {"runtime-walltime-guard-estimator-helpers", TestRuntimeWalltimeGuardEstimatorHelpers},
1475 {"runtime-walltime-guard-trigger-decision", TestRuntimeWalltimeGuardTriggerDecision},
1476 };
1477
1478 ierr = PetscInitialize(&argc, &argv, NULL, "PICurv runtime-kernel tests");
1479 if (ierr) {
1480 return (int)ierr;
1481 }
1482
1483 ierr = PicurvRunTests("unit-runtime", cases, sizeof(cases) / sizeof(cases[0]));
1484 if (ierr) {
1485 PetscFinalize();
1486 return (int)ierr;
1487 }
1488
1489 ierr = PetscFinalize();
1490 return (int)ierr;
1491}
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.
Definition BC_Handlers.c:15
PetscErrorCode SynchronizePeriodicStaggeredFields(UserCtx *user, PetscInt num_fields, const char *field_names[])
Synchronizes persistent component-staggered vector fields.
PetscErrorCode SynchronizePeriodicCellFields(UserCtx *user, PetscInt num_fields, const char *field_names[])
Synchronizes periodic endpoint cells for a list of cell-centered fields.
Header file for Particle Motion and migration related functions.
PetscErrorCode CheckAndRemoveOutOfBoundsParticles(UserCtx *user, PetscInt *removedCountLocal, PetscInt *removedCountGlobal, const BoundingBox *bboxlist)
Checks for particles outside the physical domain boundaries and removes them using DMSwarmRemovePoint...
PetscErrorCode UpdateAllParticlePositions(UserCtx *user)
Loops over all local particles in the DMSwarm, updating their positions based on velocity and the glo...
PetscErrorCode CalculateParticleCountPerCell(UserCtx *user)
Counts particles in each cell of the DMDA 'da' and stores the result in user->ParticleCount.
PetscErrorCode CalculateBrownianDisplacement(UserCtx *user, PetscReal diff_eff, Cmpnts *displacement)
Calculates the stochastic displacement vector (Brownian motion) for a single particle.
PetscErrorCode LocateAllParticlesInGrid(UserCtx *user, BoundingBox *bboxlist)
Orchestrates the complete particle location and migration process for one timestep.
PetscErrorCode UpdateParticlePosition(UserCtx *user, Particle *particle)
Updates a particle's position based on its velocity and the timestep dt (stored in user->dt).
PetscErrorCode ResetAllParticleStatuses(UserCtx *user)
Marks all local particles as NEEDS_LOCATION for the next settlement pass.
PetscErrorCode CheckAndRemoveLostParticles(UserCtx *user, PetscInt *removedCountLocal, PetscInt *removedCountGlobal)
Removes particles that have been definitively flagged as LOST by the location algorithm.
Header file for Particle related physics modules.
PetscErrorCode UpdateParticleField(const char *fieldName, PetscReal dt, PetscReal *psi_io, PetscReal diffusivity, PetscReal mean_val, PetscReal cell_vol, PetscReal C_model)
Updates a single particle's field based on its state and physics model.
Header file for Particle Swarm management functions.
PetscErrorCode UpdateParticleWeights(PetscReal *d, Particle *particle)
Updates a particle's interpolation weights based on distances to cell faces.
PetscErrorCode DistributeParticles(PetscInt numParticles, PetscMPIInt rank, PetscMPIInt size, PetscInt *particlesPerProcess, PetscInt *remainder)
Distributes particles evenly across MPI processes, handling any remainders.
PetscBool IsParticleInsideBoundingBox(const BoundingBox *bbox, const Particle *particle)
Checks if a particle's location is within a specified bounding box.
PetscErrorCode ScatterAllParticleFieldsToEulerFields(UserCtx *user)
Scatters a predefined set of particle fields to their corresponding Eulerian fields.
PetscErrorCode SetInitialInteriorField(UserCtx *user, const char *fieldName)
Sets the initial values for the INTERIOR of a specified Eulerian field.
PetscErrorCode PopulateInitialUcont(UserCtx *user)
Populate Ucont for one fresh-start block from the configured IC mode.
PetscErrorCode InterpolateAllFieldsToSwarm(UserCtx *user)
Interpolates all relevant fields from the DMDA to the DMSwarm.
PetscErrorCode WriteFieldData(UserCtx *user, const char *field_name, Vec field_vec, PetscInt ti, const char *ext)
Writes data from a specific PETSc vector to a file.
Definition io.c:1443
PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
Computes the dynamic Smagorinsky constant (Cs) for the LES model.
Definition les.c:42
PetscReal RuntimeWalltimeGuardUpdateEWMA(PetscBool has_previous, PetscReal previous_ewma_seconds, PetscReal latest_step_seconds, PetscReal alpha)
Update an EWMA estimate for timestep wall-clock duration.
Definition runloop.c:150
PetscReal RuntimeWalltimeGuardConservativeEstimate(PetscReal warmup_average_seconds, PetscReal ewma_seconds, PetscReal latest_step_seconds)
Return the conservative timestep estimate used by the walltime guard.
Definition runloop.c:162
PetscErrorCode UpdateSolverHistoryVectors(UserCtx *user)
Copies the current time step's solution fields into history vectors (e.g., U(t_n) -> U_o,...
Definition runloop.c:321
PetscReal RuntimeWalltimeGuardRequiredHeadroom(PetscReal min_seconds, PetscReal multiplier, PetscReal conservative_estimate_seconds)
Compute the required shutdown headroom from timestep estimate and floor.
Definition runloop.c:173
PetscBool RuntimeWalltimeGuardShouldTrigger(PetscInt completed_steps, PetscInt warmup_steps, PetscReal remaining_seconds, PetscReal min_seconds, PetscReal multiplier, PetscReal warmup_average_seconds, PetscReal ewma_seconds, PetscReal latest_step_seconds, PetscReal *required_headroom_seconds_out)
Decide whether the runtime walltime guard should stop before another step.
Definition runloop.c:184
PetscErrorCode GetOwnedCellRange(const DMDALocalInfo *info_nodes, PetscInt dim, PetscInt *xs_cell_global_out, PetscInt *xm_cell_local_out)
Determines the global starting index and number of CELLS owned by the current processor in a specifie...
Definition setup.c:2382
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.
Definition setup.c:2479
PetscErrorCode Cart2Contra(UserCtx *user)
Convert the ghosted Cartesian velocity field to contravariant face fluxes.
Definition setup.c:2881
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1755
PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
Parse a positive floating-point seconds value from runtime metadata.
Definition setup.c:18
static PetscErrorCode TestUpdateParticleFieldIEMRelaxation(void)
Tests IEM relaxation updates for particle-carried fields.
static PetscErrorCode TestMinimalFixtureMirrorsProductionDMLayout(void)
Tests that the shared minimal fixture mirrors the production DA contract.
static PetscErrorCode TestUpdateParticlePositionDiffusivityGradientOnly(void)
Tests particle position updates driven only by diffusivity-gradient drift.
static PetscErrorCode TestLocateAllParticlesInGridGuessPathResolvesLocalParticle(void)
Tests the guess-then-verify orchestrator path for a local particle with an unknown prior cell.
static PetscErrorCode TestWallModelVelocityHelpers(void)
Tests closed-form and iterative wall-model velocity helpers against inverse reconstructions.
static PetscErrorCode TestSetInitialInteriorFieldIgnoresNonUcontRequest(void)
Tests that non-Ucont requests do not modify interior field initialization.
static PetscErrorCode TestLocateParticleOrFindMigrationTargetCountsReSearch(void)
Verifies that later settlement passes increment re-search metrics.
static PetscErrorCode TestRuntimeWalltimeGuardParsesPositiveSeconds(void)
Tests parsing of positive runtime walltime metadata values.
static PetscErrorCode TestCart2ContraUsesFinalizedPeriodicUcat(void)
Verifies periodic Ucat endpoints are repaired before Cart2Contra terminal faces.
static PetscErrorCode TestSetInitialInteriorFieldCurvilinearConstantViaFlowDirection(void)
Tests curvilinear Constant IC: flow_direction selects the streamwise axis.
int main(int argc, char **argv)
Runs the unit-runtime PETSc test binary.
static PetscErrorCode TestInterpolateAllFieldsToSwarmConstantFields(void)
Tests direct interpolation from Eulerian fields to one localized swarm particle.
static PetscErrorCode SyncRuntimeFieldGhosts(UserCtx *user)
Synchronizes the minimal runtime fixture's global fields into their persistent local ghosts.
static PetscErrorCode TestDistributeParticlesRemainderHandling(void)
Tests particle distribution remainder handling across ranks.
static PetscErrorCode TestResetAllParticleStatusesLeavesLostParticlesUntouched(void)
Tests localized particle-status reset behavior for restart of the location workflow.
static PetscErrorCode SeedSingleParticle(UserCtx *user, PetscInt ci, PetscInt cj, PetscInt ck, PetscReal x, PetscReal y, PetscReal z, PetscReal wx, PetscReal wy, PetscReal wz, PetscInt status_value)
Seeds one localized swarm particle with the cell, position, weight, and status data used by runtime t...
static PetscErrorCode TestSetInitialInteriorFieldCartesianConstantSetsContravariantFlux(void)
Tests cartesian Constant IC: Cart2Contra sets contravariant flux via metric dot product.
static PetscErrorCode TestPopulateInitialUcontLoadsStagedUcat(void)
Tests loading a staged Ucat file IC and converting it to Ucont.
static PetscErrorCode TestPopulateInitialUcontLoadsStagedUcont(void)
Tests loading a staged Ucont file IC without Cartesian conversion.
static PetscErrorCode TestInterpolateAllFieldsToSwarmCornerAveragedConstantFields(void)
Tests the corner-averaged (legacy) interpolation path on constant fields.
static PetscErrorCode TestGetOwnedCellRangeSingleRankAccounting(void)
Tests owned-cell range accounting on a single MPI rank.
static PetscErrorCode TestUpdateParticlePositionWithoutBrownianContribution(void)
Tests particle position updates without Brownian forcing.
static PetscErrorCode TestComputeAndStoreNeighborRanksSingleRank(void)
Tests neighbor-rank discovery on a single MPI rank.
static PetscErrorCode TestCalculateParticleCountPerCellCountsGlobalCellIDs(void)
Tests particle counting by geometric cell IDs using the production +1 storage shift.
static PetscErrorCode TestValidateDrivenFlowConfigurationNoDrivenHandlers(void)
Tests driven-flow validation when no driven handlers are present.
static PetscErrorCode TestWallNoSlipAndFreeSlipHelpers(void)
Tests no-slip and free-slip wall helper kernels.
static PetscErrorCode TestMinimalFixtureRegistersProductionSwarmFields(void)
Tests that the shared swarm fixture registers the production field set.
static PetscErrorCode TestWallFunctionVectorWrappers(void)
Tests the vector wall-function wrappers on a tangential reference flow.
static PetscErrorCode TestSetInitialInteriorFieldZeroClearsInterior(void)
Tests zero IC clears physical-cell contravariant velocity.
static PetscErrorCode TestLocateAllParticlesInGridPriorCellFastPath(void)
Tests the location orchestrator fast path when a particle already carries a valid prior cell.
static PetscErrorCode TestComputeSmagorinskyConstantConstantModel(void)
Tests the constant Smagorinsky model helper path.
static PetscErrorCode TestCalculateBrownianDisplacementDeterministicSeed(void)
Tests Brownian displacement generation against a duplicated seeded RNG stream.
static PetscErrorCode TestScatterAllParticleFieldsToEulerFieldsAveragesPsi(void)
Tests particle-to-grid scattering using known cell occupancy and scalar values.
static PetscErrorCode TestUpdateSolverHistoryVectorsShiftsStates(void)
Tests solver history-vector shifting between time levels.
static PetscErrorCode TestUpdateParticleWeightsComputesExpectedRatios(void)
Tests particle weight updates against expected ratios.
static PetscErrorCode TestCheckAndRemoveOutOfBoundsParticlesRemovesEscapedParticle(void)
Tests direct removal of particles that leave every rank bounding box.
static PetscErrorCode TestSetInitialInteriorFieldPoiseuilleProfile(void)
Tests Poiseuille IC follows the discrete cross-stream profile and reaches zero at edges.
static PetscErrorCode TestCart2ContraConvertsCartesianField(void)
Tests spatially varying Cartesian-field conversion to contravariant fluxes.
static PetscErrorCode TestCheckAndRemoveLostParticlesRemovesLostEntries(void)
Tests direct removal of particles already marked LOST by the location workflow.
static PetscErrorCode TestRuntimeWalltimeGuardEstimatorHelpers(void)
Tests walltime-guard estimator helper calculations.
static PetscErrorCode TestIsParticleInsideBoundingBoxBasicCases(void)
Tests basic particle-inside-bounding-box classification cases.
static PetscErrorCode TestRuntimeWalltimeGuardTriggerDecision(void)
Tests runtime walltime-guard shutdown trigger decisions.
static PetscErrorCode TestUpdateAllParticlePositionsMovesSwarmEntries(void)
Tests swarm-wide particle position updates using the same transport path as the runtime loop.
static PetscErrorCode TestWallModelScalarHelpers(void)
Tests wall-model scalar helper kernels.
PetscErrorCode PicurvMakeTempDir(char *path, size_t path_len)
Creates a unique temporary directory for one test case.
PetscErrorCode PicurvCreateMinimalContexts(SimCtx **simCtx_out, UserCtx **user_out, PetscInt mx, PetscInt my, PetscInt mz)
Builds minimal SimCtx and UserCtx fixtures for C unit tests.
PetscErrorCode PicurvAssertRealNear(PetscReal expected, PetscReal actual, PetscReal tol, const char *context)
Asserts that two real values agree within tolerance.
PetscErrorCode PicurvDestroyMinimalContexts(SimCtx **simCtx_ptr, UserCtx **user_ptr)
Destroys minimal SimCtx/UserCtx fixtures and all owned PETSc objects.
PetscErrorCode PicurvCreateMinimalContextsWithPeriodicity(SimCtx **simCtx_out, UserCtx **user_out, PetscInt mx, PetscInt my, PetscInt mz, PetscBool x_periodic, PetscBool y_periodic, PetscBool z_periodic)
Builds minimal SimCtx and UserCtx fixtures for C unit tests with configurable periodicity.
PetscErrorCode PicurvCreateSwarmPair(UserCtx *user, PetscInt nlocal, const char *post_field_name)
Creates matched solver and post-processing swarms for tests.
PetscErrorCode PicurvDestroyRuntimeContext(SimCtx **simCtx_ptr)
Finalizes and frees a runtime context built by PicurvBuildTinyRuntimeContext.
PetscErrorCode PicurvRunTests(const char *suite_name, const PicurvTestCase *cases, size_t case_count)
Runs a named C test suite and prints pass/fail progress markers.
PetscErrorCode PicurvBuildTinyRuntimeContext(const char *bcs_contents, PetscBool enable_particles, SimCtx **simCtx_out, UserCtx **user_out, char *tmpdir, size_t tmpdir_len)
Builds a tiny runtime context through the real setup path for behavior-level tests.
PetscErrorCode PicurvAssertVecConstant(Vec vec, PetscScalar expected, PetscReal tol, const char *context)
Asserts that a PETSc vector is spatially constant within tolerance.
PetscErrorCode PicurvAssertIntEqual(PetscInt expected, PetscInt actual, const char *context)
Asserts that two integer values are equal.
PetscErrorCode PicurvAssertBool(PetscBool value, const char *context)
Asserts that one boolean condition is true.
PetscErrorCode PicurvRemoveTempDir(const char *path)
Recursively removes a temporary directory created by PicurvMakeTempDir.
Shared declarations for the PICurv C test fixture and assertion layer.
Named test case descriptor consumed by PicurvRunTests.
PetscMPIInt rank_zm
Definition variables.h:197
@ CONSTANT_SMAGORINSKY
Definition variables.h:520
PetscReal icVelocityPhysical
Definition variables.h:747
Vec lDiffusivityGradient
Definition variables.h:908
@ PERIODIC
Definition variables.h:290
Cmpnts vel
Definition variables.h:184
PetscInt ys_cell
Definition variables.h:202
PetscInt xs_cell
Definition variables.h:202
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:896
PetscMPIInt rank_yp
Definition variables.h:196
PetscInt64 searchLocatedCount
Definition variables.h:239
PetscInt64 searchLostCount
Definition variables.h:240
PetscInt cell[3]
Definition variables.h:182
InitialConditionMode initialConditionMode
Definition variables.h:742
ParticleLocationStatus
Defines the state of a particle with respect to its location and migration status during the iterativ...
Definition variables.h:135
@ LOST
Definition variables.h:139
@ NEEDS_LOCATION
Definition variables.h:136
@ ACTIVE_AND_LOCATED
Definition variables.h:137
PetscMPIInt rank_ym
Definition variables.h:196
FlowDirection flowDirection
Definition variables.h:746
PetscMPIInt rank_xp
Definition variables.h:195
PetscInt KM
Definition variables.h:885
PetscInt64 traversalStepsSum
Definition variables.h:241
PetscReal ren
Definition variables.h:732
Vec lUcont_rm1
Definition variables.h:912
PetscInt zm_cell
Definition variables.h:203
Cmpnts max_coords
Maximum x, y, z coordinates of the bounding box.
Definition variables.h:171
PetscInt zs_cell
Definition variables.h:202
Cmpnts diffusivitygradient
Definition variables.h:189
PetscInt64 searchPopulation
Definition variables.h:238
PetscReal dt
Definition variables.h:699
RankNeighbors neighbors
Definition variables.h:888
Vec lPsi
Definition variables.h:953
PetscInt currentSettlementPass
Definition variables.h:250
Vec DiffusivityGradient
Definition variables.h:908
Vec lCs
Definition variables.h:935
Vec Ucont
Definition variables.h:904
PetscInt StartStep
Definition variables.h:694
Cmpnts min_coords
Minimum x, y, z coordinates of the bounding box.
Definition variables.h:170
PetscScalar x
Definition variables.h:101
Cmpnts loc
Definition variables.h:183
PetscInt64 reSearchCount
Definition variables.h:242
char * current_io_directory
Definition variables.h:711
PetscInt xm_cell
Definition variables.h:203
Vec lUcont_o
Definition variables.h:911
PetscInt64 bboxGuessFallbackCount
Definition variables.h:248
InterpolationMethod interpolationMethod
Definition variables.h:801
RankCellInfo * RankCellInfoMap
Definition variables.h:951
PetscInt ym_cell
Definition variables.h:203
Vec Ucat_o
Definition variables.h:911
PetscInt64 bboxGuessSuccessCount
Definition variables.h:247
BoundingBox * bboxlist
Definition variables.h:799
PetscMPIInt rank_xm
Definition variables.h:195
PetscInt64 maxParticlePassDepth
Definition variables.h:249
PetscScalar z
Definition variables.h:101
@ INTERP_CORNER_AVERAGED
Definition variables.h:564
Vec Ucat
Definition variables.h:904
Vec ParticleCount
Definition variables.h:952
PetscReal Const_CS
Definition variables.h:791
Vec Ucont_o
Definition variables.h:911
PetscInt JM
Definition variables.h:885
@ FLOW_DIR_POS_ZETA
Definition variables.h:275
char initialConditionDirectory[PETSC_MAX_PATH_LEN]
Definition variables.h:744
@ IC_MODE_CONSTANT_CARTESIAN
Definition variables.h:151
@ IC_MODE_POISEUILLE
Definition variables.h:152
@ IC_MODE_CONSTANT_STREAMWISE
Definition variables.h:153
@ IC_MODE_FILE
Definition variables.h:154
@ IC_MODE_ZERO
Definition variables.h:150
Cmpnts InitialConstantContra
Definition variables.h:745
PetscMPIInt rank_zp
Definition variables.h:197
Vec Ucont_rm1
Definition variables.h:912
SearchMetricsState searchMetrics
Definition variables.h:809
PetscReal diffusivity
Definition variables.h:188
Vec lUcont
Definition variables.h:904
PetscInt step
Definition variables.h:692
Vec Diffusivity
Definition variables.h:907
PetscRandom BrownianMotionRNG
Definition variables.h:810
PetscInt GridOrientation
Definition variables.h:889
DMDALocalInfo info
Definition variables.h:883
Vec lUcat
Definition variables.h:904
PetscScalar y
Definition variables.h:101
@ IC_FIELD_UCONT
Definition variables.h:160
@ IC_FIELD_UCAT
Definition variables.h:159
PetscMPIInt size
Definition variables.h:688
PetscInt IM
Definition variables.h:885
char _io_context_buffer[PETSC_MAX_PATH_LEN]
Definition variables.h:710
Cmpnts weights
Definition variables.h:185
@ NUM_FACES
Definition variables.h:145
PetscInt les
Definition variables.h:789
Vec lDiffusivity
Definition variables.h:907
BCType mathematical_type
Definition variables.h:366
PetscInt64 searchAttempts
Definition variables.h:237
InitialConditionField initialConditionField
Definition variables.h:743
PetscInt64 PID
Definition variables.h:181
PetscInt64 maxTraversalFailCount
Definition variables.h:244
Vec Psi
Definition variables.h:953
Vec P_o
Definition variables.h:911
@ BC_FACE_NEG_X
Definition variables.h:260
@ BC_FACE_POS_X
Definition variables.h:260
Defines a 3D axis-aligned bounding box.
Definition variables.h:169
A 3D point or vector with PetscScalar components.
Definition variables.h:100
Defines a particle's core properties for Lagrangian tracking.
Definition variables.h:180
The master context for the entire simulation.
Definition variables.h:684
User-defined context containing data specific to a single computational grid level.
Definition variables.h:876
Header file for particle location functions using the walking search algorithm.
PetscErrorCode LocateParticleOrFindMigrationTarget(UserCtx *user, Particle *particle, ParticleLocationStatus *status_out)
Locates a particle's host cell or identifies its migration target using a robust walk search.
void wall_function_loglaw(UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
Applies log-law wall function with roughness correction.
double find_utau_loglaw(double velocity, double wall_distance, double roughness_length)
Solves for friction velocity using simple log-law (explicit formula)
double u_Werner(double kinematic_viscosity, double wall_distance, double friction_velocity)
Computes velocity using Werner-Wengle wall function.
void wall_function(UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
Applies standard wall function with Werner-Wengle model.
double find_utau_hydset(double kinematic_viscosity, double known_velocity, double wall_distance, double initial_guess, double roughness_height)
Solves for friction velocity using Newton-Raphson iteration.
double u_Cabot(double kinematic_viscosity, double wall_distance, double friction_velocity, double pressure_gradient_tangent, double wall_shear_stress)
Computes velocity using Cabot wall function.
double u_loglaw(double wall_distance, double friction_velocity, double roughness_length)
Computes velocity using simple log-law (smooth wall with roughness offset)
void wall_function_Cabot(UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z, double pressure_gradient_x, double pressure_gradient_y, double pressure_gradient_z, int iteration_count)
Applies Cabot non-equilibrium wall function with pressure gradients.
void noslip(UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
Applies no-slip wall boundary condition with linear interpolation.
void find_utau_Cabot(double kinematic_viscosity, double velocity, double wall_distance, double initial_guess, double pressure_gradient_tangent, double pressure_gradient_normal, double *friction_velocity, double *wall_shear_velocity, double *wall_shear_normal)
Solves for friction velocity using Cabot wall function.
double integrate_1(double kinematic_viscosity, double wall_distance, double friction_velocity, int integration_mode)
Integrates eddy viscosity profile from wall to distance y.
double f_hydset(double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
Residual function for friction velocity equation (log-law with roughness)
double E_coeff(double friction_velocity, double roughness_height, double kinematic_viscosity)
Computes roughness-modified log-law coefficient E.
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
double find_utau_Werner(double kinematic_viscosity, double velocity, double wall_distance, double initial_guess)
Solves for friction velocity using Werner-Wengle wall function.
void freeslip(UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
Applies free-slip wall boundary condition.