xref: /libCEED/examples/fluids/src/setupts.c (revision 056ea4bdb268190ce972d0582063449a424e3e57) 
1 // Copyright (c) 2017-2022, 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 "../navierstokes.h"
12 #include "../qfunctions/mass.h"
13 
14 // Compute mass matrix for explicit scheme
15 PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, CeedData ceed_data,
16                                        Vec M) {
17   Vec            M_loc;
18   CeedQFunction  qf_mass;
19   CeedOperator   op_mass;
20   CeedVector     m_ceed, ones_vec;
21   CeedInt        num_comp_q, q_data_size;
22   PetscErrorCode ierr;
23   PetscFunctionBeginUser;
24 
25   // CEED Restriction
26   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q);
27   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size);
28   CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &m_ceed, NULL);
29   CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &ones_vec, NULL);
30   CeedVectorSetValue(ones_vec, 1.0);
31 
32   // CEED QFunction
33   CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass);
34   CeedQFunctionAddInput(qf_mass, "q", num_comp_q, CEED_EVAL_INTERP);
35   CeedQFunctionAddInput(qf_mass, "qdata", q_data_size, CEED_EVAL_NONE);
36   CeedQFunctionAddOutput(qf_mass, "v", num_comp_q, CEED_EVAL_INTERP);
37 
38   // CEED Operator
39   CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass);
40   CeedOperatorSetField(op_mass, "q", ceed_data->elem_restr_q, ceed_data->basis_q,
41                        CEED_VECTOR_ACTIVE);
42   CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i,
43                        CEED_BASIS_COLLOCATED, ceed_data->q_data);
44   CeedOperatorSetField(op_mass, "v", ceed_data->elem_restr_q, ceed_data->basis_q,
45                        CEED_VECTOR_ACTIVE);
46 
47   // Place PETSc vector in CEED vector
48   CeedScalar *m;
49   PetscMemType m_mem_type;
50   ierr = DMGetLocalVector(dm, &M_loc); CHKERRQ(ierr);
51   ierr = VecGetArrayAndMemType(M_loc, (PetscScalar **)&m, &m_mem_type);
52   CHKERRQ(ierr);
53   CeedVectorSetArray(m_ceed, MemTypeP2C(m_mem_type), CEED_USE_POINTER, m);
54 
55   // Apply CEED Operator
56   CeedOperatorApply(op_mass, ones_vec, m_ceed, CEED_REQUEST_IMMEDIATE);
57 
58   // Restore vectors
59   CeedVectorTakeArray(m_ceed, MemTypeP2C(m_mem_type), NULL);
60   ierr = VecRestoreArrayReadAndMemType(M_loc, (const PetscScalar **)&m);
61   CHKERRQ(ierr);
62 
63   // Local-to-Global
64   ierr = VecZeroEntries(M); CHKERRQ(ierr);
65   ierr = DMLocalToGlobal(dm, M_loc, ADD_VALUES, M); CHKERRQ(ierr);
66   ierr = DMRestoreLocalVector(dm, &M_loc); CHKERRQ(ierr);
67 
68   // Invert diagonally lumped mass vector for RHS function
69   ierr = VecReciprocal(M); CHKERRQ(ierr);
70 
71   // Cleanup
72   CeedVectorDestroy(&ones_vec);
73   CeedVectorDestroy(&m_ceed);
74   CeedQFunctionDestroy(&qf_mass);
75   CeedOperatorDestroy(&op_mass);
76 
77   PetscFunctionReturn(0);
78 }
79 
80 // RHS (Explicit time-stepper) function setup
81 //   This is the RHS of the ODE, given as u_t = G(t,u)
82 //   This function takes in a state vector Q and writes into G
83 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) {
84 
85   User           user = *(User *)user_data;
86   PetscScalar    *q, *g;
87   Vec            Q_loc, G_loc;
88   PetscMemType   q_mem_type, g_mem_type;
89   PetscErrorCode ierr;
90   PetscFunctionBeginUser;
91 
92   // Update context field labels
93   if (user->phys->solution_time_label)
94     CeedOperatorContextSetDouble(user->op_rhs, user->phys->solution_time_label, &t);
95   if (user->phys->timestep_size_label) {
96     PetscScalar dt;
97     ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
98     CeedOperatorContextSetDouble(user->op_rhs, user->phys->timestep_size_label,
99                                  &dt);
100   }
101 
102   // Get local vectors
103   ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
104   ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
105 
106   // Global-to-local
107   ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr);
108   ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr);
109   ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t,
110                                     NULL, NULL, NULL); CHKERRQ(ierr);
111   ierr = VecZeroEntries(G_loc); CHKERRQ(ierr);
112 
113   // Place PETSc vectors in CEED vectors
114   ierr = VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type);
115   CHKERRQ(ierr);
116   ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr);
117   CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, q);
118   CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g);
119 
120   // Apply CEED operator
121   CeedOperatorApply(user->op_rhs, user->q_ceed, user->g_ceed,
122                     CEED_REQUEST_IMMEDIATE);
123 
124   // Restore vectors
125   CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL);
126   CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL);
127   ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q);
128   CHKERRQ(ierr);
129   ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr);
130 
131   // Local-to-Global
132   ierr = VecZeroEntries(G); CHKERRQ(ierr);
133   ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr);
134 
135   // Inverse of the lumped mass matrix (M is Minv)
136   ierr = VecPointwiseMult(G, G, user->M); CHKERRQ(ierr);
137 
138   // Restore vectors
139   ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
140   ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
141 
142   PetscFunctionReturn(0);
143 }
144 
145 // Implicit time-stepper function setup
146 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G,
147                             void *user_data) {
148   User              user = *(User *)user_data;
149   const PetscScalar *q, *q_dot;
150   PetscScalar       *g;
151   Vec               Q_loc, Q_dot_loc, G_loc;
152   PetscMemType      q_mem_type, q_dot_mem_type, g_mem_type;
153   PetscErrorCode    ierr;
154   PetscFunctionBeginUser;
155 
156   // Update context field labels
157   if (user->phys->solution_time_label)
158     CeedOperatorContextSetDouble(user->op_ifunction,
159                                  user->phys->solution_time_label, &t);
160   if (user->phys->timestep_size_label) {
161     PetscScalar dt;
162     ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
163     CeedOperatorContextSetDouble(user->op_ifunction,
164                                  user->phys->timestep_size_label, &dt);
165   }
166 
167   // Get local vectors
168   ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
169   ierr = DMGetLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr);
170   ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
171 
172   // Global-to-local
173   ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr);
174   ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr);
175   ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t,
176                                     NULL, NULL, NULL); CHKERRQ(ierr);
177   ierr = VecZeroEntries(Q_dot_loc); CHKERRQ(ierr);
178   ierr = DMGlobalToLocal(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc);
179   CHKERRQ(ierr);
180   ierr = VecZeroEntries(G_loc); CHKERRQ(ierr);
181 
182   // Place PETSc vectors in CEED vectors
183   ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr);
184   ierr = VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type);
185   CHKERRQ(ierr);
186   ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr);
187   CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER,
188                      (PetscScalar *)q);
189   CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type),
190                      CEED_USE_POINTER, (PetscScalar *)q_dot);
191   CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g);
192 
193   // Apply CEED operator
194   CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed,
195                     CEED_REQUEST_IMMEDIATE);
196 
197   // Restore vectors
198   CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL);
199   CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL);
200   CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL);
201   ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr);
202   ierr = VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot); CHKERRQ(ierr);
203   ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr);
204 
205   // Local-to-Global
206   ierr = VecZeroEntries(G); CHKERRQ(ierr);
207   ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr);
208 
209   // Restore vectors
210   ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
211   ierr = DMRestoreLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr);
212   ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
213 
214   PetscFunctionReturn(0);
215 }
216 
217 static PetscErrorCode MatMult_NS_IJacobian(Mat J, Vec Q, Vec G) {
218   User user;
219   const PetscScalar *q;
220   PetscScalar       *g;
221   Vec               Q_loc, G_loc;
222   PetscMemType      q_mem_type, g_mem_type;
223   PetscErrorCode    ierr;
224   PetscFunctionBeginUser;
225   MatShellGetContext(J, &user);
226   // Get local vectors
227   ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
228   ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
229 
230   // Global-to-local
231   ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr);
232   ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr);
233   ierr = VecZeroEntries(G_loc); CHKERRQ(ierr);
234 
235   // Place PETSc vectors in CEED vectors
236   ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr);
237   ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr);
238   CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER,
239                      (PetscScalar *)q);
240   CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g);
241 
242   // Apply CEED operator
243   CeedOperatorApply(user->op_ijacobian, user->q_ceed, user->g_ceed,
244                     CEED_REQUEST_IMMEDIATE);
245 
246   // Restore vectors
247   CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL);
248   CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL);
249   ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr);
250   ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr);
251 
252   // Local-to-Global
253   ierr = VecZeroEntries(G); CHKERRQ(ierr);
254   ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr);
255 
256   // Restore vectors
257   ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
258   ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr);
259   PetscFunctionReturn(0);
260 }
261 
262 PetscErrorCode MatGetDiagonal_NS_IJacobian(Mat A, Vec D) {
263   User user;
264   Vec D_loc;
265   PetscScalar *d;
266   PetscMemType mem_type;
267 
268   PetscFunctionBeginUser;
269   MatShellGetContext(A, &user);
270   PetscCall(DMGetLocalVector(user->dm, &D_loc));
271   PetscCall(VecGetArrayAndMemType(D_loc, &d, &mem_type));
272   CeedVectorSetArray(user->g_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, d);
273   CeedOperatorLinearAssembleDiagonal(user->op_ijacobian, user->g_ceed,
274                                      CEED_REQUEST_IMMEDIATE);
275   CeedVectorTakeArray(user->g_ceed, MemTypeP2C(mem_type), NULL);
276   PetscCall(VecRestoreArrayAndMemType(D_loc, &d));
277   PetscCall(VecZeroEntries(D));
278   PetscCall(DMLocalToGlobal(user->dm, D_loc, ADD_VALUES, D));
279   PetscCall(DMRestoreLocalVector(user->dm, &D_loc));
280   VecViewFromOptions(D, NULL, "-diag_vec_view");
281   PetscFunctionReturn(0);
282 }
283 
284 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot,
285                                 PetscReal shift, Mat J, Mat J_pre,
286                                 void *user_data) {
287   User user = *(User *)user_data;
288   PetscBool J_is_shell, J_pre_is_shell;
289   PetscFunctionBeginUser;
290   if (user->phys->ijacobian_time_shift_label)
291     CeedOperatorContextSetDouble(user->op_ijacobian,
292                                  user->phys->ijacobian_time_shift_label, &shift);
293   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
294   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
295   Vec coo_vec = NULL;
296   PetscCall(PetscObjectTypeCompare((PetscObject)J, MATSHELL, &J_is_shell));
297   PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATSHELL,
298                                    &J_pre_is_shell));
299   if (!user->matrices_set_up) {
300     if (J_is_shell) {
301       PetscCall(MatShellSetContext(J, user));
302       PetscCall(MatShellSetOperation(J, MATOP_MULT,
303                                      (void (*)(void))MatMult_NS_IJacobian));
304       PetscCall(MatShellSetOperation(J, MATOP_GET_DIAGONAL,
305                                      (void (*)(void))MatGetDiagonal_NS_IJacobian));
306       PetscCall(MatSetUp(J));
307     }
308     if (!J_pre_is_shell) {
309       PetscCount ncoo;
310       PetscInt *rows, *cols;
311       PetscCall(CeedOperatorLinearAssembleSymbolic(user->op_ijacobian, &ncoo, &rows,
312                 &cols));
313       PetscCall(MatSetPreallocationCOOLocal(J_pre, ncoo, rows, cols));
314       free(rows);
315       free(cols);
316       CeedVectorCreate(user->ceed, ncoo, &user->coo_values);
317       user->matrices_set_up = true;
318       VecCreateSeq(PETSC_COMM_WORLD, ncoo, &coo_vec);
319     }
320   }
321   if (!J_pre_is_shell) {
322     CeedMemType mem_type = CEED_MEM_HOST;
323     const PetscScalar *values;
324     MatType mat_type;
325     PetscCall(MatGetType(J_pre, &mat_type));
326     //if (strstr(mat_type, "kokkos") || strstr(mat_type, "cusparse")) mem_type = CEED_MEM_DEVICE;
327     CeedOperatorLinearAssemble(user->op_ijacobian, user->coo_values);
328     CeedVectorGetArrayRead(user->coo_values, mem_type, &values);
329     if (coo_vec) {
330       VecPlaceArray(coo_vec, values);
331       VecViewFromOptions(coo_vec, NULL, "-coo_vec_view");
332       VecDestroy(&coo_vec);
333     }
334     PetscCall(MatSetValuesCOO(J_pre, values, INSERT_VALUES));
335     CeedVectorRestoreArrayRead(user->coo_values, &values);
336   }
337   PetscFunctionReturn(0);
338 }
339 
340 // User provided TS Monitor
341 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time,
342                             Vec Q, void *ctx) {
343   User           user = ctx;
344   Vec            Q_loc;
345   char           file_path[PETSC_MAX_PATH_LEN];
346   PetscViewer    viewer;
347   PetscErrorCode ierr;
348   PetscFunctionBeginUser;
349 
350   // Print every 'output_freq' steps
351   if (step_no % user->app_ctx->output_freq != 0)
352     PetscFunctionReturn(0);
353 
354   // Set up output
355   ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
356   ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr);
357   ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr);
358   ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr);
359 
360   // Output
361   ierr = PetscSNPrintf(file_path, sizeof file_path,
362                        "%s/ns-%03" PetscInt_FMT ".vtu",
363                        user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps);
364   CHKERRQ(ierr);
365   ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path,
366                             FILE_MODE_WRITE, &viewer); CHKERRQ(ierr);
367   ierr = VecView(Q_loc, viewer); CHKERRQ(ierr);
368   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
369   if (user->dm_viz) {
370     Vec         Q_refined, Q_refined_loc;
371     char        file_path_refined[PETSC_MAX_PATH_LEN];
372     PetscViewer viewer_refined;
373 
374     ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr);
375     ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr);
376     ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined");
377     CHKERRQ(ierr);
378     ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr);
379     ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr);
380     ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc);
381     CHKERRQ(ierr);
382     ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined,
383                          "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir,
384                          step_no + user->app_ctx->cont_steps);
385     CHKERRQ(ierr);
386     ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined),
387                               file_path_refined, FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr);
388     ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr);
389     ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr);
390     ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr);
391     ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr);
392   }
393   ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr);
394 
395   // Save data in a binary file for continuation of simulations
396   ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin",
397                        user->app_ctx->output_dir); CHKERRQ(ierr);
398   ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer);
399   CHKERRQ(ierr);
400   ierr = VecView(Q, viewer); CHKERRQ(ierr);
401   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
402 
403   // Save time stamp
404   // Dimensionalize time back
405   time /= user->units->second;
406   ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin",
407                        user->app_ctx->output_dir); CHKERRQ(ierr);
408   ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer);
409   CHKERRQ(ierr);
410   #if PETSC_VERSION_GE(3,13,0)
411   ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL);
412   #else
413   ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true);
414   #endif
415   CHKERRQ(ierr);
416   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
417 
418   PetscFunctionReturn(0);
419 }
420 
421 // TS: Create, setup, and solve
422 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys,
423                           Vec *Q, PetscScalar *f_time, TS *ts) {
424   MPI_Comm       comm = user->comm;
425   TSAdapt        adapt;
426   PetscScalar    final_time;
427   PetscErrorCode ierr;
428   PetscFunctionBeginUser;
429 
430   ierr = TSCreate(comm, ts); CHKERRQ(ierr);
431   ierr = TSSetDM(*ts, dm); CHKERRQ(ierr);
432   if (phys->implicit) {
433     ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr);
434     if (user->op_ifunction) {
435       ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr);
436     } else { // Implicit integrators can fall back to using an RHSFunction
437       ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr);
438     }
439     if (user->op_ijacobian) {
440       ierr = DMTSSetIJacobian(dm, FormIJacobian_NS, &user); CHKERRQ(ierr);
441     }
442   } else {
443     if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL,
444                                  "Problem does not provide RHSFunction");
445     ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr);
446     ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr);
447     ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr);
448   }
449   ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr);
450   ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr);
451   ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr);
452   if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);}
453   ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr);
454   ierr = TSAdaptSetStepLimits(adapt,
455                               1.e-12 * user->units->second,
456                               1.e2 * user->units->second); CHKERRQ(ierr);
457   ierr = TSSetFromOptions(*ts); CHKERRQ(ierr);
458   if (!app_ctx->cont_steps) { // print initial condition
459     if (!app_ctx->test_mode) {
460       ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr);
461     }
462   } else { // continue from time of last output
463     PetscReal   time;
464     PetscInt    count;
465     PetscViewer viewer;
466     char        file_path[PETSC_MAX_PATH_LEN];
467     ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin",
468                          app_ctx->output_dir); CHKERRQ(ierr);
469     ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer);
470     CHKERRQ(ierr);
471     ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL);
472     CHKERRQ(ierr);
473     ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
474     ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr);
475   }
476   if (!app_ctx->test_mode) {
477     ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr);
478   }
479 
480   // Solve
481   double start, cpu_time_used;
482   start = MPI_Wtime();
483   ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr);
484   ierr = TSSolve(*ts, *Q); CHKERRQ(ierr);
485   cpu_time_used = MPI_Wtime() - start;
486   ierr = TSGetSolveTime(*ts, &final_time); CHKERRQ(ierr);
487   *f_time = final_time;
488   ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN,
489                        comm); CHKERRQ(ierr);
490   if (!app_ctx->test_mode) {
491     ierr = PetscPrintf(PETSC_COMM_WORLD,
492                        "Time taken for solution (sec): %g\n",
493                        (double)cpu_time_used); CHKERRQ(ierr);
494   }
495   PetscFunctionReturn(0);
496 }
497