1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Time-stepping functions for Navier-Stokes example using PETSc
10
11 #include <ceed.h>
12 #include <petscdmplex.h>
13 #include <petscts.h>
14
15 #include "../navierstokes.h"
16 #include "../qfunctions/newtonian_state.h"
17
18 // Insert Boundary values if it's a new time
UpdateBoundaryValues(User user,Vec Q_loc,PetscReal t)19 PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t) {
20 PetscFunctionBeginUser;
21 if (user->time_bc_set != t) {
22 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, NULL, NULL, NULL));
23 user->time_bc_set = t;
24 }
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
28 // RHS (Explicit time-stepper) function setup
29 // This is the RHS of the ODE, given as u_t = G(t,u)
30 // This function takes in a state vector Q and writes into G
RHS_NS(TS ts,PetscReal t,Vec Q,Vec G,void * user_data)31 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) {
32 User user = *(User *)user_data;
33 Ceed ceed = user->ceed;
34 PetscScalar dt;
35 Vec Q_loc = user->Q_loc;
36 PetscMemType q_mem_type;
37
38 PetscFunctionBeginUser;
39 // Update time dependent data
40 PetscCall(UpdateBoundaryValues(user, Q_loc, t));
41 if (user->phys->solution_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_rhs_ctx->op, user->phys->solution_time_label, &t));
42 PetscCall(TSGetTimeStep(ts, &dt));
43 if (user->phys->timestep_size_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_rhs_ctx->op, user->phys->timestep_size_label, &dt));
44
45 PetscCall(ApplyCeedOperatorGlobalToGlobal(Q, G, user->op_rhs_ctx));
46
47 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, user->q_ceed));
48
49 // Inverse of the mass matrix
50 PetscCall(KSPSolve(user->mass_ksp, G, G));
51
52 PetscCall(VecReadCeedToPetsc(user->q_ceed, q_mem_type, Q_loc));
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55
56 // Surface forces function setup
Surface_Forces_NS(DM dm,Vec G_loc,PetscInt num_walls,const PetscInt walls[],PetscScalar * reaction_force)57 static PetscErrorCode Surface_Forces_NS(DM dm, Vec G_loc, PetscInt num_walls, const PetscInt walls[], PetscScalar *reaction_force) {
58 DMLabel face_label;
59 const PetscScalar *g;
60 PetscInt dof, dim = 3;
61 MPI_Comm comm;
62 PetscSection s;
63
64 PetscFunctionBeginUser;
65 PetscCall(PetscArrayzero(reaction_force, num_walls * dim));
66 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
67 PetscCall(DMGetLabel(dm, "Face Sets", &face_label));
68 PetscCall(VecGetArrayRead(G_loc, &g));
69 for (PetscInt w = 0; w < num_walls; w++) {
70 const PetscInt wall = walls[w];
71 IS wall_is;
72 PetscCall(DMGetLocalSection(dm, &s));
73 PetscCall(DMLabelGetStratumIS(face_label, wall, &wall_is));
74 if (wall_is) { // There exist such points on this process
75 PetscInt num_points;
76 PetscInt num_comp = 0;
77 const PetscInt *points;
78 PetscCall(PetscSectionGetFieldComponents(s, 0, &num_comp));
79 PetscCall(ISGetSize(wall_is, &num_points));
80 PetscCall(ISGetIndices(wall_is, &points));
81 for (PetscInt i = 0; i < num_points; i++) {
82 const PetscInt p = points[i];
83 const StateConservative *r;
84 PetscCall(DMPlexPointLocalRead(dm, p, g, &r));
85 PetscCall(PetscSectionGetDof(s, p, &dof));
86 for (PetscInt node = 0; node < dof / num_comp; node++) {
87 for (PetscInt j = 0; j < 3; j++) {
88 reaction_force[w * dim + j] -= r[node].momentum[j];
89 }
90 }
91 }
92 PetscCall(ISRestoreIndices(wall_is, &points));
93 }
94 PetscCall(ISDestroy(&wall_is));
95 }
96 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, reaction_force, dim * num_walls, MPIU_SCALAR, MPI_SUM, comm));
97 // Restore Vectors
98 PetscCall(VecRestoreArrayRead(G_loc, &g));
99 PetscFunctionReturn(PETSC_SUCCESS);
100 }
101
102 // Implicit time-stepper function setup
IFunction_NS(TS ts,PetscReal t,Vec Q,Vec Q_dot,Vec G,void * user_data)103 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, void *user_data) {
104 User user = *(User *)user_data;
105 Ceed ceed = user->ceed;
106 PetscScalar dt;
107 Vec Q_loc = user->Q_loc, Q_dot_loc = user->Q_dot_loc, G_loc;
108 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type;
109
110 PetscFunctionBeginUser;
111 // Get local vectors
112 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
113
114 // Update time dependent data
115 PetscCall(UpdateBoundaryValues(user, Q_loc, t));
116 if (user->phys->solution_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_ifunction, user->phys->solution_time_label, &t));
117 PetscCall(TSGetTimeStep(ts, &dt));
118 if (user->phys->timestep_size_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_ifunction, user->phys->timestep_size_label, &dt));
119
120 // Global-to-local
121 PetscCall(DMGlobalToLocalBegin(user->dm, Q, INSERT_VALUES, Q_loc));
122 PetscCall(DMGlobalToLocalBegin(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc));
123 PetscCall(DMGlobalToLocalEnd(user->dm, Q, INSERT_VALUES, Q_loc));
124 PetscCall(DMGlobalToLocalEnd(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc));
125
126 // Place PETSc vectors in CEED vectors
127 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, user->q_ceed));
128 PetscCall(VecReadPetscToCeed(Q_dot_loc, &q_dot_mem_type, user->q_dot_ceed));
129 PetscCall(VecPetscToCeed(G_loc, &g_mem_type, user->g_ceed));
130
131 // Apply CEED operator
132 PetscCall(PetscLogEventBegin(FLUIDS_CeedOperatorApply, Q, G, 0, 0));
133 PetscCall(PetscLogGpuTimeBegin());
134 PetscCallCeed(user->ceed, CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE));
135 PetscCall(PetscLogGpuTimeEnd());
136 PetscCall(PetscLogEventEnd(FLUIDS_CeedOperatorApply, Q, G, 0, 0));
137
138 // Restore vectors
139 PetscCall(VecReadCeedToPetsc(user->q_ceed, q_mem_type, Q_loc));
140 PetscCall(VecReadCeedToPetsc(user->q_dot_ceed, q_dot_mem_type, Q_dot_loc));
141 PetscCall(VecCeedToPetsc(user->g_ceed, g_mem_type, G_loc));
142
143 // Local-to-Global
144 PetscCall(VecZeroEntries(G));
145 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G));
146
147 // Restore vectors
148 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
FormIJacobian_NS(TS ts,PetscReal t,Vec Q,Vec Q_dot,PetscReal shift,Mat J,Mat J_pre,void * user_data)152 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) {
153 User user = *(User *)user_data;
154 PetscBool J_is_matceed, J_is_mffd, J_pre_is_matceed, J_pre_is_mffd;
155
156 PetscFunctionBeginUser;
157 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd));
158 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATCEED, &J_is_matceed));
159 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATMFFD, &J_pre_is_mffd));
160 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATCEED, &J_pre_is_matceed));
161
162 PetscCall(MatCeedSetContextReal(user->mat_ijacobian, "ijacobian time shift", shift));
163
164 if (J_is_matceed || J_is_mffd) {
165 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
166 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
167 } else PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J));
168
169 if (J_pre_is_matceed && J != J_pre) {
170 PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY));
171 PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY));
172 } else if (!J_pre_is_matceed && !J_pre_is_mffd && J != J_pre) {
173 PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J_pre));
174 }
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
WriteOutput(User user,Vec Q,PetscInt step_no,PetscScalar time)178 PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) {
179 Vec Q_loc;
180 char file_path[PETSC_MAX_PATH_LEN];
181 PetscViewer viewer;
182
183 PetscFunctionBeginUser;
184 if (user->app_ctx->checkpoint_vtk) {
185 // Set up output
186 PetscCall(DMGetLocalVector(user->dm, &Q_loc));
187 PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec"));
188 PetscCall(VecZeroEntries(Q_loc));
189 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc));
190
191 // Output
192 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no));
193
194 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer));
195 PetscCall(VecView(Q_loc, viewer));
196 PetscCall(PetscViewerDestroy(&viewer));
197 if (user->dm_viz) {
198 Vec Q_refined, Q_refined_loc;
199 char file_path_refined[PETSC_MAX_PATH_LEN];
200 PetscViewer viewer_refined;
201
202 PetscCall(DMGetGlobalVector(user->dm_viz, &Q_refined));
203 PetscCall(DMGetLocalVector(user->dm_viz, &Q_refined_loc));
204 PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"));
205
206 PetscCall(MatInterpolate(user->interp_viz, Q, Q_refined));
207 PetscCall(VecZeroEntries(Q_refined_loc));
208 PetscCall(DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc));
209
210 PetscCall(PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir,
211 step_no));
212
213 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined));
214 PetscCall(VecView(Q_refined_loc, viewer_refined));
215 PetscCall(DMRestoreLocalVector(user->dm_viz, &Q_refined_loc));
216 PetscCall(DMRestoreGlobalVector(user->dm_viz, &Q_refined));
217 PetscCall(PetscViewerDestroy(&viewer_refined));
218 }
219 PetscCall(DMRestoreLocalVector(user->dm, &Q_loc));
220 }
221
222 // Save data in a binary file for continuation of simulations
223 if (user->app_ctx->add_stepnum2bin) {
224 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", user->app_ctx->output_dir, step_no));
225 } else {
226 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", user->app_ctx->output_dir));
227 }
228 PetscCall(PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer));
229
230 PetscInt32 token = PetscDefined(USE_64BIT_INDICES) ? FLUIDS_FILE_TOKEN_64 : FLUIDS_FILE_TOKEN_32;
231 PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT32));
232 PetscCall(PetscViewerBinaryWrite(viewer, &step_no, 1, PETSC_INT));
233 time /= user->units->second; // Dimensionalize time back
234 PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL));
235 PetscCall(VecView(Q, viewer));
236 PetscCall(PetscViewerDestroy(&viewer));
237 PetscFunctionReturn(PETSC_SUCCESS);
238 }
239
240 // CSV Monitor
TSMonitor_WallForce(TS ts,PetscInt step_no,PetscReal time,Vec Q,void * ctx)241 PetscErrorCode TSMonitor_WallForce(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
242 User user = ctx;
243 Vec G_loc;
244 PetscInt num_wall = user->app_ctx->wall_forces.num_wall, dim = 3;
245 const PetscInt *walls = user->app_ctx->wall_forces.walls;
246 PetscViewer viewer = user->app_ctx->wall_forces.viewer;
247 PetscViewerFormat format = user->app_ctx->wall_forces.viewer_format;
248 PetscScalar *reaction_force;
249 PetscBool iascii;
250
251 PetscFunctionBeginUser;
252 if (!viewer) PetscFunctionReturn(PETSC_SUCCESS);
253 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
254 PetscCall(PetscMalloc1(num_wall * dim, &reaction_force));
255 PetscCall(Surface_Forces_NS(user->dm, G_loc, num_wall, walls, reaction_force));
256 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
257
258 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
259
260 if (iascii) {
261 if (format == PETSC_VIEWER_ASCII_CSV && !user->app_ctx->wall_forces.header_written) {
262 PetscCall(PetscViewerASCIIPrintf(viewer, "Step,Time,Wall,ForceX,ForceY,ForceZ\n"));
263 user->app_ctx->wall_forces.header_written = PETSC_TRUE;
264 }
265 for (PetscInt w = 0; w < num_wall; w++) {
266 PetscInt wall = walls[w];
267 if (format == PETSC_VIEWER_ASCII_CSV) {
268 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ",%g,%" PetscInt_FMT ",%g,%g,%g\n", step_no, time, wall,
269 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
270
271 } else {
272 PetscCall(PetscViewerASCIIPrintf(viewer, "Wall %" PetscInt_FMT " Forces: Force_x = %12g, Force_y = %12g, Force_z = %12g\n", wall,
273 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
274 }
275 }
276 }
277 PetscCall(PetscFree(reaction_force));
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
281 // User provided TS Monitor
TSMonitor_NS(TS ts,PetscInt step_no,PetscReal time,Vec Q,void * ctx)282 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
283 User user = ctx;
284
285 PetscFunctionBeginUser;
286 // Print every 'checkpoint_interval' steps
287 if (user->app_ctx->checkpoint_interval <= 0 || step_no % user->app_ctx->checkpoint_interval != 0 ||
288 (user->app_ctx->cont_steps == step_no && step_no != 0)) {
289 PetscFunctionReturn(PETSC_SUCCESS);
290 }
291
292 PetscCall(WriteOutput(user, Q, step_no, time));
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
296 // TS: Create, setup, and solve
TSSolve_NS(DM dm,User user,AppCtx app_ctx,Physics phys,ProblemData problem,Vec * Q,PetscScalar * f_time,TS * ts)297 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, ProblemData problem, Vec *Q, PetscScalar *f_time, TS *ts) {
298 MPI_Comm comm = user->comm;
299 TSAdapt adapt;
300 PetscScalar final_time;
301
302 PetscFunctionBeginUser;
303 PetscCall(TSCreate(comm, ts));
304 PetscCall(TSSetDM(*ts, dm));
305 PetscCall(TSSetApplicationContext(*ts, user));
306 if (phys->implicit) {
307 PetscCall(TSSetType(*ts, TSBDF));
308 if (user->op_ifunction) {
309 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user));
310 } else { // Implicit integrators can fall back to using an RHSFunction
311 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
312 }
313 if (user->mat_ijacobian) {
314 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user));
315 }
316 } else {
317 PetscCheck(user->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction");
318 PetscCall(TSSetType(*ts, TSRK));
319 PetscCall(TSRKSetType(*ts, TSRK5F));
320 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
321 }
322 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second));
323 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
324 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE));
325 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second));
326 PetscCall(TSGetAdapt(*ts, &adapt));
327 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second));
328 PetscCall(TSSetFromOptions(*ts));
329 if (user->mat_ijacobian) {
330 if (app_ctx->amat_type && !strcmp(app_ctx->amat_type, MATSHELL)) {
331 SNES snes;
332 KSP ksp;
333 Mat Pmat, Amat;
334
335 PetscCall(TSGetSNES(*ts, &snes));
336 PetscCall(SNESGetKSP(snes, &ksp));
337 PetscCall(CreateSolveOperatorsFromMatCeed(ksp, user->mat_ijacobian, PETSC_FALSE, &Amat, &Pmat));
338 PetscCall(TSSetIJacobian(*ts, user->mat_ijacobian, Pmat, NULL, NULL));
339 PetscCall(MatDestroy(&Amat));
340 PetscCall(MatDestroy(&Pmat));
341 }
342 }
343 user->time_bc_set = -1.0; // require all BCs be updated
344 if (app_ctx->cont_steps) { // continue from previous timestep data
345 PetscInt count;
346 PetscViewer viewer;
347
348 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time
349 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer));
350 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL));
351 PetscCall(PetscViewerDestroy(&viewer));
352 PetscCheck(app_ctx->cont_steps != -1, comm, PETSC_ERR_ARG_INCOMP,
353 "-continue step number not specified, but checkpoint file does not contain a step number (likely written by older code version)");
354 }
355 PetscCall(TSSetTime(*ts, app_ctx->cont_time * user->units->second));
356 PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps));
357 }
358 if (app_ctx->test_type == TESTTYPE_NONE) {
359 PetscCall(TSMonitorSet(*ts, TSMonitor_NS, user, NULL));
360 }
361 if (app_ctx->wall_forces.viewer) {
362 PetscCall(TSMonitorSet(*ts, TSMonitor_WallForce, user, NULL));
363 }
364 if (app_ctx->turb_spanstats_enable) {
365 PetscCall(TSMonitorSet(*ts, TSMonitor_TurbulenceStatistics, user, NULL));
366 CeedScalar previous_time = app_ctx->cont_time * user->units->second;
367 PetscCallCeed(user->ceed,
368 CeedOperatorSetContextDouble(user->spanstats.op_stats_collect_ctx->op, user->spanstats.previous_time_label, &previous_time));
369 }
370 if (app_ctx->diff_filter_monitor) PetscCall(TSMonitorSet(*ts, TSMonitor_DifferentialFilter, user, NULL));
371
372 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(PrintRunInfo(user, user->phys, problem, *ts));
373 // Solve
374 PetscReal start_time;
375 PetscInt start_step;
376 PetscCall(TSGetTime(*ts, &start_time));
377 PetscCall(TSGetStepNumber(*ts, &start_step));
378
379 PetscCall(PetscLogDefaultBegin()); // So we can use PetscLogStageGetPerfInfo without -log_view
380 PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve");
381 PetscCall(TSSetTime(*ts, start_time));
382 PetscCall(TSSetStepNumber(*ts, start_step));
383 if (PetscPreLoadingOn) {
384 // LCOV_EXCL_START
385 SNES snes;
386 Vec Q_preload;
387 PetscReal rtol;
388 PetscCall(VecDuplicate(*Q, &Q_preload));
389 PetscCall(VecCopy(*Q, Q_preload));
390 PetscCall(TSGetSNES(*ts, &snes));
391 PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL));
392 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
393 PetscCall(TSSetSolution(*ts, Q_preload));
394 PetscCall(TSStep(*ts));
395 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
396 PetscCall(VecDestroy(&Q_preload));
397 // LCOV_EXCL_STOP
398 } else {
399 PetscCall(PetscBarrier((PetscObject)*ts));
400 PetscCall(TSSolve(*ts, *Q));
401 }
402 PetscPreLoadEnd();
403
404 PetscCall(TSGetSolveTime(*ts, &final_time));
405 *f_time = final_time;
406
407 if (app_ctx->test_type == TESTTYPE_NONE) {
408 PetscInt step_no;
409 PetscCall(TSGetStepNumber(*ts, &step_no));
410 if (user->app_ctx->checkpoint_interval > 0 || user->app_ctx->checkpoint_interval == -1) {
411 PetscCall(WriteOutput(user, *Q, step_no, final_time));
412 }
413
414 PetscLogStage stage_id;
415 PetscEventPerfInfo stage_perf;
416
417 PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id));
418 PetscCall(PetscLogStageGetPerfInfo(stage_id, &stage_perf));
419 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_perf.time));
420 }
421 PetscFunctionReturn(PETSC_SUCCESS);
422 }
423