PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
test_solver_kernels.c
Go to the documentation of this file.
1/**
2 * @file test_solver_kernels.c
3 * @brief C unit tests for solver-side analytical and LES helper kernels.
4 */
5
6#include "test_support.h"
7
8#include "AnalyticalSolutions.h"
9#include "BodyForces.h"
10#include "Filter.h"
11#include "les.h"
12#include "solvers.h"
13/**
14 * @brief Tests LES test-filter helper paths for representative cases.
15 */
16
17static PetscErrorCode TestLESTestFilterPaths(void)
18{
19 SimCtx simCtx;
20 double values[3][3][3];
21 double weights[3][3][3];
22
23 PetscFunctionBeginUser;
24 PetscCall(PetscMemzero(&simCtx, sizeof(simCtx)));
25 for (PetscInt k = 0; k < 3; ++k) {
26 for (PetscInt j = 0; j < 3; ++j) {
27 for (PetscInt i = 0; i < 3; ++i) {
28 values[k][j][i] = 2.0;
29 weights[k][j][i] = 1.0;
30 }
31 }
32 }
33
34 simCtx.testfilter_ik = 1;
35 PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
36 "Simpson-rule filter should preserve a constant field"));
37
38 simCtx.testfilter_ik = 0;
39 PetscCall(PicurvAssertRealNear(2.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
40 "box filter should preserve a constant field"));
41
42 for (PetscInt k = 0; k < 3; ++k) {
43 for (PetscInt j = 0; j < 3; ++j) {
44 for (PetscInt i = 0; i < 3; ++i) {
45 weights[k][j][i] = 0.0;
46 }
47 }
48 }
49 PetscCall(PicurvAssertRealNear(0.0, ApplyLESTestFilter(&simCtx, values, weights), 1.0e-12,
50 "box filter should return zero when all weights are zero"));
51 PetscFunctionReturn(0);
52}
53/**
54 * @brief Tests analytical geometry selection for supported analytical solutions.
55 */
56
57static PetscErrorCode TestAnalyticalGeometrySelection(void)
58{
59 SimCtx *simCtx = NULL;
60 UserCtx *user = NULL;
61 PetscErrorCode ierr_grid = 0;
62 PetscErrorCode ierr_non_square = 0;
63
64 PetscFunctionBeginUser;
65 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 8, 8, 8));
66 PetscCall(PicurvAssertBool(AnalyticalTypeRequiresCustomGeometry("TGV3D"),
67 "TGV3D should require custom geometry"));
68 PetscCall(PicurvAssertBool((PetscBool)!AnalyticalTypeRequiresCustomGeometry("ZERO_FLOW"),
69 "ZERO_FLOW should not require custom geometry"));
70
71 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
72 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
73 ierr_grid = SetAnalyticalGridInfo(user);
74 PetscCall(PetscPopErrorHandler());
75 PetscCall(PicurvAssertBool((PetscBool)(ierr_grid != 0),
76 "SetAnalyticalGridInfo should reject analytical types without custom geometry"));
77
78 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "TGV3D", sizeof(simCtx->AnalyticalSolutionType)));
79 simCtx->block_number = 2;
80 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
81 ierr_non_square = SetAnalyticalGridInfo(user);
82 PetscCall(PetscPopErrorHandler());
83 PetscCall(PicurvAssertBool((PetscBool)(ierr_non_square != 0),
84 "TGV3D multi-block setup should reject non-square block counts"));
85
86 simCtx->block_number = 1;
87 PetscCall(SetAnalyticalGridInfo(user));
88 PetscCall(PicurvAssertRealNear(0.0, user->Min_X, 1.0e-12, "TGV3D single-block xmin"));
89 PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_X, 1.0e-12, "TGV3D single-block xmax"));
90
91 simCtx->block_number = 4;
92 user->_this = 3;
93 PetscCall(SetAnalyticalGridInfo(user));
94 PetscCall(PicurvAssertRealNear(PETSC_PI, user->Min_X, 1.0e-12, "TGV3D multi-block xmin should reflect the block column"));
95 PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_X, 1.0e-12, "TGV3D multi-block xmax should reflect the block column"));
96 PetscCall(PicurvAssertRealNear(PETSC_PI, user->Min_Y, 1.0e-12, "TGV3D multi-block ymin should reflect the block row"));
97 PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_Y, 1.0e-12, "TGV3D multi-block ymax should reflect the block row"));
98 PetscCall(PicurvAssertRealNear(2.0 * PETSC_PI, user->Max_Z, 1.0e-12, "TGV3D multi-block zmax should span the full domain"));
99
100 simCtx->block_number = 1;
101 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
102 PetscFunctionReturn(0);
103}
104/**
105 * @brief Tests analytical scalar verification helper routines.
106 */
107
108static PetscErrorCode TestAnalyticalScalarVerificationHelpers(void)
109{
110 SimCtx *simCtx = NULL;
111 UserCtx *user = NULL;
112 Vec target = NULL;
113 PetscReal value = 0.0;
114 PetscReal ***target_arr = NULL;
115 PetscReal *positions = NULL;
116 PetscReal *psi = NULL;
117
118 PetscFunctionBeginUser;
119 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
120 PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
121 simCtx->verificationScalar.enabled = PETSC_TRUE;
122 PetscCall(PetscStrncpy(simCtx->verificationScalar.mode,
123 "analytical",
124 sizeof(simCtx->verificationScalar.mode)));
125
126 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
127 "CONSTANT",
128 sizeof(simCtx->verificationScalar.profile)));
129 simCtx->verificationScalar.value = 2.5;
130 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.2, 0.3, 0.4, 0.0, &value));
131 PetscCall(PicurvAssertRealNear(2.5, value, 1.0e-12,
132 "constant scalar profile should evaluate to the configured value"));
133
134 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
135 positions[0] = 0.1; positions[1] = 0.2; positions[2] = 0.3;
136 positions[3] = 0.8; positions[4] = 0.6; positions[5] = 0.4;
137 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
138
139 PetscCall(SetAnalyticalScalarFieldOnParticles(user, "Psi"));
140 PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
141 PetscCall(PicurvAssertRealNear(2.5, psi[0], 1.0e-12,
142 "SetAnalyticalScalarFieldOnParticles should overwrite the first particle scalar"));
143 PetscCall(PicurvAssertRealNear(2.5, psi[1], 1.0e-12,
144 "SetAnalyticalScalarFieldOnParticles should overwrite the second particle scalar"));
145 PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
146
147 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
148 "LINEAR_X",
149 sizeof(simCtx->verificationScalar.profile)));
150 simCtx->verificationScalar.phi0 = 1.0;
151 simCtx->verificationScalar.slope_x = 2.0;
152 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.25, 0.0, 0.0, 0.0, &value));
153 PetscCall(PicurvAssertRealNear(1.5, value, 1.0e-12,
154 "linear-x scalar profile should evaluate phi0 + slope_x * x"));
155
156 PetscCall(VecDuplicate(user->Psi, &target));
157 PetscCall(SetAnalyticalScalarFieldAtCellCenters(user, target));
158 PetscCall(DMDAVecGetArrayRead(user->da, target, &target_arr));
159 PetscCall(PicurvAssertRealNear(1.25, target_arr[1][1][1], 1.0e-12,
160 "cell-center scalar fill should use the physical x center coordinate"));
161 PetscCall(DMDAVecRestoreArrayRead(user->da, target, &target_arr));
162 PetscCall(VecDestroy(&target));
163
164 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
165 "SIN_PRODUCT",
166 sizeof(simCtx->verificationScalar.profile)));
167 simCtx->verificationScalar.amplitude = 3.0;
168 simCtx->verificationScalar.kx = PETSC_PI;
169 simCtx->verificationScalar.ky = PETSC_PI;
170 simCtx->verificationScalar.kz = PETSC_PI;
171 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.5, 0.5, 0.5, 0.0, &value));
172 PetscCall(PicurvAssertRealNear(3.0, value, 1.0e-12,
173 "sin-product scalar profile should peak at pi/2 in each coordinate"));
174
175 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
176 PetscFunctionReturn(0);
177}
178/**
179 * @brief Tests analytical solution engine ZERO_FLOW, UNIFORM_FLOW, and unknown-type dispatch.
180 */
181
182static PetscErrorCode TestAnalyticalSolutionEngineDispatch(void)
183{
184 SimCtx *simCtx = NULL;
185 UserCtx *user = NULL;
186 PetscErrorCode ierr_unknown = 0;
187 Cmpnts ***ucat = NULL;
188 Cmpnts ***ubcs = NULL;
189
190 PetscFunctionBeginUser;
191 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
192 PetscCall(VecSet(user->Ucat, 3.0));
193 PetscCall(VecSet(user->P, 5.0));
194 PetscCall(VecSet(user->Bcs.Ubcs, 7.0));
195
196 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
197 PetscCall(AnalyticalSolutionEngine(simCtx));
198 PetscCall(PicurvAssertVecConstant(user->Ucat, 0.0, 1.0e-12, "ZERO_FLOW should zero the Eulerian velocity field"));
199 PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "ZERO_FLOW should zero the pressure field"));
200 PetscCall(PicurvAssertVecConstant(user->Bcs.Ubcs, 0.0, 1.0e-12, "ZERO_FLOW should zero boundary-condition velocity data"));
201
202 /* UNIFORM_FLOW now works in curvilinear form (sets Ucont via metric dot products,
203 derives Ucat via Contra2Cart). The test grid needs identity metrics so the
204 transformation is invertible and Ucat recovers the physical velocity exactly. */
205 {
206 Cmpnts ***l_csi, ***l_eta, ***l_zet;
207 DMDALocalInfo linfo;
208 PetscCall(DMDAGetLocalInfo(user->fda, &linfo));
209 PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
210 PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
211 PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
212 for (PetscInt k = linfo.zs; k < linfo.zs + linfo.zm; k++)
213 for (PetscInt j = linfo.ys; j < linfo.ys + linfo.ym; j++)
214 for (PetscInt i = linfo.xs; i < linfo.xs + linfo.xm; i++) {
215 l_csi[k][j][i] = (Cmpnts){1.0, 0.0, 0.0};
216 l_eta[k][j][i] = (Cmpnts){0.0, 1.0, 0.0};
217 l_zet[k][j][i] = (Cmpnts){0.0, 0.0, 1.0};
218 }
219 PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
220 PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
221 PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
222 }
223
224 simCtx->AnalyticalUniformVelocity.x = 1.25;
225 simCtx->AnalyticalUniformVelocity.y = -0.5;
226 simCtx->AnalyticalUniformVelocity.z = 0.75;
227 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
228 PetscCall(AnalyticalSolutionEngine(simCtx));
229 PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "UNIFORM_FLOW should keep the pressure field zero"));
230
231 Cmpnts ***ucont = NULL;
232 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
233 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
234 PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
235 PetscCall(PicurvAssertRealNear(1.25, ucat[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should impose the configured x velocity"));
236 PetscCall(PicurvAssertRealNear(-0.5, ucat[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should impose the configured y velocity"));
237 PetscCall(PicurvAssertRealNear(0.75, ucat[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should impose the configured z velocity"));
238 PetscCall(PicurvAssertRealNear(1.25, ucont[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should set contravariant x flux (identity metric)"));
239 PetscCall(PicurvAssertRealNear(-0.5, ucont[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should set contravariant y flux (identity metric)"));
240 PetscCall(PicurvAssertRealNear(0.75, ucont[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should set contravariant z flux (identity metric)"));
241 PetscCall(PicurvAssertRealNear(1.25, ubcs[0][1][1].x, 1.0e-12, "UNIFORM_FLOW should populate boundary x velocity"));
242 PetscCall(PicurvAssertRealNear(-0.5, ubcs[0][1][1].y, 1.0e-12, "UNIFORM_FLOW should populate boundary y velocity"));
243 PetscCall(PicurvAssertRealNear(0.75, ubcs[0][1][1].z, 1.0e-12, "UNIFORM_FLOW should populate boundary z velocity"));
244 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
245 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
246 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
247
248 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "NOT_A_REAL_ANALYTICAL_TYPE", sizeof(simCtx->AnalyticalSolutionType)));
249 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
250 ierr_unknown = AnalyticalSolutionEngine(simCtx);
251 PetscCall(PetscPopErrorHandler());
252 PetscCall(PicurvAssertBool((PetscBool)(ierr_unknown != 0),
253 "AnalyticalSolutionEngine should reject unknown analytical type strings"));
254
255 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
256 PetscFunctionReturn(0);
257}
258/**
259 * @brief Tests exact Taylor-Green samples on selected Eulerian interior and boundary points.
260 */
261static PetscErrorCode TestAnalyticalSolutionEngineTaylorGreenSamples(void)
262{
263 SimCtx *simCtx = NULL;
264 UserCtx *user = NULL;
265 Cmpnts ***cent = NULL;
266 Cmpnts ***cent_x = NULL;
267 Cmpnts ***cent_y = NULL;
268 Cmpnts ***cent_z = NULL;
269 Cmpnts ***ucat = NULL;
270 Cmpnts ***ubcs = NULL;
271 PetscReal ***p = NULL;
272 const PetscReal vel_decay = PetscExpReal(-0.5);
273
274 PetscFunctionBeginUser;
275 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
276 simCtx->ren = 2.0;
277 simCtx->ti = 0.5;
278 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "TGV3D", sizeof(simCtx->AnalyticalSolutionType)));
279
280 PetscCall(DMDAVecGetArray(user->fda, user->Cent, &cent));
281 PetscCall(DMDAVecGetArray(user->fda, user->Centx, &cent_x));
282 PetscCall(DMDAVecGetArray(user->fda, user->Centy, &cent_y));
283 PetscCall(DMDAVecGetArray(user->fda, user->Centz, &cent_z));
284 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
285 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
286 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
287 cent[k][j][i].x = 0.0;
288 cent[k][j][i].y = 0.0;
289 cent[k][j][i].z = 0.0;
290 cent_x[k][j][i] = cent[k][j][i];
291 cent_y[k][j][i] = cent[k][j][i];
292 cent_z[k][j][i] = cent[k][j][i];
293 }
294 }
295 }
296 cent[1][1][1].x = 0.5 * PETSC_PI;
297 cent[1][1][1].y = 0.0;
298 cent[1][1][1].z = 0.0;
299 cent[1][2][1].x = 0.0;
300 cent[1][2][1].y = 0.5 * PETSC_PI;
301 cent[1][2][1].z = 0.0;
302 cent_z[0][1][1].x = 0.5 * PETSC_PI;
303 cent_z[0][1][1].y = 0.0;
304 cent_z[0][1][1].z = 0.0;
305 PetscCall(DMDAVecRestoreArray(user->fda, user->Centz, &cent_z));
306 PetscCall(DMDAVecRestoreArray(user->fda, user->Centy, &cent_y));
307 PetscCall(DMDAVecRestoreArray(user->fda, user->Centx, &cent_x));
308 PetscCall(DMDAVecRestoreArray(user->fda, user->Cent, &cent));
309
310 PetscCall(AnalyticalSolutionEngine(simCtx));
311
312 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
313 PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
314 PetscCall(DMDAVecGetArrayRead(user->da, user->P, &p));
315 PetscCall(PicurvAssertRealNear(vel_decay, ucat[1][1][1].x, 1.0e-12, "TGV sample should set the expected interior x velocity"));
316 PetscCall(PicurvAssertRealNear(0.0, ucat[1][1][1].y, 1.0e-12, "TGV sample should keep the paired interior y velocity at zero"));
317 PetscCall(PicurvAssertRealNear(-vel_decay, ucat[1][2][1].y, 1.0e-12, "TGV sample should set the expected interior y velocity"));
318 PetscCall(PicurvAssertRealNear(0.0, p[1][1][1], 1.0e-12, "Chosen TGV sample should produce zero pressure"));
319 PetscCall(PicurvAssertRealNear(vel_decay, ubcs[0][1][1].x, 1.0e-12, "TGV sample should set the expected boundary x velocity"));
320 PetscCall(PicurvAssertRealNear(0.0, ubcs[0][1][1].y, 1.0e-12, "TGV boundary sample should keep the paired y velocity at zero"));
321 PetscCall(DMDAVecRestoreArrayRead(user->da, user->P, &p));
322 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
323 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
324
325 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
326 PetscFunctionReturn(0);
327}
328/**
329 * @brief Tests particle analytical-solution dispatch for TGV3D, UNIFORM_FLOW, and non-analytical no-op paths.
330 */
331
332static PetscErrorCode TestAnalyticalSolutionForParticlesDispatch(void)
333{
334 SimCtx simCtx;
335 Vec tempVec = NULL;
336 PetscReal *data = NULL;
337 const PetscReal vel_decay = PetscExpReal(-0.5);
338
339 PetscFunctionBeginUser;
340 PetscCall(PetscMemzero(&simCtx, sizeof(simCtx)));
341 simCtx.ren = 2.0;
342 simCtx.ti = 0.5;
343 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "TGV3D", sizeof(simCtx.AnalyticalSolutionType)));
344
345 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 6, &tempVec));
346 PetscCall(VecGetArray(tempVec, &data));
347 data[0] = 0.5 * PETSC_PI; data[1] = 0.0; data[2] = 0.0;
348 data[3] = 0.0; data[4] = 0.5 * PETSC_PI; data[5] = 0.0;
349 PetscCall(VecRestoreArray(tempVec, &data));
350
351 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
352 PetscCall(VecGetArray(tempVec, &data));
353 PetscCall(PicurvAssertRealNear(vel_decay, data[0], 1.0e-12,
354 "TGV3D particle dispatch should populate the x velocity at x=pi/2"));
355 PetscCall(PicurvAssertRealNear(0.0, data[1], 1.0e-12,
356 "TGV3D particle dispatch should leave the first particle y velocity at zero"));
357 PetscCall(PicurvAssertRealNear(0.0, data[2], 1.0e-12,
358 "TGV3D particle dispatch should leave the first particle z velocity at zero"));
359 PetscCall(PicurvAssertRealNear(0.0, data[3], 1.0e-12,
360 "TGV3D particle dispatch should leave the second particle x velocity at zero"));
361 PetscCall(PicurvAssertRealNear(-vel_decay, data[4], 1.0e-12,
362 "TGV3D particle dispatch should populate the y velocity at y=pi/2"));
363 PetscCall(PicurvAssertRealNear(0.0, data[5], 1.0e-12,
364 "TGV3D particle dispatch should leave the second particle z velocity at zero"));
365 PetscCall(VecRestoreArray(tempVec, &data));
366
367 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
368 PetscCall(VecGetArray(tempVec, &data));
369 data[0] = 3.0;
370 data[1] = 4.0;
371 data[2] = 5.0;
372 PetscCall(VecRestoreArray(tempVec, &data));
373 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
374 PetscCall(VecGetArray(tempVec, &data));
375 PetscCall(PicurvAssertRealNear(3.0, data[0], 1.0e-12,
376 "Non-TGV particle dispatch should leave the vector untouched"));
377 PetscCall(PicurvAssertRealNear(4.0, data[1], 1.0e-12,
378 "Non-TGV particle dispatch should preserve the y component"));
379 PetscCall(PicurvAssertRealNear(5.0, data[2], 1.0e-12,
380 "Non-TGV particle dispatch should preserve the z component"));
381 PetscCall(VecRestoreArray(tempVec, &data));
382
383 simCtx.AnalyticalUniformVelocity.x = 0.125;
384 simCtx.AnalyticalUniformVelocity.y = -0.25;
385 simCtx.AnalyticalUniformVelocity.z = 0.375;
386 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
387 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
388 PetscCall(VecGetArray(tempVec, &data));
389 PetscCall(PicurvAssertRealNear(0.125, data[0], 1.0e-12,
390 "UNIFORM_FLOW particle dispatch should populate the x velocity"));
391 PetscCall(PicurvAssertRealNear(-0.25, data[1], 1.0e-12,
392 "UNIFORM_FLOW particle dispatch should populate the y velocity"));
393 PetscCall(PicurvAssertRealNear(0.375, data[2], 1.0e-12,
394 "UNIFORM_FLOW particle dispatch should populate the z velocity"));
395 PetscCall(PicurvAssertRealNear(0.125, data[3], 1.0e-12,
396 "UNIFORM_FLOW particle dispatch should use the same x velocity for each particle"));
397 PetscCall(PicurvAssertRealNear(-0.25, data[4], 1.0e-12,
398 "UNIFORM_FLOW particle dispatch should use the same y velocity for each particle"));
399 PetscCall(PicurvAssertRealNear(0.375, data[5], 1.0e-12,
400 "UNIFORM_FLOW particle dispatch should use the same z velocity for each particle"));
401 PetscCall(VecRestoreArray(tempVec, &data));
402
403 PetscCall(VecDestroy(&tempVec));
404 PetscFunctionReturn(0);
405}
406/**
407 * @brief Tests deterministic LES eddy-viscosity computation on a linear velocity field.
408 */
409
410static PetscErrorCode TestComputeEddyViscosityLESDeterministicField(void)
411{
412 SimCtx *simCtx = NULL;
413 UserCtx *user = NULL;
414 Cmpnts ***ucat = NULL;
415 PetscReal ***nu_t = NULL;
416 const PetscReal expected_nu_t = 0.25 * PetscSqrtReal(2.0);
417
418 PetscFunctionBeginUser;
419 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
420 PetscCall(DMCreateGlobalVector(user->da, &user->Nu_t));
421 PetscCall(DMCreateLocalVector(user->da, &user->lNu_t));
422 PetscCall(DMCreateGlobalVector(user->da, &user->CS));
423 PetscCall(DMCreateLocalVector(user->da, &user->lCs));
424
425 PetscCall(VecSet(user->Aj, 1.0));
426 PetscCall(VecSet(user->Nu_t, 0.0));
427 PetscCall(VecSet(user->CS, 0.5));
428 PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
429 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
430 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
431 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
432 ucat[k][j][i].x = (PetscReal)i;
433 ucat[k][j][i].y = 0.0;
434 ucat[k][j][i].z = 0.0;
435 }
436 }
437 }
438 PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
439 PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
440 PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
441 PetscCall(DMGlobalToLocalBegin(user->da, user->Aj, INSERT_VALUES, user->lAj));
442 PetscCall(DMGlobalToLocalEnd(user->da, user->Aj, INSERT_VALUES, user->lAj));
443 PetscCall(DMGlobalToLocalBegin(user->da, user->CS, INSERT_VALUES, user->lCs));
444 PetscCall(DMGlobalToLocalEnd(user->da, user->CS, INSERT_VALUES, user->lCs));
445
446 PetscCall(ComputeEddyViscosityLES(user));
447 PetscCall(DMDAVecGetArrayRead(user->da, user->Nu_t, &nu_t));
448 PetscCall(PicurvAssertRealNear(expected_nu_t, nu_t[2][2][2], 1.0e-6,
449 "linear velocity field should yield deterministic LES eddy viscosity"));
450 PetscCall(PicurvAssertRealNear(0.0, nu_t[0][2][2], 1.0e-12,
451 "boundary cells should remain untouched by the interior LES loop"));
452 PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nu_t, &nu_t));
453
454 PetscCall(VecDestroy(&user->CS));
455 PetscCall(VecDestroy(&user->lCs));
456 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
457 PetscFunctionReturn(0);
458}
459/**
460 * @brief Tests FlowSolver guardrails for unsupported momentum solver selections.
461 */
462
463static PetscErrorCode TestFlowSolverRejectsUnsupportedMomentumSolverType(void)
464{
465 SimCtx *simCtx = NULL;
466 UserCtx *user = NULL;
467 PetscErrorCode ierr_flow = 0;
468
469 PetscFunctionBeginUser;
470 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
471 simCtx->mom_solver_type = (MomentumSolverType)999;
472
473 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
474 ierr_flow = FlowSolver(simCtx);
475 PetscCall(PetscPopErrorHandler());
476 PetscCall(PicurvAssertBool((PetscBool)(ierr_flow != 0),
477 "FlowSolver should reject unsupported momentum solver selectors"));
478
479 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
480 PetscFunctionReturn(0);
481}
482/**
483 * @brief Tests driven-channel flow source-term evaluation.
484 */
485
486static PetscErrorCode TestDrivenChannelFlowSource(void)
487{
488 SimCtx *simCtx = NULL;
489 UserCtx *user = NULL;
490 Vec rct = NULL;
491 Cmpnts ***rct_arr = NULL;
492 Cmpnts ***l_csi = NULL;
493 Cmpnts ***l_eta = NULL;
494 Cmpnts ***l_zet = NULL;
495
496 PetscFunctionBeginUser;
497 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
498 PetscCall(VecDuplicate(user->Ucont, &rct));
499 PetscCall(VecZeroEntries(rct));
500
501 user->boundary_faces[BC_FACE_NEG_X].face_id = BC_FACE_NEG_X;
502 user->boundary_faces[BC_FACE_NEG_X].mathematical_type = PERIODIC;
503 user->boundary_faces[BC_FACE_NEG_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
504 user->boundary_faces[BC_FACE_POS_X].face_id = BC_FACE_POS_X;
505 user->boundary_faces[BC_FACE_POS_X].mathematical_type = PERIODIC;
506 user->boundary_faces[BC_FACE_POS_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
507 simCtx->bulkVelocityCorrection = 2.0;
508 simCtx->dt = 1.0;
509 simCtx->forceScalingFactor = 1.0;
510 simCtx->drivingForceMagnitude = 0.0;
511 PetscCall(VecSet(user->Nvert, 0.0));
512 PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
513 PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
514 PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
515 PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
516 PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
517 l_csi[1][1][1].x = 1.0; l_csi[1][1][1].y = 0.0; l_csi[1][1][1].z = 0.0;
518 l_eta[1][1][1].x = 0.0; l_eta[1][1][1].y = 1.0; l_eta[1][1][1].z = 0.0;
519 l_zet[1][1][1].x = 0.0; l_zet[1][1][1].y = 0.0; l_zet[1][1][1].z = 1.0;
520 PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
521 PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
522 PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
523
524 PetscCall(ComputeDrivenChannelFlowSource(user, rct));
525 PetscCall(DMDAVecGetArrayRead(user->fda, rct, &rct_arr));
526 PetscCall(PicurvAssertRealNear(1.5, simCtx->drivingForceMagnitude, 1.0e-12,
527 "driven flow source should update the controller magnitude"));
528 PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].y, 1.0e-12,
529 "driven flow source should leave the y component unchanged"));
530 PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].z, 1.0e-12,
531 "driven flow source should leave the z component unchanged"));
532 PetscCall(DMDAVecRestoreArrayRead(user->fda, rct, &rct_arr));
533
534 PetscCall(VecDestroy(&rct));
535 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
536 PetscFunctionReturn(0);
537}
538/**
539 * @brief Runs the unit-solver PETSc test binary.
540 */
541
542int main(int argc, char **argv)
543{
544 PetscErrorCode ierr;
545 const PicurvTestCase cases[] = {
546 {"les-filter-paths", TestLESTestFilterPaths},
547 {"analytical-geometry-selection", TestAnalyticalGeometrySelection},
548 {"analytical-scalar-verification-helpers", TestAnalyticalScalarVerificationHelpers},
549 {"analytical-solution-engine-dispatch", TestAnalyticalSolutionEngineDispatch},
550 {"analytical-solution-engine-taylor-green-samples", TestAnalyticalSolutionEngineTaylorGreenSamples},
551 {"analytical-solution-for-particles-dispatch", TestAnalyticalSolutionForParticlesDispatch},
552 {"compute-eddy-viscosity-les-deterministic-field", TestComputeEddyViscosityLESDeterministicField},
553 {"flow-solver-rejects-unsupported-momentum-solver-type", TestFlowSolverRejectsUnsupportedMomentumSolverType},
554 {"driven-channel-flow-source", TestDrivenChannelFlowSource},
555 };
556
557 ierr = PetscInitialize(&argc, &argv, NULL, "PICurv solver utility tests");
558 if (ierr) {
559 return (int)ierr;
560 }
561
562 ierr = PicurvRunTests("unit-solver", cases, sizeof(cases) / sizeof(cases[0]));
563 if (ierr) {
564 PetscFinalize();
565 return (int)ierr;
566 }
567
568 ierr = PetscFinalize();
569 return (int)ierr;
570}
SetAnalyticalScalarFieldOnParticles
PetscErrorCode SetAnalyticalScalarFieldOnParticles(UserCtx *user, const char *swarm_field_name)
Writes the configured verification scalar profile onto a particle swarm scalar field.
Definition AnalyticalSolutions.c:653
EvaluateAnalyticalScalarProfile
PetscErrorCode EvaluateAnalyticalScalarProfile(const SimCtx *simCtx, PetscReal x, PetscReal y, PetscReal z, PetscReal t, PetscReal *value)
Evaluates the configured verification scalar profile at one physical point.
Definition AnalyticalSolutions.c:627
SetAnalyticalScalarFieldAtCellCenters
PetscErrorCode SetAnalyticalScalarFieldAtCellCenters(UserCtx *user, Vec targetVec)
Writes the configured verification scalar profile at physical cell centers into a scalar Vec.
Definition AnalyticalSolutions.c:695
AnalyticalSolutionEngine
PetscErrorCode AnalyticalSolutionEngine(SimCtx *simCtx)
Dispatches to the appropriate analytical solution function based on simulation settings.
Definition AnalyticalSolutions.c:184
AnalyticalTypeRequiresCustomGeometry
PetscBool AnalyticalTypeRequiresCustomGeometry(const char *analytical_type)
Reports whether an analytical type requires custom geometry/decomposition logic.
Definition AnalyticalSolutions.c:55
SetAnalyticalSolutionForParticles
PetscErrorCode SetAnalyticalSolutionForParticles(Vec tempVec, SimCtx *simCtx)
Applies the analytical solution to particle velocity vector.
Definition AnalyticalSolutions.c:564
SetAnalyticalGridInfo
PetscErrorCode SetAnalyticalGridInfo(UserCtx *user)
Sets the grid domain and resolution for analytical solution cases.
Definition AnalyticalSolutions.c:80
ComputeDrivenChannelFlowSource
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:14
ApplyLESTestFilter
double ApplyLESTestFilter(const SimCtx *simCtx, double values[3][3][3], double weights[3][3][3])
Applies a numerical "test filter" to a 3x3x3 stencil of data points.
Definition Filter.c:123
ComputeEddyViscosityLES
PetscErrorCode ComputeEddyViscosityLES(UserCtx *user)
Computes the turbulent eddy viscosity (Nu_t) for the LES model.
Definition les.c:315
FlowSolver
PetscErrorCode FlowSolver(SimCtx *simCtx)
Orchestrates a single time step of the Eulerian fluid solver.
Definition solvers.c:11
TestAnalyticalSolutionEngineTaylorGreenSamples
static PetscErrorCode TestAnalyticalSolutionEngineTaylorGreenSamples(void)
Tests exact Taylor-Green samples on selected Eulerian interior and boundary points.
Definition test_solver_kernels.c:261
TestComputeEddyViscosityLESDeterministicField
static PetscErrorCode TestComputeEddyViscosityLESDeterministicField(void)
Tests deterministic LES eddy-viscosity computation on a linear velocity field.
Definition test_solver_kernels.c:410
TestFlowSolverRejectsUnsupportedMomentumSolverType
static PetscErrorCode TestFlowSolverRejectsUnsupportedMomentumSolverType(void)
Tests FlowSolver guardrails for unsupported momentum solver selections.
Definition test_solver_kernels.c:463
main
int main(int argc, char **argv)
Runs the unit-solver PETSc test binary.
Definition test_solver_kernels.c:542
TestAnalyticalScalarVerificationHelpers
static PetscErrorCode TestAnalyticalScalarVerificationHelpers(void)
Tests analytical scalar verification helper routines.
Definition test_solver_kernels.c:108
TestDrivenChannelFlowSource
static PetscErrorCode TestDrivenChannelFlowSource(void)
Tests driven-channel flow source-term evaluation.
Definition test_solver_kernels.c:486
TestAnalyticalGeometrySelection
static PetscErrorCode TestAnalyticalGeometrySelection(void)
Tests analytical geometry selection for supported analytical solutions.
Definition test_solver_kernels.c:57
TestAnalyticalSolutionForParticlesDispatch
static PetscErrorCode TestAnalyticalSolutionForParticlesDispatch(void)
Tests particle analytical-solution dispatch for TGV3D, UNIFORM_FLOW, and non-analytical no-op paths.
Definition test_solver_kernels.c:332
TestLESTestFilterPaths
static PetscErrorCode TestLESTestFilterPaths(void)
Tests LES test-filter helper paths for representative cases.
Definition test_solver_kernels.c:17
TestAnalyticalSolutionEngineDispatch
static PetscErrorCode TestAnalyticalSolutionEngineDispatch(void)
Tests analytical solution engine ZERO_FLOW, UNIFORM_FLOW, and unknown-type dispatch.
Definition test_solver_kernels.c:182
PicurvCreateMinimalContexts
PetscErrorCode PicurvCreateMinimalContexts(SimCtx **simCtx_out, UserCtx **user_out, PetscInt mx, PetscInt my, PetscInt mz)
Builds minimal SimCtx and UserCtx fixtures for C unit tests.
Definition test_support.c:531
PicurvAssertRealNear
PetscErrorCode PicurvAssertRealNear(PetscReal expected, PetscReal actual, PetscReal tol, const char *context)
Asserts that two real values agree within tolerance.
Definition test_support.c:906
PicurvDestroyMinimalContexts
PetscErrorCode PicurvDestroyMinimalContexts(SimCtx **simCtx_ptr, UserCtx **user_ptr)
Destroys minimal SimCtx/UserCtx fixtures and all owned PETSc objects.
Definition test_support.c:783
PicurvCreateSwarmPair
PetscErrorCode PicurvCreateSwarmPair(UserCtx *user, PetscInt nlocal, const char *post_field_name)
Creates matched solver and post-processing swarms for tests.
Definition test_support.c:699
PicurvRunTests
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.
Definition test_support.c:138
PicurvAssertVecConstant
PetscErrorCode PicurvAssertVecConstant(Vec vec, PetscScalar expected, PetscReal tol, const char *context)
Asserts that a PETSc vector is spatially constant within tolerance.
Definition test_support.c:965
PicurvAssertBool
PetscErrorCode PicurvAssertBool(PetscBool value, const char *context)
Asserts that one boolean condition is true.
Definition test_support.c:936
test_support.h
Shared declarations for the PICurv C test fixture and assertion layer.
PicurvTestCase
Named test case descriptor consumed by PicurvRunTests.
Definition test_support.h:17
PERIODIC
@ PERIODIC
Definition variables.h:260
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:848
SimCtx::block_number
PetscInt block_number
Definition variables.h:726
UserCtx::lNvert
Vec lNvert
Definition variables.h:856
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:737
UserCtx::Centz
Vec Centz
Definition variables.h:880
UserCtx::Min_X
PetscReal Min_X
Definition variables.h:840
UserCtx::lZet
Vec lZet
Definition variables.h:879
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
SimCtx::ren
PetscReal ren
Definition variables.h:699
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
UserCtx::_this
PetscInt _this
Definition variables.h:843
SimCtx::dt
PetscReal dt
Definition variables.h:666
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:739
UserCtx::Max_Y
PetscReal Max_Y
Definition variables.h:840
UserCtx::lCs
Vec lCs
Definition variables.h:886
UserCtx::Ucont
Vec Ucont
Definition variables.h:856
BCS::Ubcs
Vec Ubcs
Physical Cartesian velocity at boundary faces. Full 3D array but only boundary-face entries are meani...
Definition variables.h:121
MomentumSolverType
MomentumSolverType
Enumerator to identify the implemented momentum solver strategies.
Definition variables.h:502
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Centx
Vec Centx
Definition variables.h:880
UserCtx::Bcs
BCS Bcs
Definition variables.h:851
SimCtx::verificationScalar
VerificationScalarConfig verificationScalar
Definition variables.h:714
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:886
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:886
UserCtx::lCsi
Vec lCsi
Definition variables.h:879
SimCtx::AnalyticalUniformVelocity
Cmpnts AnalyticalUniformVelocity
Definition variables.h:706
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:856
SimCtx::AnalyticalSolutionType
char AnalyticalSolutionType[PETSC_MAX_PATH_LEN]
Definition variables.h:684
UserCtx::Max_X
PetscReal Max_X
Definition variables.h:840
UserCtx::Min_Y
PetscReal Min_Y
Definition variables.h:840
UserCtx::lAj
Vec lAj
Definition variables.h:879
SimCtx::testfilter_ik
PetscInt testfilter_ik
Definition variables.h:749
UserCtx::info
DMDALocalInfo info
Definition variables.h:837
UserCtx::lUcat
Vec lUcat
Definition variables.h:856
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:879
UserCtx::Cent
Vec Cent
Definition variables.h:879
UserCtx::Nvert
Vec Nvert
Definition variables.h:856
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
UserCtx::Centy
Vec Centy
Definition variables.h:880
SimCtx::ti
PetscReal ti
Definition variables.h:660
UserCtx::Max_Z
PetscReal Max_Z
Definition variables.h:840
SimCtx::mom_solver_type
MomentumSolverType mom_solver_type
Definition variables.h:691
SimCtx::drivingForceMagnitude
PetscReal drivingForceMagnitude
Definition variables.h:737
UserCtx::Psi
Vec Psi
Definition variables.h:904
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:651
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:830
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:303
Generated by doxygen 1.9.8