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