xref: /honee/src/setupts.c (revision 51bb547f138f3ee4bfebc6254d262ff074a6b0ca)
1 // Copyright (c) 2017-2024, 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
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
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
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
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   if (user->app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) {
144     PetscCall(SgsDDApplyIFunction(user, Q_loc, G_loc));
145   }
146 
147   // Local-to-Global
148   PetscCall(VecZeroEntries(G));
149   PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G));
150 
151   // Restore vectors
152   PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
153   PetscFunctionReturn(PETSC_SUCCESS);
154 }
155 
156 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) {
157   User      user = *(User *)user_data;
158   PetscBool J_is_matceed, J_is_mffd, J_pre_is_matceed, J_pre_is_mffd;
159 
160   PetscFunctionBeginUser;
161   PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd));
162   PetscCall(PetscObjectTypeCompare((PetscObject)J, MATCEED, &J_is_matceed));
163   PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATMFFD, &J_pre_is_mffd));
164   PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATCEED, &J_pre_is_matceed));
165 
166   PetscCall(MatCeedSetContextReal(user->mat_ijacobian, "ijacobian time shift", shift));
167 
168   if (J_is_matceed || J_is_mffd) {
169     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
170     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
171   } else PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J));
172 
173   if (J_pre_is_matceed && J != J_pre) {
174     PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY));
175     PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY));
176   } else if (!J_pre_is_matceed && !J_pre_is_mffd && J != J_pre) {
177     PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J_pre));
178   }
179   PetscFunctionReturn(PETSC_SUCCESS);
180 }
181 
182 PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) {
183   Vec         Q_loc;
184   char        file_path[PETSC_MAX_PATH_LEN];
185   PetscViewer viewer;
186 
187   PetscFunctionBeginUser;
188   if (user->app_ctx->checkpoint_vtk) {
189     // Set up output
190     PetscCall(DMGetLocalVector(user->dm, &Q_loc));
191     PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec"));
192     PetscCall(VecZeroEntries(Q_loc));
193     PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc));
194 
195     // Output
196     PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no));
197 
198     PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer));
199     PetscCall(VecView(Q_loc, viewer));
200     PetscCall(PetscViewerDestroy(&viewer));
201     if (user->dm_viz) {
202       Vec         Q_refined, Q_refined_loc;
203       char        file_path_refined[PETSC_MAX_PATH_LEN];
204       PetscViewer viewer_refined;
205 
206       PetscCall(DMGetGlobalVector(user->dm_viz, &Q_refined));
207       PetscCall(DMGetLocalVector(user->dm_viz, &Q_refined_loc));
208       PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"));
209 
210       PetscCall(MatInterpolate(user->interp_viz, Q, Q_refined));
211       PetscCall(VecZeroEntries(Q_refined_loc));
212       PetscCall(DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc));
213 
214       PetscCall(
215           PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no));
216 
217       PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined));
218       PetscCall(VecView(Q_refined_loc, viewer_refined));
219       PetscCall(DMRestoreLocalVector(user->dm_viz, &Q_refined_loc));
220       PetscCall(DMRestoreGlobalVector(user->dm_viz, &Q_refined));
221       PetscCall(PetscViewerDestroy(&viewer_refined));
222     }
223     PetscCall(DMRestoreLocalVector(user->dm, &Q_loc));
224   }
225 
226   // Save data in a binary file for continuation of simulations
227   if (user->app_ctx->add_stepnum2bin) {
228     PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", user->app_ctx->output_dir, step_no));
229   } else {
230     PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", user->app_ctx->output_dir));
231   }
232   PetscCall(PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer));
233 
234   PetscInt32 token = PetscDefined(USE_64BIT_INDICES) ? FLUIDS_FILE_TOKEN_64 : FLUIDS_FILE_TOKEN_32;
235   PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT32));
236   PetscCall(PetscViewerBinaryWrite(viewer, &step_no, 1, PETSC_INT));
237   time /= user->units->second;  // Dimensionalize time back
238   PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL));
239   PetscCall(VecView(Q, viewer));
240   PetscCall(PetscViewerDestroy(&viewer));
241   PetscFunctionReturn(PETSC_SUCCESS);
242 }
243 
244 // CSV Monitor
245 PetscErrorCode TSMonitor_WallForce(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
246   User              user = ctx;
247   Vec               G_loc;
248   PetscInt          num_wall = user->app_ctx->wall_forces.num_wall, dim = 3;
249   const PetscInt   *walls  = user->app_ctx->wall_forces.walls;
250   PetscViewer       viewer = user->app_ctx->wall_forces.viewer;
251   PetscViewerFormat format = user->app_ctx->wall_forces.viewer_format;
252   PetscScalar      *reaction_force;
253   PetscBool         iascii;
254 
255   PetscFunctionBeginUser;
256   if (!viewer) PetscFunctionReturn(PETSC_SUCCESS);
257   PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
258   PetscCall(PetscMalloc1(num_wall * dim, &reaction_force));
259   PetscCall(Surface_Forces_NS(user->dm, G_loc, num_wall, walls, reaction_force));
260   PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
261 
262   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
263 
264   if (iascii) {
265     if (format == PETSC_VIEWER_ASCII_CSV && !user->app_ctx->wall_forces.header_written) {
266       PetscCall(PetscViewerASCIIPrintf(viewer, "Step,Time,Wall,ForceX,ForceY,ForceZ\n"));
267       user->app_ctx->wall_forces.header_written = PETSC_TRUE;
268     }
269     for (PetscInt w = 0; w < num_wall; w++) {
270       PetscInt wall = walls[w];
271       if (format == PETSC_VIEWER_ASCII_CSV) {
272         PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ",%g,%" PetscInt_FMT ",%g,%g,%g\n", step_no, time, wall,
273                                          reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
274 
275       } else {
276         PetscCall(PetscViewerASCIIPrintf(viewer, "Wall %" PetscInt_FMT " Forces: Force_x = %12g, Force_y = %12g, Force_z = %12g\n", wall,
277                                          reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
278       }
279     }
280   }
281   PetscCall(PetscFree(reaction_force));
282   PetscFunctionReturn(PETSC_SUCCESS);
283 }
284 
285 // User provided TS Monitor
286 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
287   User user = ctx;
288 
289   PetscFunctionBeginUser;
290   // Print every 'checkpoint_interval' steps
291   if (user->app_ctx->checkpoint_interval <= 0 || step_no % user->app_ctx->checkpoint_interval != 0 ||
292       (user->app_ctx->cont_steps == step_no && step_no != 0)) {
293     PetscFunctionReturn(PETSC_SUCCESS);
294   }
295 
296   PetscCall(WriteOutput(user, Q, step_no, time));
297   PetscFunctionReturn(PETSC_SUCCESS);
298 }
299 
300 // TS: Create, setup, and solve
301 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, Vec *Q, PetscScalar *f_time, TS *ts) {
302   MPI_Comm    comm = user->comm;
303   TSAdapt     adapt;
304   PetscScalar final_time;
305 
306   PetscFunctionBeginUser;
307   PetscCall(TSCreate(comm, ts));
308   PetscCall(TSSetDM(*ts, dm));
309   PetscCall(TSSetApplicationContext(*ts, user));
310   if (phys->implicit) {
311     PetscCall(TSSetType(*ts, TSBDF));
312     if (user->op_ifunction) {
313       PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user));
314     } else {  // Implicit integrators can fall back to using an RHSFunction
315       PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
316     }
317     if (user->mat_ijacobian) {
318       PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user));
319     }
320   } else {
321     PetscCheck(user->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction");
322     PetscCall(TSSetType(*ts, TSRK));
323     PetscCall(TSRKSetType(*ts, TSRK5F));
324     PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
325   }
326   PetscCall(TSSetMaxTime(*ts, 500. * user->units->second));
327   PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
328   if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE));
329   PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second));
330   PetscCall(TSGetAdapt(*ts, &adapt));
331   PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second));
332   PetscCall(TSSetFromOptions(*ts));
333   if (user->mat_ijacobian) {
334     if (app_ctx->amat_type && !strcmp(app_ctx->amat_type, MATSHELL)) {
335       SNES snes;
336       KSP  ksp;
337       Mat  Pmat, Amat;
338 
339       PetscCall(TSGetSNES(*ts, &snes));
340       PetscCall(SNESGetKSP(snes, &ksp));
341       PetscCall(CreateSolveOperatorsFromMatCeed(ksp, user->mat_ijacobian, PETSC_FALSE, &Amat, &Pmat));
342       PetscCall(TSSetIJacobian(*ts, user->mat_ijacobian, Pmat, NULL, NULL));
343       PetscCall(MatDestroy(&Amat));
344       PetscCall(MatDestroy(&Pmat));
345     }
346   }
347   user->time_bc_set = -1.0;   // require all BCs be updated
348   if (app_ctx->cont_steps) {  // continue from previous timestep data
349     PetscInt    count;
350     PetscViewer viewer;
351 
352     if (app_ctx->cont_time <= 0) {  // Legacy files did not include step number and time
353       PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer));
354       PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL));
355       PetscCall(PetscViewerDestroy(&viewer));
356       PetscCheck(app_ctx->cont_steps != -1, comm, PETSC_ERR_ARG_INCOMP,
357                  "-continue step number not specified, but checkpoint file does not contain a step number (likely written by older code version)");
358     }
359     PetscCall(TSSetTime(*ts, app_ctx->cont_time * user->units->second));
360     PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps));
361   }
362   if (app_ctx->test_type == TESTTYPE_NONE) {
363     PetscCall(TSMonitorSet(*ts, TSMonitor_NS, user, NULL));
364   }
365   if (app_ctx->wall_forces.viewer) {
366     PetscCall(TSMonitorSet(*ts, TSMonitor_WallForce, user, NULL));
367   }
368   if (app_ctx->turb_spanstats_enable) {
369     PetscCall(TSMonitorSet(*ts, TSMonitor_TurbulenceStatistics, user, NULL));
370     CeedScalar previous_time = app_ctx->cont_time * user->units->second;
371     PetscCallCeed(user->ceed,
372                   CeedOperatorSetContextDouble(user->spanstats.op_stats_collect_ctx->op, user->spanstats.previous_time_label, &previous_time));
373   }
374   if (app_ctx->diff_filter_monitor) PetscCall(TSMonitorSet(*ts, TSMonitor_DifferentialFilter, user, NULL));
375 
376   if (app_ctx->sgs_train_enable) {
377     PetscCall(TSMonitorSet(*ts, TSMonitor_SGS_DD_Training, user, NULL));
378     PetscCall(TSSetPostStep(*ts, TSPostStep_SGS_DD_Training));
379   }
380   // Solve
381   PetscReal start_time;
382   PetscInt  start_step;
383   PetscCall(TSGetTime(*ts, &start_time));
384   PetscCall(TSGetStepNumber(*ts, &start_step));
385 
386   PetscCall(PetscLogDefaultBegin());  // So we can use PetscLogStageGetPerfInfo without -log_view
387   PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve");
388   PetscCall(TSSetTime(*ts, start_time));
389   PetscCall(TSSetStepNumber(*ts, start_step));
390   if (PetscPreLoadingOn) {
391     // LCOV_EXCL_START
392     SNES      snes;
393     Vec       Q_preload;
394     PetscReal rtol;
395     PetscCall(VecDuplicate(*Q, &Q_preload));
396     PetscCall(VecCopy(*Q, Q_preload));
397     PetscCall(TSGetSNES(*ts, &snes));
398     PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL));
399     PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
400     PetscCall(TSSetSolution(*ts, Q_preload));
401     PetscCall(TSStep(*ts));
402     PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
403     PetscCall(VecDestroy(&Q_preload));
404     // LCOV_EXCL_STOP
405   } else {
406     PetscCall(PetscBarrier((PetscObject)*ts));
407     PetscCall(TSSolve(*ts, *Q));
408   }
409   PetscPreLoadEnd();
410 
411   PetscCall(TSGetSolveTime(*ts, &final_time));
412   *f_time = final_time;
413 
414   if (app_ctx->test_type == TESTTYPE_NONE) {
415     PetscInt step_no;
416     PetscCall(TSGetStepNumber(*ts, &step_no));
417     if (user->app_ctx->checkpoint_interval > 0 || user->app_ctx->checkpoint_interval == -1) {
418       PetscCall(WriteOutput(user, *Q, step_no, final_time));
419     }
420 
421     PetscLogStage      stage_id;
422     PetscEventPerfInfo stage_perf;
423 
424     PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id));
425     PetscCall(PetscLogStageGetPerfInfo(stage_id, &stage_perf));
426     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_perf.time));
427   }
428   PetscFunctionReturn(PETSC_SUCCESS);
429 }
430