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 HONEE
6
7 #include <ceed.h>
8 #include <petscdmplex.h>
9 #include <petscts.h>
10
11 #include <differential_filter.h>
12 #include <navierstokes.h>
13 #include <smartsim.h>
14 #include <spanstats.h>
15 #include "../qfunctions/newtonian_state.h"
16
17 // @brief Insert Boundary values if it's a new time
UpdateBoundaryValues(Honee honee,Vec Q_loc,PetscReal t)18 PetscErrorCode UpdateBoundaryValues(Honee honee, Vec Q_loc, PetscReal t) {
19 PetscFunctionBeginUser;
20 if (honee->time_bc_set != t) {
21 PetscCall(DMPlexInsertBoundaryValues(honee->dm, PETSC_TRUE, Q_loc, t, NULL, NULL, NULL));
22 honee->time_bc_set = t;
23 }
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
27 // RHS (Explicit time-stepper) function setup
28 // This is the RHS of the ODE, given as u_t = G(t,u)
29 // 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)30 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) {
31 Honee honee = *(Honee *)user_data;
32 Ceed ceed = honee->ceed;
33 PetscScalar dt;
34 Vec Q_loc = honee->Q_loc, R;
35 PetscMemType q_mem_type;
36
37 PetscFunctionBeginUser;
38 // Update time dependent data
39 PetscCall(UpdateBoundaryValues(honee, Q_loc, t));
40 if (honee->phys->solution_time_label)
41 PetscCallCeed(ceed, CeedOperatorSetContextDouble(honee->op_rhs_ctx->op, honee->phys->solution_time_label, &t));
42 PetscCall(TSGetTimeStep(ts, &dt));
43 if (honee->phys->timestep_size_label)
44 PetscCallCeed(ceed, CeedOperatorSetContextDouble(honee->op_rhs_ctx->op, honee->phys->timestep_size_label, &dt));
45
46 PetscCall(DMGetNamedGlobalVector(honee->dm, "RHS Residual", &R));
47 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, Q_loc));
48 if (honee->app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) PetscCall(DivDiffFluxProjectionApply(honee->diff_flux_proj, Q_loc));
49 PetscCall(ApplyCeedOperatorLocalToGlobal(Q_loc, R, honee->op_rhs_ctx));
50
51 // Inverse of the mass matrix
52 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, honee->q_ceed));
53 {
54 // Run PCApply manually if using ksp_type preonly -pc_type jacobi
55 // This is to avoid an AllReduce call in KSPSolve_Preonly, which causes significant slowdowns for lumped mass matrix solves.
56 // See https://gitlab.com/petsc/petsc/-/merge_requests/8048 for more details and a possible fix
57 PC pc;
58 PetscBool ispreonly, isjacobi;
59 PetscCall(KSPGetPC(honee->mass_ksp, &pc));
60 PetscCall(PetscObjectTypeCompare((PetscObject)honee->mass_ksp, KSPPREONLY, &ispreonly));
61 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &isjacobi));
62 if (ispreonly && isjacobi) PetscCall(PCApply(pc, R, G));
63 else PetscCall(KSPSolve(honee->mass_ksp, R, G));
64 }
65 PetscCall(VecReadCeedToPetsc(honee->q_ceed, q_mem_type, Q_loc));
66
67 PetscCall(DMRestoreNamedGlobalVector(honee->dm, "RHS Residual", &R));
68 PetscFunctionReturn(PETSC_SUCCESS);
69 }
70
71 // Surface forces function setup
Surface_Forces_NS(DM dm,Vec G_loc,PetscInt num_walls,const PetscInt walls[],PetscScalar * reaction_force)72 static PetscErrorCode Surface_Forces_NS(DM dm, Vec G_loc, PetscInt num_walls, const PetscInt walls[], PetscScalar *reaction_force) {
73 DMLabel face_label;
74 const PetscScalar *g_array;
75 PetscInt dim = 3;
76 MPI_Comm comm = PetscObjectComm((PetscObject)dm);
77 PetscSection section;
78
79 PetscFunctionBeginUser;
80 PetscCall(DMGetLabel(dm, "Face Sets", &face_label));
81 PetscCall(VecGetArrayRead(G_loc, &g_array));
82 for (PetscInt w = 0; w < num_walls; w++) {
83 const PetscInt wall = walls[w], *points;
84 IS wall_is;
85 PetscInt num_points, num_comp = 0;
86
87 PetscCall(DMLabelGetStratumIS(face_label, wall, &wall_is));
88 if (!wall_is) continue; // No wall points on this process, skip
89
90 PetscCall(DMGetLocalSection(dm, §ion));
91 PetscCall(PetscSectionGetFieldComponents(section, 0, &num_comp));
92 PetscCall(ISGetSize(wall_is, &num_points));
93 PetscCall(ISGetIndices(wall_is, &points));
94 for (PetscInt i = 0; i < num_points; i++) {
95 const PetscInt p = points[i];
96 const StateConservative *r;
97 PetscInt dof;
98
99 PetscCall(DMPlexPointLocalRead(dm, p, g_array, &r));
100 PetscCall(PetscSectionGetDof(section, p, &dof));
101 for (PetscInt node = 0; node < dof / num_comp; node++) {
102 for (PetscInt j = 0; j < dim; j++) {
103 reaction_force[w * dim + j] -= r[node].momentum[j];
104 }
105 }
106 }
107 PetscCall(ISRestoreIndices(wall_is, &points));
108 PetscCall(ISDestroy(&wall_is));
109 }
110 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, reaction_force, dim * num_walls, MPIU_SCALAR, MPI_SUM, comm));
111 PetscCall(VecRestoreArrayRead(G_loc, &g_array));
112 PetscFunctionReturn(PETSC_SUCCESS);
113 }
114
115 // Implicit time-stepper function setup
IFunction_NS(TS ts,PetscReal t,Vec Q,Vec Q_dot,Vec G,void * user_data)116 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, void *user_data) {
117 Honee honee = *(Honee *)user_data;
118 Ceed ceed = honee->ceed;
119 PetscScalar dt;
120 Vec Q_loc = honee->Q_loc, Q_dot_loc = honee->Q_dot_loc, G_loc;
121 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type;
122
123 PetscFunctionBeginUser;
124 PetscCall(DMGlobalToLocalBegin(honee->dm, Q_dot, INSERT_VALUES, Q_dot_loc));
125 PetscCall(DMGetNamedLocalVector(honee->dm, "ResidualLocal", &G_loc));
126
127 // Update time dependent data
128 PetscCall(UpdateBoundaryValues(honee, Q_loc, t));
129 if (honee->phys->solution_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(honee->op_ifunction, honee->phys->solution_time_label, &t));
130 PetscCall(TSGetTimeStep(ts, &dt));
131 if (honee->phys->timestep_size_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(honee->op_ifunction, honee->phys->timestep_size_label, &dt));
132
133 // Global-to-local
134 PetscCall(DMGlobalToLocalBegin(honee->dm, Q, INSERT_VALUES, Q_loc));
135 PetscCall(DMGlobalToLocalEnd(honee->dm, Q, INSERT_VALUES, Q_loc));
136 if (honee->app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) PetscCall(DivDiffFluxProjectionApply(honee->diff_flux_proj, Q_loc));
137 PetscCall(DMGlobalToLocalEnd(honee->dm, Q_dot, INSERT_VALUES, Q_dot_loc));
138
139 // Place PETSc vectors in CEED vectors
140 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, honee->q_ceed));
141 PetscCall(VecReadPetscToCeed(Q_dot_loc, &q_dot_mem_type, honee->q_dot_ceed));
142 PetscCall(VecPetscToCeed(G_loc, &g_mem_type, honee->g_ceed));
143
144 // Apply CEED operator
145 PetscCall(PetscLogEventBegin(HONEE_CeedOperatorApply, Q, G, 0, 0));
146 PetscCall(PetscLogGpuTimeBegin());
147 PetscCallCeed(honee->ceed, CeedOperatorApply(honee->op_ifunction, honee->q_ceed, honee->g_ceed, CEED_REQUEST_IMMEDIATE));
148 PetscCall(PetscLogGpuTimeEnd());
149 PetscCall(PetscLogEventEnd(HONEE_CeedOperatorApply, Q, G, 0, 0));
150
151 // Restore vectors
152 PetscCall(VecReadCeedToPetsc(honee->q_ceed, q_mem_type, Q_loc));
153 PetscCall(VecReadCeedToPetsc(honee->q_dot_ceed, q_dot_mem_type, Q_dot_loc));
154 PetscCall(VecCeedToPetsc(honee->g_ceed, g_mem_type, G_loc));
155
156 if (honee->app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) {
157 PetscCall(SgsDDApplyIFunction(honee, Q_loc, G_loc));
158 }
159
160 // Local-to-Global
161 PetscCall(VecZeroEntries(G));
162 PetscCall(DMLocalToGlobal(honee->dm, G_loc, ADD_VALUES, G));
163
164 // Restore vectors
165 PetscCall(DMRestoreNamedLocalVector(honee->dm, "ResidualLocal", &G_loc));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
FormIJacobian_NS(TS ts,PetscReal t,Vec Q,Vec Q_dot,PetscReal shift,Mat J,Mat J_pre,void * user_data)169 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) {
170 Honee honee = *(Honee *)user_data;
171 PetscBool J_is_matceed, J_is_mffd, J_pre_is_matceed, J_pre_is_mffd;
172
173 PetscFunctionBeginUser;
174 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd));
175 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATCEED, &J_is_matceed));
176 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATMFFD, &J_pre_is_mffd));
177 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATCEED, &J_pre_is_matceed));
178
179 PetscCall(MatCeedSetContextReal(honee->mat_ijacobian, "ijacobian time shift", shift));
180
181 if (J_is_matceed || J_is_mffd) {
182 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
183 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
184 } else PetscCall(MatCeedAssembleCOO(honee->mat_ijacobian, J));
185
186 if (J_pre_is_matceed && J != J_pre) {
187 PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY));
188 PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY));
189 } else if (!J_pre_is_matceed && !J_pre_is_mffd && J != J_pre) {
190 PetscCall(MatCeedAssembleCOO(honee->mat_ijacobian, J_pre));
191 }
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
WriteOutput(Honee honee,Vec Q,PetscInt step_no,PetscScalar time)195 PetscErrorCode WriteOutput(Honee honee, Vec Q, PetscInt step_no, PetscScalar time) {
196 Vec Q_loc;
197 char file_path[PETSC_MAX_PATH_LEN];
198 PetscViewer viewer;
199
200 PetscFunctionBeginUser;
201 if (honee->app_ctx->checkpoint_vtk) {
202 // Set up output
203 PetscCall(DMGetLocalVector(honee->dm, &Q_loc));
204 PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec"));
205 PetscCall(VecZeroEntries(Q_loc));
206 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, Q_loc));
207
208 // Output
209 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", honee->app_ctx->output_dir, step_no));
210
211 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer));
212 PetscCall(VecView(Q_loc, viewer));
213 PetscCall(PetscViewerDestroy(&viewer));
214 if (honee->dm_viz) {
215 Vec Q_refined, Q_refined_loc;
216 char file_path_refined[PETSC_MAX_PATH_LEN];
217 PetscViewer viewer_refined;
218
219 PetscCall(DMGetGlobalVector(honee->dm_viz, &Q_refined));
220 PetscCall(DMGetLocalVector(honee->dm_viz, &Q_refined_loc));
221 PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"));
222
223 PetscCall(MatInterpolate(honee->interp_viz, Q, Q_refined));
224 PetscCall(VecZeroEntries(Q_refined_loc));
225 PetscCall(DMGlobalToLocal(honee->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc));
226
227 PetscCall(PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", honee->app_ctx->output_dir,
228 step_no));
229
230 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined));
231 PetscCall(VecView(Q_refined_loc, viewer_refined));
232 PetscCall(DMRestoreLocalVector(honee->dm_viz, &Q_refined_loc));
233 PetscCall(DMRestoreGlobalVector(honee->dm_viz, &Q_refined));
234 PetscCall(PetscViewerDestroy(&viewer_refined));
235 }
236 PetscCall(DMRestoreLocalVector(honee->dm, &Q_loc));
237 }
238
239 // Save data in a binary file for continuation of simulations
240 if (honee->app_ctx->add_stepnum2bin) {
241 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", honee->app_ctx->output_dir, step_no));
242 } else {
243 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", honee->app_ctx->output_dir));
244 }
245 PetscCall(PetscViewerBinaryOpen(honee->comm, file_path, FILE_MODE_WRITE, &viewer));
246
247 time /= honee->units->second; // Dimensionalize time back
248 PetscCall(HoneeWriteBinaryVec(viewer, Q, time, step_no));
249 PetscCall(PetscViewerDestroy(&viewer));
250 PetscFunctionReturn(PETSC_SUCCESS);
251 }
252
253 // CSV Monitor
TSMonitor_WallForce(TS ts,PetscInt step_no,PetscReal time,Vec Q,void * ctx)254 PetscErrorCode TSMonitor_WallForce(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
255 Honee honee = ctx;
256 Vec G_loc;
257 PetscInt num_wall = honee->app_ctx->wall_forces.num_wall, dim = 3;
258 const PetscInt *walls = honee->app_ctx->wall_forces.walls;
259 PetscViewer viewer = honee->app_ctx->wall_forces.viewer;
260 PetscViewerFormat format = honee->app_ctx->wall_forces.viewer_format;
261 PetscScalar *reaction_force;
262 PetscBool is_ascii;
263
264 PetscFunctionBeginUser;
265 if (!viewer) PetscFunctionReturn(PETSC_SUCCESS);
266 PetscCall(DMGetNamedLocalVector(honee->dm, "ResidualLocal", &G_loc));
267 PetscCall(PetscCalloc1(num_wall * dim, &reaction_force));
268 PetscCall(Surface_Forces_NS(honee->dm, G_loc, num_wall, walls, reaction_force));
269 PetscCall(DMRestoreNamedLocalVector(honee->dm, "ResidualLocal", &G_loc));
270
271 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &is_ascii));
272
273 if (is_ascii) {
274 if (format == PETSC_VIEWER_ASCII_CSV && !honee->app_ctx->wall_forces.header_written) {
275 PetscCall(PetscViewerASCIIPrintf(viewer, "Step,Time,Wall,ForceX,ForceY,ForceZ\n"));
276 honee->app_ctx->wall_forces.header_written = PETSC_TRUE;
277 }
278 for (PetscInt w = 0; w < num_wall; w++) {
279 PetscInt wall = walls[w];
280 if (format == PETSC_VIEWER_ASCII_CSV) {
281 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ",%g,%" PetscInt_FMT ",%g,%g,%g\n", step_no, time, wall,
282 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
283
284 } else {
285 PetscCall(PetscViewerASCIIPrintf(viewer, "Wall %" PetscInt_FMT " Forces: Force_x = %12g, Force_y = %12g, Force_z = %12g\n", wall,
286 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2]));
287 }
288 }
289 }
290 PetscCall(PetscFree(reaction_force));
291 PetscFunctionReturn(PETSC_SUCCESS);
292 }
293
294 // User provided TS Monitor
TSMonitor_NS(TS ts,PetscInt step_no,PetscReal time,Vec Q,void * ctx)295 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) {
296 Honee honee = ctx;
297
298 PetscFunctionBeginUser;
299 // Print every 'checkpoint_interval' steps
300 if (honee->app_ctx->checkpoint_interval <= 0 || step_no % honee->app_ctx->checkpoint_interval != 0 ||
301 (honee->app_ctx->cont_steps == step_no && step_no != 0)) {
302 PetscFunctionReturn(PETSC_SUCCESS);
303 }
304
305 PetscCall(WriteOutput(honee, Q, step_no, time));
306 PetscFunctionReturn(PETSC_SUCCESS);
307 }
308
TSPostStep_CheckStep(TS ts)309 PetscErrorCode TSPostStep_CheckStep(TS ts) {
310 Honee honee;
311 PetscReal norm;
312 PetscInt step;
313 Vec Q;
314
315 PetscFunctionBeginUser;
316 PetscCall(TSGetApplicationContext(ts, &honee));
317 PetscCall(TSGetStepNumber(ts, &step));
318 PetscCall(TSGetSolution(ts, &Q));
319 if (step % honee->app_ctx->check_step_interval) PetscFunctionReturn(PETSC_SUCCESS);
320 PetscCall(VecNorm(Q, NORM_1, &norm));
321 if (PetscIsInfOrNanReal(norm)) {
322 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Solution diverged: Nans found in solution\n"));
323 PetscCall(TSSetConvergedReason(ts, TS_DIVERGED_NONLINEAR_SOLVE));
324 }
325 PetscFunctionReturn(PETSC_SUCCESS);
326 }
327
TSPostStep_MaxWallTime(TS ts)328 PetscErrorCode TSPostStep_MaxWallTime(TS ts) {
329 Honee honee;
330 PetscInt step;
331 PetscMPIInt rank;
332 MPI_Comm comm;
333 PetscBool is_wall_time_exceeded = PETSC_FALSE;
334
335 PetscFunctionBeginUser;
336 PetscCall(TSGetApplicationContext(ts, &honee));
337 PetscCall(TSGetStepNumber(ts, &step));
338 if (step % honee->max_wall_time_interval) PetscFunctionReturn(PETSC_SUCCESS);
339 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
340 PetscCallMPI(MPI_Comm_rank(comm, &rank));
341 if (rank == 0) is_wall_time_exceeded = time(NULL) > honee->max_wall_time ? PETSC_TRUE : PETSC_FALSE;
342 // Must broadcast to avoid race condition
343 PetscCallMPI(MPI_Bcast(&is_wall_time_exceeded, 1, MPI_C_BOOL, 0, comm));
344 if (PetscUnlikely(is_wall_time_exceeded)) {
345 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Stopping TSSolve: Set max wall time exceeded\n"));
346 PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER));
347 }
348 PetscFunctionReturn(PETSC_SUCCESS);
349 }
350
351 /**
352 @brief TSPostStep for HONEE
353
354 `TSSetPostStep()` only accepts a single function argument, so this function groups together all post-step
355 functionality needed for HONEE features
356
357 @param[in] ts TS object
358 **/
TSPostStep_Honee(TS ts)359 PetscErrorCode TSPostStep_Honee(TS ts) {
360 Honee honee;
361
362 PetscFunctionBeginUser;
363 PetscCall(TSGetApplicationContext(ts, &honee));
364 if (honee->max_wall_time > 0) PetscCall(TSPostStep_MaxWallTime(ts));
365 if (honee->app_ctx->sgs_train_enable) PetscCall(TSPostStep_SGS_DD_Training(ts));
366 if (honee->app_ctx->check_step_interval > 0) PetscCall(TSPostStep_CheckStep(ts));
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
370 /**
371 @brief Save `TSEvaluationSolutions()` to file defined by `-ts_eval_solutions_view`
372
373 Intended to be used with `TSSetEvaluationTimes()`/`-ts_eval_times`
374
375 @param[in] ts TS object that has the evaluation solutions
376 **/
HoneeTSEvaluationSolutions(TS ts)377 static PetscErrorCode HoneeTSEvaluationSolutions(TS ts) {
378 MPI_Comm comm = PetscObjectComm((PetscObject)ts);
379 const PetscReal *sol_times;
380 PetscReal orig_time;
381 PetscInt num_sols, orig_step;
382 Vec *sol_vecs;
383 DM dm;
384 PetscBool is_viewer_set;
385 PetscViewer viewer;
386 PetscViewerFormat format;
387
388 PetscFunctionBeginUser;
389 PetscCall(TSGetEvaluationSolutions(ts, &num_sols, &sol_times, &sol_vecs));
390 if (num_sols == 0) PetscFunctionReturn(PETSC_SUCCESS);
391 PetscCall(TSGetDM(ts, &dm));
392 PetscCall(DMGetOutputSequenceNumber(dm, &orig_step, &orig_time));
393
394 const char *option = "-ts_eval_solutions_view";
395 PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, option, &viewer, &format, &is_viewer_set));
396 if (!is_viewer_set) {
397 PetscCall(PetscPrintf(comm, "\n\nWARNING: Viewer not set for TSEvaluationSolutions. Set %s to save this data. Now throwing it away!\n", option));
398 PetscFunctionReturn(PETSC_SUCCESS);
399 }
400 for (PetscInt i = 0; i < num_sols; i++) {
401 PetscCall(DMSetOutputSequenceNumber(dm, i, sol_times[i]));
402 PetscCall(VecView(sol_vecs[i], viewer));
403 }
404 PetscCall(PetscViewerPopFormat(viewer));
405 PetscCall(PetscViewerDestroy(&viewer));
406 PetscFunctionReturn(PETSC_SUCCESS);
407 }
408
409 // TS: Create, setup, and solve
TSSolve_NS(DM dm,Honee honee,AppCtx app_ctx,Physics phys,ProblemData problem,Vec Q,PetscScalar * f_time,TS * ts)410 PetscErrorCode TSSolve_NS(DM dm, Honee honee, AppCtx app_ctx, Physics phys, ProblemData problem, Vec Q, PetscScalar *f_time, TS *ts) {
411 MPI_Comm comm = honee->comm;
412 TSAdapt adapt;
413 PetscScalar final_time;
414
415 PetscFunctionBeginUser;
416 PetscCall(TSCreate(comm, ts));
417 PetscCall(TSSetDM(*ts, dm));
418 PetscCall(TSSetApplicationContext(*ts, honee));
419 if (phys->implicit) {
420 PetscCall(TSSetType(*ts, TSBDF));
421 if (honee->op_ifunction) {
422 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &honee));
423 } else { // Implicit integrators can fall back to using an RHSFunction
424 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &honee));
425 }
426 if (honee->mat_ijacobian) {
427 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &honee));
428 }
429 } else {
430 PetscCheck(honee->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction");
431 PetscCall(TSSetType(*ts, TSRK));
432 PetscCall(TSRKSetType(*ts, TSRK5F));
433 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &honee));
434 }
435 PetscCall(TSSetMaxTime(*ts, 500. * honee->units->second));
436 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
437 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE));
438 PetscCall(TSSetTimeStep(*ts, 1.e-2 * honee->units->second));
439 PetscCall(TSGetAdapt(*ts, &adapt));
440 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * honee->units->second, 1.e2 * honee->units->second));
441 PetscCall(TSSetFromOptions(*ts));
442 if (honee->mat_ijacobian) {
443 if (app_ctx->amat_type && !strcmp(app_ctx->amat_type, MATSHELL)) {
444 SNES snes;
445 KSP ksp;
446 Mat Pmat, Amat;
447
448 PetscCall(TSGetSNES(*ts, &snes));
449 PetscCall(SNESGetKSP(snes, &ksp));
450 PetscCall(CreateSolveOperatorsFromMatCeed(ksp, honee->mat_ijacobian, PETSC_FALSE, &Amat, &Pmat));
451 PetscCall(TSSetIJacobian(*ts, honee->mat_ijacobian, Pmat, NULL, NULL));
452 PetscCall(MatDestroy(&Amat));
453 PetscCall(MatDestroy(&Pmat));
454 }
455 }
456 honee->time_bc_set = -1.0; // require all BCs be updated
457 if (app_ctx->cont_steps) { // continue from previous timestep data
458 PetscCall(TSSetTime(*ts, app_ctx->cont_time * honee->units->second));
459 PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps));
460 }
461
462 PetscBool add_ksp_postsolve_residual = PETSC_FALSE;
463 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_post_solve_residual", &add_ksp_postsolve_residual, NULL));
464 if (add_ksp_postsolve_residual) {
465 SNES snes;
466 KSP ksp;
467
468 PetscCall(TSGetSNES(*ts, &snes));
469 PetscCall(SNESGetKSP(snes, &ksp));
470 PetscCall(KSPSetResidualHistory(ksp, NULL, PETSC_DECIDE, PETSC_TRUE));
471 PetscCall(KSPSetPostSolve(ksp, KSPPostSolve_Honee, NULL));
472 }
473 if (honee->set_poststep) PetscCall(TSSetPostStep(*ts, TSPostStep_Honee));
474 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSMonitorSet(*ts, TSMonitor_NS, honee, NULL));
475 if (app_ctx->wall_forces.viewer) PetscCall(TSMonitorSet(*ts, TSMonitor_WallForce, honee, NULL));
476
477 PetscOptionsBegin(comm, NULL, "HONEE TS Monitor Options", NULL);
478 PetscCall(TSMonitorSetFromOptions(*ts, "-ts_monitor_total_kinetic_energy", "Monitor total kinetic energy balance terms in the domain", NULL,
479 TSMonitor_TotalKineticEnergy, SetupMontiorTotalKineticEnergy));
480 PetscCall(TSMonitorSetFromOptions(*ts, "-ts_monitor_cfl", "Monitor element CFL statistics", NULL, TSMonitor_Cfl, SetupMontiorCfl));
481
482 PetscCall(PetscOptionsDeprecated("-ts_monitor_turbulence_spanstats_viewer_interval", "-ts_monitor_spanstats_turbulence_interval", "HONEE 0.0",
483 NULL));
484 PetscCall(PetscOptionsDeprecated("-ts_monitor_turbulence_spanstats_viewer", "-ts_monitor_spanstats_turbulence", "HONEE 0.0", NULL));
485 PetscCall(TSMonitorSetFromOptions(*ts, "-ts_monitor_spanstats_turbulence", "Setup spanwise statistics collection", NULL,
486 TSMonitor_SpanwiseStatisticsTurbulence, SpanwiseStatisticsSetup_Turbulence));
487 PetscCall(TSMonitorSetFromOptions(*ts, "-ts_monitor_spanstats_cflpe", "Setup spanwise statistics collection of CFL and Pe", NULL,
488 TSMonitor_SpanwiseStatisticsCflPe, SpanwiseStatisticsSetup_CflPe));
489 PetscCall(TSMonitorSetFromOptions(*ts, "-diff_filter_monitor", "Run differential filtering for every timestep (used for testing)", NULL,
490 TSMonitor_DifferentialFilter, TSMonitor_DifferentialFilterSetup));
491 PetscCall(TSMonitorSetFromOptions(*ts, "-ts_monitor_smartsim_solution", "Write solution to SmartSim database", NULL, TSMonitor_SmartSimSolution,
492 TSMonitor_SmartSimSolutionSetup));
493 PetscOptionsEnd();
494
495 if (app_ctx->sgs_train_enable) PetscCall(TSMonitorSet(*ts, TSMonitor_SGS_DD_Training, honee, NULL));
496
497 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(PrintRunInfo(honee, honee->phys, problem, *ts));
498 if (honee->mass_ksp) PetscCall(KSPViewFromOptions(honee->mass_ksp, NULL, "-ksp_view_pre_ts_solve"));
499
500 // Solve
501 PetscReal start_time;
502 PetscInt start_step;
503 PetscCall(TSGetTime(*ts, &start_time));
504 PetscCall(TSGetStepNumber(*ts, &start_step));
505
506 PetscCall(PetscLogDefaultBegin()); // So we can use PetscLogStageGetPerfInfo without -log_view
507 PetscPreLoadBegin(PETSC_FALSE, "HONEE Solve");
508 PetscCall(TSSetTime(*ts, start_time));
509 PetscCall(TSSetStepNumber(*ts, start_step));
510 if (PetscPreLoadingOn) {
511 // LCOV_EXCL_START
512 SNES snes;
513 KSP ksp;
514 Vec Q_preload;
515 PetscReal rtol_snes, rtol_ksp;
516 PetscCall(VecDuplicate(Q, &Q_preload));
517 PetscCall(VecCopy(Q, Q_preload));
518 PetscCall(TSGetSNES(*ts, &snes));
519 PetscCall(SNESGetTolerances(snes, NULL, &rtol_snes, NULL, NULL, NULL));
520 PetscCall(SNESSetTolerances(snes, PETSC_CURRENT, .99, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
521 PetscCall(SNESGetKSP(snes, &ksp));
522 PetscCall(KSPGetTolerances(ksp, &rtol_ksp, NULL, NULL, NULL));
523 PetscCall(KSPSetTolerances(ksp, .99, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
524 PetscCall(TSSetSolution(*ts, Q_preload));
525 PetscCall(TSStep(*ts));
526 PetscCall(SNESSetTolerances(snes, PETSC_CURRENT, rtol_snes, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
527 PetscCall(KSPSetTolerances(ksp, rtol_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
528 PetscCall(VecDestroy(&Q_preload));
529 // LCOV_EXCL_STOP
530 } else {
531 PetscCall(PetscBarrier((PetscObject)*ts));
532 PetscCall(TSSolve(*ts, Q));
533 }
534 PetscPreLoadEnd();
535
536 PetscCall(TSGetSolveTime(*ts, &final_time));
537 *f_time = final_time;
538
539 if (app_ctx->test_type == TESTTYPE_NONE) {
540 PetscInt step_no;
541 PetscLogStage stage_id;
542 PetscEventPerfInfo stage_perf;
543
544 PetscCall(TSGetStepNumber(*ts, &step_no));
545 if (honee->app_ctx->checkpoint_interval > 0 || honee->app_ctx->checkpoint_interval == -1) {
546 PetscCall(WriteOutput(honee, Q, step_no, final_time));
547 }
548 PetscCall(HoneeTSEvaluationSolutions(*ts));
549
550 PetscCall(PetscLogStageGetId("HONEE Solve", &stage_id));
551 PetscCall(PetscLogStageGetPerfInfo(stage_id, &stage_perf));
552 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_perf.time));
553 }
554 PetscFunctionReturn(PETSC_SUCCESS);
555 }
556