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 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
101 PetscFunctionReturn(0);
102}
103/**
104 * @brief Tests analytical scalar verification helper routines.
105 */
106
107static PetscErrorCode TestAnalyticalScalarVerificationHelpers(void)
108{
109 SimCtx *simCtx = NULL;
110 UserCtx *user = NULL;
111 Vec target = NULL;
112 PetscReal value = 0.0;
113 PetscReal ***target_arr = NULL;
114 PetscReal *positions = NULL;
115 PetscReal *psi = NULL;
116
117 PetscFunctionBeginUser;
118 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
119 PetscCall(PicurvCreateSwarmPair(user, 2, "ske"));
120 simCtx->verificationScalar.enabled = PETSC_TRUE;
121 PetscCall(PetscStrncpy(simCtx->verificationScalar.mode,
122 "analytical",
123 sizeof(simCtx->verificationScalar.mode)));
124
125 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
126 "CONSTANT",
127 sizeof(simCtx->verificationScalar.profile)));
128 simCtx->verificationScalar.value = 2.5;
129 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.2, 0.3, 0.4, 0.0, &value));
130 PetscCall(PicurvAssertRealNear(2.5, value, 1.0e-12,
131 "constant scalar profile should evaluate to the configured value"));
132
133 PetscCall(DMSwarmGetField(user->swarm, "position", NULL, NULL, (void **)&positions));
134 positions[0] = 0.1; positions[1] = 0.2; positions[2] = 0.3;
135 positions[3] = 0.8; positions[4] = 0.6; positions[5] = 0.4;
136 PetscCall(DMSwarmRestoreField(user->swarm, "position", NULL, NULL, (void **)&positions));
137
138 PetscCall(SetAnalyticalScalarFieldOnParticles(user, "Psi"));
139 PetscCall(DMSwarmGetField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
140 PetscCall(PicurvAssertRealNear(2.5, psi[0], 1.0e-12,
141 "SetAnalyticalScalarFieldOnParticles should overwrite the first particle scalar"));
142 PetscCall(PicurvAssertRealNear(2.5, psi[1], 1.0e-12,
143 "SetAnalyticalScalarFieldOnParticles should overwrite the second particle scalar"));
144 PetscCall(DMSwarmRestoreField(user->swarm, "Psi", NULL, NULL, (void **)&psi));
145
146 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
147 "LINEAR_X",
148 sizeof(simCtx->verificationScalar.profile)));
149 simCtx->verificationScalar.phi0 = 1.0;
150 simCtx->verificationScalar.slope_x = 2.0;
151 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.25, 0.0, 0.0, 0.0, &value));
152 PetscCall(PicurvAssertRealNear(1.5, value, 1.0e-12,
153 "linear-x scalar profile should evaluate phi0 + slope_x * x"));
154
155 PetscCall(VecDuplicate(user->Psi, &target));
156 PetscCall(SetAnalyticalScalarFieldAtCellCenters(user, target));
157 PetscCall(DMDAVecGetArrayRead(user->da, target, &target_arr));
158 PetscCall(PicurvAssertRealNear(1.25, target_arr[1][1][1], 1.0e-12,
159 "cell-center scalar fill should use the physical x center coordinate"));
160 PetscCall(DMDAVecRestoreArrayRead(user->da, target, &target_arr));
161 PetscCall(VecDestroy(&target));
162
163 PetscCall(PetscStrncpy(simCtx->verificationScalar.profile,
164 "SIN_PRODUCT",
165 sizeof(simCtx->verificationScalar.profile)));
166 simCtx->verificationScalar.amplitude = 3.0;
167 simCtx->verificationScalar.kx = PETSC_PI;
168 simCtx->verificationScalar.ky = PETSC_PI;
169 simCtx->verificationScalar.kz = PETSC_PI;
170 PetscCall(EvaluateAnalyticalScalarProfile(simCtx, 0.5, 0.5, 0.5, 0.0, &value));
171 PetscCall(PicurvAssertRealNear(3.0, value, 1.0e-12,
172 "sin-product scalar profile should peak at pi/2 in each coordinate"));
173
174 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
175 PetscFunctionReturn(0);
176}
177/**
178 * @brief Tests analytical solution engine ZERO_FLOW, UNIFORM_FLOW, and unknown-type dispatch.
179 */
180
181static PetscErrorCode TestAnalyticalSolutionEngineDispatch(void)
182{
183 SimCtx *simCtx = NULL;
184 UserCtx *user = NULL;
185 PetscErrorCode ierr_unknown = 0;
186 Cmpnts ***ucat = NULL;
187 Cmpnts ***ubcs = NULL;
188
189 PetscFunctionBeginUser;
190 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
191 PetscCall(VecSet(user->Ucat, 3.0));
192 PetscCall(VecSet(user->P, 5.0));
193 PetscCall(VecSet(user->Bcs.Ubcs, 7.0));
194
195 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
196 PetscCall(AnalyticalSolutionEngine(simCtx));
197 PetscCall(PicurvAssertVecConstant(user->Ucat, 0.0, 1.0e-12, "ZERO_FLOW should zero the Eulerian velocity field"));
198 PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "ZERO_FLOW should zero the pressure field"));
199 PetscCall(PicurvAssertVecConstant(user->Bcs.Ubcs, 0.0, 1.0e-12, "ZERO_FLOW should zero boundary-condition velocity data"));
200
201 /* UNIFORM_FLOW now works in curvilinear form (sets Ucont via metric dot products,
202 derives Ucat via Contra2Cart). The test grid needs identity metrics so the
203 transformation is invertible and Ucat recovers the physical velocity exactly. */
204 {
205 Cmpnts ***l_csi, ***l_eta, ***l_zet;
206 DMDALocalInfo linfo;
207 PetscCall(DMDAGetLocalInfo(user->fda, &linfo));
208 PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
209 PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
210 PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
211 for (PetscInt k = linfo.zs; k < linfo.zs + linfo.zm; k++)
212 for (PetscInt j = linfo.ys; j < linfo.ys + linfo.ym; j++)
213 for (PetscInt i = linfo.xs; i < linfo.xs + linfo.xm; i++) {
214 l_csi[k][j][i] = (Cmpnts){1.0, 0.0, 0.0};
215 l_eta[k][j][i] = (Cmpnts){0.0, 1.0, 0.0};
216 l_zet[k][j][i] = (Cmpnts){0.0, 0.0, 1.0};
217 }
218 PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
219 PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
220 PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
221 }
222
223 simCtx->AnalyticalUniformVelocity.x = 1.25;
224 simCtx->AnalyticalUniformVelocity.y = -0.5;
225 simCtx->AnalyticalUniformVelocity.z = 0.75;
226 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx->AnalyticalSolutionType)));
227 PetscCall(AnalyticalSolutionEngine(simCtx));
228 PetscCall(PicurvAssertVecConstant(user->P, 0.0, 1.0e-12, "UNIFORM_FLOW should keep the pressure field zero"));
229
230 Cmpnts ***ucont = NULL;
231 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
232 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucont, &ucont));
233 PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
234 PetscCall(PicurvAssertRealNear(1.25, ucat[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should impose the configured x velocity"));
235 PetscCall(PicurvAssertRealNear(-0.5, ucat[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should impose the configured y velocity"));
236 PetscCall(PicurvAssertRealNear(0.75, ucat[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should impose the configured z velocity"));
237 PetscCall(PicurvAssertRealNear(1.25, ucont[1][1][1].x, 1.0e-12, "UNIFORM_FLOW should set contravariant x flux (identity metric)"));
238 PetscCall(PicurvAssertRealNear(-0.5, ucont[1][1][1].y, 1.0e-12, "UNIFORM_FLOW should set contravariant y flux (identity metric)"));
239 PetscCall(PicurvAssertRealNear(0.75, ucont[1][1][1].z, 1.0e-12, "UNIFORM_FLOW should set contravariant z flux (identity metric)"));
240 PetscCall(PicurvAssertRealNear(1.25, ubcs[0][1][1].x, 1.0e-12, "UNIFORM_FLOW should populate boundary x velocity"));
241 PetscCall(PicurvAssertRealNear(-0.5, ubcs[0][1][1].y, 1.0e-12, "UNIFORM_FLOW should populate boundary y velocity"));
242 PetscCall(PicurvAssertRealNear(0.75, ubcs[0][1][1].z, 1.0e-12, "UNIFORM_FLOW should populate boundary z velocity"));
243 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
244 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucont, &ucont));
245 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
246
247 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "NOT_A_REAL_ANALYTICAL_TYPE", sizeof(simCtx->AnalyticalSolutionType)));
248 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
249 ierr_unknown = AnalyticalSolutionEngine(simCtx);
250 PetscCall(PetscPopErrorHandler());
251 PetscCall(PicurvAssertBool((PetscBool)(ierr_unknown != 0),
252 "AnalyticalSolutionEngine should reject unknown analytical type strings"));
253
254 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
255 PetscFunctionReturn(0);
256}
257/**
258 * @brief Tests exact Taylor-Green samples on selected Eulerian interior and boundary points.
259 */
260static PetscErrorCode TestAnalyticalSolutionEngineTaylorGreenSamples(void)
261{
262 SimCtx *simCtx = NULL;
263 UserCtx *user = NULL;
264 Cmpnts ***cent = NULL;
265 Cmpnts ***cent_x = NULL;
266 Cmpnts ***cent_y = NULL;
267 Cmpnts ***cent_z = NULL;
268 Cmpnts ***ucat = NULL;
269 Cmpnts ***ubcs = NULL;
270 PetscReal ***p = NULL;
271 const PetscReal vel_decay = PetscExpReal(-0.5);
272
273 PetscFunctionBeginUser;
274 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
275 simCtx->ren = 2.0;
276 simCtx->ti = 0.5;
277 PetscCall(PetscStrncpy(simCtx->AnalyticalSolutionType, "TGV3D", sizeof(simCtx->AnalyticalSolutionType)));
278
279 PetscCall(DMDAVecGetArray(user->fda, user->Cent, &cent));
280 PetscCall(DMDAVecGetArray(user->fda, user->Centx, &cent_x));
281 PetscCall(DMDAVecGetArray(user->fda, user->Centy, &cent_y));
282 PetscCall(DMDAVecGetArray(user->fda, user->Centz, &cent_z));
283 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
284 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
285 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
286 cent[k][j][i].x = 0.0;
287 cent[k][j][i].y = 0.0;
288 cent[k][j][i].z = 0.0;
289 cent_x[k][j][i] = cent[k][j][i];
290 cent_y[k][j][i] = cent[k][j][i];
291 cent_z[k][j][i] = cent[k][j][i];
292 }
293 }
294 }
295 cent[1][1][1].x = 0.5 * PETSC_PI;
296 cent[1][1][1].y = 0.0;
297 cent[1][1][1].z = 0.0;
298 cent[1][2][1].x = 0.0;
299 cent[1][2][1].y = 0.5 * PETSC_PI;
300 cent[1][2][1].z = 0.0;
301 cent_z[0][1][1].x = 0.5 * PETSC_PI;
302 cent_z[0][1][1].y = 0.0;
303 cent_z[0][1][1].z = 0.0;
304 PetscCall(DMDAVecRestoreArray(user->fda, user->Centz, &cent_z));
305 PetscCall(DMDAVecRestoreArray(user->fda, user->Centy, &cent_y));
306 PetscCall(DMDAVecRestoreArray(user->fda, user->Centx, &cent_x));
307 PetscCall(DMDAVecRestoreArray(user->fda, user->Cent, &cent));
308
309 PetscCall(AnalyticalSolutionEngine(simCtx));
310
311 PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
312 PetscCall(DMDAVecGetArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
313 PetscCall(DMDAVecGetArrayRead(user->da, user->P, &p));
314 PetscCall(PicurvAssertRealNear(vel_decay, ucat[1][1][1].x, 1.0e-12, "TGV sample should set the expected interior x velocity"));
315 PetscCall(PicurvAssertRealNear(0.0, ucat[1][1][1].y, 1.0e-12, "TGV sample should keep the paired interior y velocity at zero"));
316 PetscCall(PicurvAssertRealNear(-vel_decay, ucat[1][2][1].y, 1.0e-12, "TGV sample should set the expected interior y velocity"));
317 PetscCall(PicurvAssertRealNear(0.0, p[1][1][1], 1.0e-12, "Chosen TGV sample should produce zero pressure"));
318 PetscCall(PicurvAssertRealNear(vel_decay, ubcs[0][1][1].x, 1.0e-12, "TGV sample should set the expected boundary x velocity"));
319 PetscCall(PicurvAssertRealNear(0.0, ubcs[0][1][1].y, 1.0e-12, "TGV boundary sample should keep the paired y velocity at zero"));
320 PetscCall(DMDAVecRestoreArrayRead(user->da, user->P, &p));
321 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Bcs.Ubcs, &ubcs));
322 PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
323
324 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
325 PetscFunctionReturn(0);
326}
327/**
328 * @brief Tests particle analytical-solution dispatch for TGV3D, UNIFORM_FLOW, and non-analytical no-op paths.
329 */
330
331static PetscErrorCode TestAnalyticalSolutionForParticlesDispatch(void)
332{
333 SimCtx simCtx;
334 Vec tempVec = NULL;
335 PetscReal *data = NULL;
336 const PetscReal vel_decay = PetscExpReal(-0.5);
337
338 PetscFunctionBeginUser;
339 PetscCall(PetscMemzero(&simCtx, sizeof(simCtx)));
340 simCtx.ren = 2.0;
341 simCtx.ti = 0.5;
342 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "TGV3D", sizeof(simCtx.AnalyticalSolutionType)));
343
344 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 6, &tempVec));
345 PetscCall(VecGetArray(tempVec, &data));
346 data[0] = 0.5 * PETSC_PI; data[1] = 0.0; data[2] = 0.0;
347 data[3] = 0.0; data[4] = 0.5 * PETSC_PI; data[5] = 0.0;
348 PetscCall(VecRestoreArray(tempVec, &data));
349
350 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
351 PetscCall(VecGetArray(tempVec, &data));
352 PetscCall(PicurvAssertRealNear(vel_decay, data[0], 1.0e-12,
353 "TGV3D particle dispatch should populate the x velocity at x=pi/2"));
354 PetscCall(PicurvAssertRealNear(0.0, data[1], 1.0e-12,
355 "TGV3D particle dispatch should leave the first particle y velocity at zero"));
356 PetscCall(PicurvAssertRealNear(0.0, data[2], 1.0e-12,
357 "TGV3D particle dispatch should leave the first particle z velocity at zero"));
358 PetscCall(PicurvAssertRealNear(0.0, data[3], 1.0e-12,
359 "TGV3D particle dispatch should leave the second particle x velocity at zero"));
360 PetscCall(PicurvAssertRealNear(-vel_decay, data[4], 1.0e-12,
361 "TGV3D particle dispatch should populate the y velocity at y=pi/2"));
362 PetscCall(PicurvAssertRealNear(0.0, data[5], 1.0e-12,
363 "TGV3D particle dispatch should leave the second particle z velocity at zero"));
364 PetscCall(VecRestoreArray(tempVec, &data));
365
366 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "ZERO_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
367 PetscCall(VecGetArray(tempVec, &data));
368 data[0] = 3.0;
369 data[1] = 4.0;
370 data[2] = 5.0;
371 PetscCall(VecRestoreArray(tempVec, &data));
372 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
373 PetscCall(VecGetArray(tempVec, &data));
374 PetscCall(PicurvAssertRealNear(3.0, data[0], 1.0e-12,
375 "Non-TGV particle dispatch should leave the vector untouched"));
376 PetscCall(PicurvAssertRealNear(4.0, data[1], 1.0e-12,
377 "Non-TGV particle dispatch should preserve the y component"));
378 PetscCall(PicurvAssertRealNear(5.0, data[2], 1.0e-12,
379 "Non-TGV particle dispatch should preserve the z component"));
380 PetscCall(VecRestoreArray(tempVec, &data));
381
382 simCtx.AnalyticalUniformVelocity.x = 0.125;
383 simCtx.AnalyticalUniformVelocity.y = -0.25;
384 simCtx.AnalyticalUniformVelocity.z = 0.375;
385 PetscCall(PetscStrncpy(simCtx.AnalyticalSolutionType, "UNIFORM_FLOW", sizeof(simCtx.AnalyticalSolutionType)));
386 PetscCall(SetAnalyticalSolutionForParticles(tempVec, &simCtx));
387 PetscCall(VecGetArray(tempVec, &data));
388 PetscCall(PicurvAssertRealNear(0.125, data[0], 1.0e-12,
389 "UNIFORM_FLOW particle dispatch should populate the x velocity"));
390 PetscCall(PicurvAssertRealNear(-0.25, data[1], 1.0e-12,
391 "UNIFORM_FLOW particle dispatch should populate the y velocity"));
392 PetscCall(PicurvAssertRealNear(0.375, data[2], 1.0e-12,
393 "UNIFORM_FLOW particle dispatch should populate the z velocity"));
394 PetscCall(PicurvAssertRealNear(0.125, data[3], 1.0e-12,
395 "UNIFORM_FLOW particle dispatch should use the same x velocity for each particle"));
396 PetscCall(PicurvAssertRealNear(-0.25, data[4], 1.0e-12,
397 "UNIFORM_FLOW particle dispatch should use the same y velocity for each particle"));
398 PetscCall(PicurvAssertRealNear(0.375, data[5], 1.0e-12,
399 "UNIFORM_FLOW particle dispatch should use the same z velocity for each particle"));
400 PetscCall(VecRestoreArray(tempVec, &data));
401
402 PetscCall(VecDestroy(&tempVec));
403 PetscFunctionReturn(0);
404}
405/**
406 * @brief Tests deterministic LES eddy-viscosity computation on a linear velocity field.
407 */
408
409static PetscErrorCode TestComputeEddyViscosityLESDeterministicField(void)
410{
411 SimCtx *simCtx = NULL;
412 UserCtx *user = NULL;
413 Cmpnts ***ucat = NULL;
414 PetscReal ***nu_t = NULL;
415 const PetscReal expected_nu_t = 0.25 * PetscSqrtReal(2.0);
416
417 PetscFunctionBeginUser;
418 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 5, 5, 5));
419 PetscCall(DMCreateGlobalVector(user->da, &user->Nu_t));
420 PetscCall(DMCreateLocalVector(user->da, &user->lNu_t));
421 PetscCall(DMCreateGlobalVector(user->da, &user->CS));
422 PetscCall(DMCreateLocalVector(user->da, &user->lCs));
423
424 PetscCall(VecSet(user->Aj, 1.0));
425 PetscCall(VecSet(user->Nu_t, 0.0));
426 PetscCall(VecSet(user->CS, 0.5));
427 PetscCall(DMDAVecGetArray(user->fda, user->Ucat, &ucat));
428 for (PetscInt k = user->info.zs; k < user->info.zs + user->info.zm; ++k) {
429 for (PetscInt j = user->info.ys; j < user->info.ys + user->info.ym; ++j) {
430 for (PetscInt i = user->info.xs; i < user->info.xs + user->info.xm; ++i) {
431 ucat[k][j][i].x = (PetscReal)i;
432 ucat[k][j][i].y = 0.0;
433 ucat[k][j][i].z = 0.0;
434 }
435 }
436 }
437 PetscCall(DMDAVecRestoreArray(user->fda, user->Ucat, &ucat));
438 PetscCall(DMGlobalToLocalBegin(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
439 PetscCall(DMGlobalToLocalEnd(user->fda, user->Ucat, INSERT_VALUES, user->lUcat));
440 PetscCall(DMGlobalToLocalBegin(user->da, user->Aj, INSERT_VALUES, user->lAj));
441 PetscCall(DMGlobalToLocalEnd(user->da, user->Aj, INSERT_VALUES, user->lAj));
442 PetscCall(DMGlobalToLocalBegin(user->da, user->CS, INSERT_VALUES, user->lCs));
443 PetscCall(DMGlobalToLocalEnd(user->da, user->CS, INSERT_VALUES, user->lCs));
444
445 PetscCall(ComputeEddyViscosityLES(user));
446 PetscCall(DMDAVecGetArrayRead(user->da, user->Nu_t, &nu_t));
447 PetscCall(PicurvAssertRealNear(expected_nu_t, nu_t[2][2][2], 1.0e-6,
448 "linear velocity field should yield deterministic LES eddy viscosity"));
449 PetscCall(PicurvAssertRealNear(0.0, nu_t[0][2][2], 1.0e-12,
450 "boundary cells should remain untouched by the interior LES loop"));
451 PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nu_t, &nu_t));
452
453 PetscCall(VecDestroy(&user->CS));
454 PetscCall(VecDestroy(&user->lCs));
455 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
456 PetscFunctionReturn(0);
457}
458/**
459 * @brief Tests FlowSolver guardrails for unsupported momentum solver selections.
460 */
461
462static PetscErrorCode TestFlowSolverRejectsUnsupportedMomentumSolverType(void)
463{
464 SimCtx *simCtx = NULL;
465 UserCtx *user = NULL;
466 PetscErrorCode ierr_flow = 0;
467
468 PetscFunctionBeginUser;
469 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
470 simCtx->mom_solver_type = (MomentumSolverType)999;
471
472 PetscCall(PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL));
473 ierr_flow = FlowSolver(simCtx);
474 PetscCall(PetscPopErrorHandler());
475 PetscCall(PicurvAssertBool((PetscBool)(ierr_flow != 0),
476 "FlowSolver should reject unsupported momentum solver selectors"));
477
478 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
479 PetscFunctionReturn(0);
480}
481/**
482 * @brief Tests driven-channel flow source-term evaluation.
483 */
484
485static PetscErrorCode TestDrivenChannelFlowSource(void)
486{
487 SimCtx *simCtx = NULL;
488 UserCtx *user = NULL;
489 Vec rct = NULL;
490 Cmpnts ***rct_arr = NULL;
491 Cmpnts ***l_csi = NULL;
492 Cmpnts ***l_eta = NULL;
493 Cmpnts ***l_zet = NULL;
494
495 PetscFunctionBeginUser;
496 PetscCall(PicurvCreateMinimalContexts(&simCtx, &user, 4, 4, 4));
497 PetscCall(VecDuplicate(user->Ucont, &rct));
498 PetscCall(VecZeroEntries(rct));
499
500 user->boundary_faces[BC_FACE_NEG_X].face_id = BC_FACE_NEG_X;
501 user->boundary_faces[BC_FACE_NEG_X].mathematical_type = PERIODIC;
502 user->boundary_faces[BC_FACE_NEG_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
503 user->boundary_faces[BC_FACE_POS_X].face_id = BC_FACE_POS_X;
504 user->boundary_faces[BC_FACE_POS_X].mathematical_type = PERIODIC;
505 user->boundary_faces[BC_FACE_POS_X].handler_type = BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX;
506 simCtx->bulkVelocityCorrection = 2.0;
507 simCtx->dt = 1.0;
508 simCtx->forceScalingFactor = 1.0;
509 simCtx->drivingForceMagnitude = 0.0;
510 PetscCall(VecSet(user->Nvert, 0.0));
511 PetscCall(DMGlobalToLocalBegin(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
512 PetscCall(DMGlobalToLocalEnd(user->da, user->Nvert, INSERT_VALUES, user->lNvert));
513 PetscCall(DMDAVecGetArray(user->fda, user->lCsi, &l_csi));
514 PetscCall(DMDAVecGetArray(user->fda, user->lEta, &l_eta));
515 PetscCall(DMDAVecGetArray(user->fda, user->lZet, &l_zet));
516 l_csi[1][1][1].x = 1.0; l_csi[1][1][1].y = 0.0; l_csi[1][1][1].z = 0.0;
517 l_eta[1][1][1].x = 0.0; l_eta[1][1][1].y = 1.0; l_eta[1][1][1].z = 0.0;
518 l_zet[1][1][1].x = 0.0; l_zet[1][1][1].y = 0.0; l_zet[1][1][1].z = 1.0;
519 PetscCall(DMDAVecRestoreArray(user->fda, user->lZet, &l_zet));
520 PetscCall(DMDAVecRestoreArray(user->fda, user->lEta, &l_eta));
521 PetscCall(DMDAVecRestoreArray(user->fda, user->lCsi, &l_csi));
522
523 PetscCall(ComputeDrivenChannelFlowSource(user, rct));
524 PetscCall(DMDAVecGetArrayRead(user->fda, rct, &rct_arr));
525 PetscCall(PicurvAssertRealNear(1.5, simCtx->drivingForceMagnitude, 1.0e-12,
526 "driven flow source should update the controller magnitude"));
527 PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].y, 1.0e-12,
528 "driven flow source should leave the y component unchanged"));
529 PetscCall(PicurvAssertRealNear(0.0, rct_arr[1][1][1].z, 1.0e-12,
530 "driven flow source should leave the z component unchanged"));
531 PetscCall(DMDAVecRestoreArrayRead(user->fda, rct, &rct_arr));
532
533 PetscCall(VecDestroy(&rct));
534 PetscCall(PicurvDestroyMinimalContexts(&simCtx, &user));
535 PetscFunctionReturn(0);
536}
537/**
538 * @brief Runs the unit-solver PETSc test binary.
539 */
540
541int main(int argc, char **argv)
542{
543 PetscErrorCode ierr;
544 const PicurvTestCase cases[] = {
545 {"les-filter-paths", TestLESTestFilterPaths},
546 {"analytical-geometry-selection", TestAnalyticalGeometrySelection},
547 {"analytical-scalar-verification-helpers", TestAnalyticalScalarVerificationHelpers},
548 {"analytical-solution-engine-dispatch", TestAnalyticalSolutionEngineDispatch},
549 {"analytical-solution-engine-taylor-green-samples", TestAnalyticalSolutionEngineTaylorGreenSamples},
550 {"analytical-solution-for-particles-dispatch", TestAnalyticalSolutionForParticlesDispatch},
551 {"compute-eddy-viscosity-les-deterministic-field", TestComputeEddyViscosityLESDeterministicField},
552 {"flow-solver-rejects-unsupported-momentum-solver-type", TestFlowSolverRejectsUnsupportedMomentumSolverType},
553 {"driven-channel-flow-source", TestDrivenChannelFlowSource},
554 };
555
556 ierr = PetscInitialize(&argc, &argv, NULL, "PICurv solver utility tests");
557 if (ierr) {
558 return (int)ierr;
559 }
560
561 ierr = PicurvRunTests("unit-solver", cases, sizeof(cases) / sizeof(cases[0]));
562 if (ierr) {
563 PetscFinalize();
564 return (int)ierr;
565 }
566
567 ierr = PetscFinalize();
568 return (int)ierr;
569}
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:260
TestComputeEddyViscosityLESDeterministicField
static PetscErrorCode TestComputeEddyViscosityLESDeterministicField(void)
Tests deterministic LES eddy-viscosity computation on a linear velocity field.
Definition test_solver_kernels.c:409
TestFlowSolverRejectsUnsupportedMomentumSolverType
static PetscErrorCode TestFlowSolverRejectsUnsupportedMomentumSolverType(void)
Tests FlowSolver guardrails for unsupported momentum solver selections.
Definition test_solver_kernels.c:462
main
int main(int argc, char **argv)
Runs the unit-solver PETSc test binary.
Definition test_solver_kernels.c:541
TestAnalyticalScalarVerificationHelpers
static PetscErrorCode TestAnalyticalScalarVerificationHelpers(void)
Tests analytical scalar verification helper routines.
Definition test_solver_kernels.c:107
TestDrivenChannelFlowSource
static PetscErrorCode TestDrivenChannelFlowSource(void)
Tests driven-channel flow source-term evaluation.
Definition test_solver_kernels.c:485
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:331
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:181
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:530
PicurvAssertRealNear
PetscErrorCode PicurvAssertRealNear(PetscReal expected, PetscReal actual, PetscReal tol, const char *context)
Asserts that two real values agree within tolerance.
Definition test_support.c:902
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:698
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:961
PicurvAssertBool
PetscErrorCode PicurvAssertBool(PetscBool value, const char *context)
Asserts that one boolean condition is true.
Definition test_support.c:932
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:829
SimCtx::block_number
PetscInt block_number
Definition variables.h:712
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
SimCtx::forceScalingFactor
PetscReal forceScalingFactor
Definition variables.h:723
UserCtx::Centz
Vec Centz
Definition variables.h:859
UserCtx::Min_X
PetscReal Min_X
Definition variables.h:821
UserCtx::lZet
Vec lZet
Definition variables.h:858
BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
@ BC_HANDLER_PERIODIC_DRIVEN_CONSTANT_FLUX
Definition variables.h:286
SimCtx::ren
PetscReal ren
Definition variables.h:691
BoundaryFaceConfig::handler_type
BCHandlerType handler_type
Definition variables.h:337
UserCtx::_this
PetscInt _this
Definition variables.h:824
SimCtx::dt
PetscReal dt
Definition variables.h:658
SimCtx::bulkVelocityCorrection
PetscReal bulkVelocityCorrection
Definition variables.h:725
UserCtx::Max_Y
PetscReal Max_Y
Definition variables.h:821
UserCtx::lCs
Vec lCs
Definition variables.h:865
UserCtx::Ucont
Vec Ucont
Definition variables.h:837
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:859
UserCtx::Bcs
BCS Bcs
Definition variables.h:832
SimCtx::verificationScalar
VerificationScalarConfig verificationScalar
Definition variables.h:700
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:865
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:865
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
SimCtx::AnalyticalUniformVelocity
Cmpnts AnalyticalUniformVelocity
Definition variables.h:698
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:837
SimCtx::AnalyticalSolutionType
char AnalyticalSolutionType[PETSC_MAX_PATH_LEN]
Definition variables.h:676
UserCtx::Max_X
PetscReal Max_X
Definition variables.h:821
UserCtx::Min_Y
PetscReal Min_Y
Definition variables.h:821
UserCtx::lAj
Vec lAj
Definition variables.h:858
SimCtx::testfilter_ik
PetscInt testfilter_ik
Definition variables.h:735
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:858
UserCtx::Cent
Vec Cent
Definition variables.h:858
UserCtx::Nvert
Vec Nvert
Definition variables.h:837
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
UserCtx::Centy
Vec Centy
Definition variables.h:859
SimCtx::ti
PetscReal ti
Definition variables.h:652
UserCtx::Max_Z
PetscReal Max_Z
Definition variables.h:821
SimCtx::mom_solver_type
MomentumSolverType mom_solver_type
Definition variables.h:683
SimCtx::drivingForceMagnitude
PetscReal drivingForceMagnitude
Definition variables.h:723
UserCtx::Psi
Vec Psi
Definition variables.h:883
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:643
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:811
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:303
Generated by doxygen 1.9.8