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));
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));
38 PetscFunctionReturn(0);
54 PetscInt status_value)
56 PetscReal *positions = NULL;
57 PetscReal *weights = NULL;
58 PetscInt *cell_ids = NULL;
59 PetscInt *status = NULL;
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));
76 status[0] = status_value;
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);
90 PetscInt particles_per_rank = -1;
91 PetscInt remainder = -1;
93 PetscFunctionBeginUser;
95 PetscCall(
PicurvAssertIntEqual(4, particles_per_rank,
"rank 0 should receive one remainder particle"));
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);
112 PetscFunctionBeginUser;
113 PetscCall(PetscMemzero(&bbox,
sizeof(bbox)));
114 PetscCall(PetscMemzero(&particle,
sizeof(particle)));
123 particle.
loc.
x = 0.25;
124 particle.
loc.
y = 1.0;
125 particle.
loc.
z = 2.5;
128 particle.
loc.
x = 1.5;
130 PetscFunctionReturn(0);
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};
142 PetscFunctionBeginUser;
143 PetscCall(PetscMemzero(&particle,
sizeof(particle)));
153 PetscFunctionReturn(0);
165 PetscFunctionBeginUser;
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;
187 PetscFunctionReturn(0);
199 PetscFunctionBeginUser;
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;
221 PetscFunctionReturn(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;
237 PetscFunctionBeginUser;
240 mean_val + (1.0 - mean_val) * PetscExpReal(-(c_model * diffusivity / PetscPowReal(cell_vol, 0.6666667)) * dt),
243 "IEM update should match analytical relaxation"));
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);
258 PetscFunctionBeginUser;
260 PetscCall(VecSet(user->
Ucont, 7.0));
266 PetscFunctionReturn(0);
280 PetscFunctionBeginUser;
286 PetscCall(VecSet(user->
Ucont, 0.0));
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));
296 PetscFunctionReturn(0);
312 PetscFunctionBeginUser;
322 PetscCall(VecSet(user->
Ucont, 0.0));
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));
332 PetscFunctionReturn(0);
344 PetscFunctionBeginUser;
347 PetscCall(VecSet(user->
Ucont, 7.0));
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));
357 PetscFunctionReturn(0);
369 PetscFunctionBeginUser;
375 PetscCall(VecZeroEntries(user->
Ucont));
378 PetscCall(DMDAVecGetArrayRead(user->
fda, user->
Ucont, &ucont));
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));
387 PetscFunctionReturn(0);
400 PetscFunctionBeginUser;
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));
411 PetscCall(DMDAVecGetArrayRead(user->
fda, user->
Ucont, &ucont));
413 "Xi flux uses face-interpolated Ucat"));
415 "Eta flux uses Cartesian y velocity"));
417 "Zeta flux uses Cartesian z velocity"));
418 PetscCall(DMDAVecRestoreArrayRead(user->
fda, user->
Ucont, &ucont));
421 PetscFunctionReturn(0);
428 const char *cell_fields[]={
"Ucat"}, *staggered_fields[]={
"Ucont"};
429 PetscFunctionBeginUser;
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));
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));
449 PetscFunctionReturn(0);
460 char tmpdir[PETSC_MAX_PATH_LEN];
462 PetscFunctionBeginUser;
465 PetscCall(VecSet(user->
Ucat, 2.5));
471 PetscCall(VecZeroEntries(user->
Ucat));
472 PetscCall(VecZeroEntries(user->
Ucont));
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));
486 PetscFunctionReturn(0);
496 char tmpdir[PETSC_MAX_PATH_LEN];
498 PetscFunctionBeginUser;
501 PetscCall(VecSet(user->
Ucont, 4.5));
507 PetscCall(VecZeroEntries(user->
Ucont));
516 PetscFunctionReturn(0);
527 PetscReal *velocity = NULL;
528 PetscReal *diffusivity = NULL;
529 PetscReal *diffusivity_gradient = NULL;
531 PetscFunctionBeginUser;
534 PetscCall(VecSet(user->
Ucat, 2.0));
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;
549 PetscCall(
SeedSingleParticle(user, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
ACTIVE_AND_LOCATED));
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));
568 PetscFunctionReturn(0);
579 PetscReal *velocity = NULL;
580 PetscReal *diffusivity = NULL;
581 PetscReal *diffusivity_gradient = NULL;
583 PetscFunctionBeginUser;
587 PetscCall(VecSet(user->
Ucat, 2.0));
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;
602 PetscCall(
SeedSingleParticle(user, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
ACTIVE_AND_LOCATED));
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));
621 PetscFunctionReturn(0);
631 PetscInt *cell_ids = NULL;
632 PetscReal *psi = NULL;
633 PetscReal ***psi_grid = NULL;
635 PetscFunctionBeginUser;
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));
644 PetscCall(DMSwarmGetField(user->
swarm,
"Psi", NULL, NULL, (
void **)&psi));
647 PetscCall(DMSwarmRestoreField(user->
swarm,
"Psi", NULL, NULL, (
void **)&psi));
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));
656 PetscFunctionReturn(0);
666 PetscInt *cell_ids = NULL;
667 PetscReal ***counts = NULL;
669 PetscFunctionBeginUser;
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));
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));
687 PetscFunctionReturn(0);
697 PetscInt *status = NULL;
699 PetscFunctionBeginUser;
703 PetscCall(DMSwarmGetField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
707 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
711 PetscCall(DMSwarmGetField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
715 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
718 PetscFunctionReturn(0);
728 PetscReal *positions = NULL;
729 PetscInt removed_local = 0;
730 PetscInt removed_global = 0;
733 PetscFunctionBeginUser;
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));
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"));
749 PetscFunctionReturn(0);
759 PetscInt *status = NULL;
760 PetscInt removed_local = 0;
761 PetscInt removed_global = 0;
764 PetscFunctionBeginUser;
768 PetscCall(DMSwarmGetField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
772 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
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"));
781 PetscFunctionReturn(0);
791 char tmpdir[PETSC_MAX_PATH_LEN];
795 PetscFunctionBeginUser;
799 "runtime setup path should initialize the Brownian RNG"));
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"));
815 PetscFunctionReturn(0);
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;
834 PetscFunctionBeginUser;
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));
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;
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;
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));
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));
880 PetscFunctionReturn(0);
889 PetscReal *positions = NULL;
890 PetscReal *weights = NULL;
891 PetscInt *cell_ids = NULL;
892 PetscInt *status = NULL;
893 PetscInt64 *pid = NULL;
895 PetscFunctionBeginUser;
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;
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));
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"));
934 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
935 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_CellID", NULL, NULL, (
void **)&cell_ids));
938 PetscFunctionReturn(0);
947 PetscReal *positions = NULL;
948 PetscReal *weights = NULL;
949 PetscInt *cell_ids = NULL;
950 PetscInt *status = NULL;
951 PetscInt64 *pid = NULL;
953 PetscFunctionBeginUser;
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;
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));
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"));
992 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_location_status", NULL, NULL, (
void **)&status));
993 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_CellID", NULL, NULL, (
void **)&cell_ids));
996 PetscFunctionReturn(0);
1008 PetscFunctionBeginUser;
1009 PetscCall(PetscMemzero(&particle,
sizeof(particle)));
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;
1036 PetscFunctionReturn(0);
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};
1049 PetscFunctionBeginUser;
1050 noslip(NULL, 2.0, 1.0, wall_velocity, reference_velocity, &boundary_velocity, 1.0, 0.0, 0.0);
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);
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;
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"));
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"));
1083 PetscFunctionReturn(0);
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;
1101 PetscFunctionBeginUser;
1102 utau_loglaw =
find_utau_loglaw(target_velocity, wall_distance, roughness_length);
1104 "simple log-law inversion should reconstruct the target velocity"));
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"));
1109 "Werner-Wengle inversion should reconstruct the target velocity"));
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"));
1117 "zero normal pressure gradient should keep Cabot normal wall shear at zero"));
1118 PetscFunctionReturn(0);
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;
1132 PetscFunctionBeginUser;
1134 simCtx->
ren = 1000.0;
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"));
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"));
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"));
1153 PetscFunctionReturn(0);
1164 PetscFunctionBeginUser;
1168 PetscFunctionReturn(0);
1178 PetscReal ***lcs = NULL;
1180 PetscFunctionBeginUser;
1182 PetscCall(DMCreateGlobalVector(user->
da, &user->
CS));
1183 PetscCall(DMCreateLocalVector(user->
da, &user->
lCs));
1191 PetscCall(DMDAVecGetArrayRead(user->
da, user->
lCs, &lcs));
1193 "constant Smagorinsky branch should refresh local CS"));
1194 PetscCall(DMDAVecRestoreArrayRead(user->
da, user->
lCs, &lcs));
1197 PetscFunctionReturn(0);
1208 PetscInt mx = 0, my = 0, mz = 0;
1210 PetscFunctionBeginUser;
1213 PetscCall(DMDAGetInfo(user->
da, NULL, &mx, &my, &mz, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1221 PetscCall(DMGetCoordinateDM(user->
da, &coord_dm));
1223 "minimal fixture should derive fda from the coordinate-DM path"));
1226 PetscFunctionReturn(0);
1237 void *field_ptr = NULL;
1239 PetscFunctionBeginUser;
1243 PetscCall(DMSwarmGetField(user->
swarm,
"position", &bs, NULL, &field_ptr));
1245 PetscCall(
PicurvAssertBool((PetscBool)(field_ptr != NULL),
"position field should be retrievable"));
1246 PetscCall(DMSwarmRestoreField(user->
swarm,
"position", &bs, NULL, &field_ptr));
1248 PetscCall(DMSwarmGetField(user->
swarm,
"velocity", &bs, NULL, &field_ptr));
1250 PetscCall(DMSwarmRestoreField(user->
swarm,
"velocity", &bs, NULL, &field_ptr));
1252 PetscCall(DMSwarmGetField(user->
swarm,
"DMSwarm_CellID", &bs, NULL, &field_ptr));
1254 PetscCall(DMSwarmRestoreField(user->
swarm,
"DMSwarm_CellID", &bs, NULL, &field_ptr));
1256 PetscCall(DMSwarmGetField(user->
swarm,
"weight", &bs, NULL, &field_ptr));
1258 PetscCall(DMSwarmRestoreField(user->
swarm,
"weight", &bs, NULL, &field_ptr));
1260 PetscCall(DMSwarmGetField(user->
swarm,
"Diffusivity", &bs, NULL, &field_ptr));
1262 PetscCall(DMSwarmRestoreField(user->
swarm,
"Diffusivity", &bs, NULL, &field_ptr));
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));
1268 PetscCall(DMSwarmGetField(user->
swarm,
"Psi", &bs, NULL, &field_ptr));
1270 PetscCall(DMSwarmRestoreField(user->
swarm,
"Psi", &bs, NULL, &field_ptr));
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));
1277 PetscFunctionReturn(0);
1288 PetscFunctionBeginUser;
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));
1309 PetscFunctionReturn(0);
1319 PetscInt xs = -1, ys = -1, zs = -1;
1320 PetscInt xm = -1, ym = -1, zm = -1;
1322 PetscFunctionBeginUser;
1337 PetscFunctionReturn(0);
1348 PetscFunctionBeginUser;
1360 PetscFunctionReturn(0);
1368 PetscReal seconds = 0.0;
1370 PetscFunctionBeginUser;
1377 PetscFunctionReturn(0);
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;
1391 PetscFunctionBeginUser;
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);
1411 PetscBool should_trigger = PETSC_FALSE;
1412 PetscReal required_headroom = 0.0;
1414 PetscFunctionBeginUser;
1416 PetscCall(
PicurvAssertBool((PetscBool)!should_trigger,
"guard should not trigger before warmup completes"));
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"));
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);
1433 PetscErrorCode ierr;
1478 ierr = PetscInitialize(&argc, &argv, NULL,
"PICurv runtime-kernel tests");
1483 ierr =
PicurvRunTests(
"unit-runtime", cases,
sizeof(cases) /
sizeof(cases[0]));
1489 ierr = PetscFinalize();
PetscErrorCode Validate_DrivenFlowConfiguration(UserCtx *user)
(Private) Validates all consistency rules for a driven flow (channel/pipe) setup.
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.
PetscErrorCode ComputeSmagorinskyConstant(UserCtx *user)
Computes the dynamic Smagorinsky constant (Cs) for the LES model.
PetscReal RuntimeWalltimeGuardUpdateEWMA(PetscBool has_previous, PetscReal previous_ewma_seconds, PetscReal latest_step_seconds, PetscReal alpha)
Update an EWMA estimate for timestep wall-clock duration.
PetscReal RuntimeWalltimeGuardConservativeEstimate(PetscReal warmup_average_seconds, PetscReal ewma_seconds, PetscReal latest_step_seconds)
Return the conservative timestep estimate used by the walltime guard.
PetscErrorCode UpdateSolverHistoryVectors(UserCtx *user)
Copies the current time step's solution fields into history vectors (e.g., U(t_n) -> U_o,...
PetscReal RuntimeWalltimeGuardRequiredHeadroom(PetscReal min_seconds, PetscReal multiplier, PetscReal conservative_estimate_seconds)
Compute the required shutdown headroom from timestep estimate and floor.
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.
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...
PetscErrorCode ComputeAndStoreNeighborRanks(UserCtx *user)
Computes and stores the Cartesian neighbor ranks for the DMDA decomposition.
PetscErrorCode Cart2Contra(UserCtx *user)
Convert the ghosted Cartesian velocity field to contravariant face fluxes.
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
PetscBool RuntimeWalltimeGuardParsePositiveSeconds(const char *text, PetscReal *seconds_out)
Parse a positive floating-point seconds value from runtime metadata.
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.
PetscReal icVelocityPhysical
BoundaryFaceConfig boundary_faces[6]
PetscInt64 searchLocatedCount
PetscInt64 searchLostCount
InitialConditionMode initialConditionMode
ParticleLocationStatus
Defines the state of a particle with respect to its location and migration status during the iterativ...
FlowDirection flowDirection
PetscInt64 traversalStepsSum
Cmpnts max_coords
Maximum x, y, z coordinates of the bounding box.
Cmpnts diffusivitygradient
PetscInt64 searchPopulation
PetscInt currentSettlementPass
Cmpnts min_coords
Minimum x, y, z coordinates of the bounding box.
char * current_io_directory
PetscInt64 bboxGuessFallbackCount
InterpolationMethod interpolationMethod
RankCellInfo * RankCellInfoMap
PetscInt64 bboxGuessSuccessCount
PetscInt64 maxParticlePassDepth
char initialConditionDirectory[PETSC_MAX_PATH_LEN]
@ IC_MODE_CONSTANT_CARTESIAN
@ IC_MODE_CONSTANT_STREAMWISE
Cmpnts InitialConstantContra
SearchMetricsState searchMetrics
PetscRandom BrownianMotionRNG
char _io_context_buffer[PETSC_MAX_PATH_LEN]
PetscInt64 searchAttempts
InitialConditionField initialConditionField
PetscInt64 maxTraversalFailCount
Defines a 3D axis-aligned bounding box.
A 3D point or vector with PetscScalar components.
Defines a particle's core properties for Lagrangian tracking.
The master context for the entire simulation.
User-defined context containing data specific to a single computational grid level.
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.