setupts.c (ed9ed3dea34fef5853dd3f23fdf30a7807a29b2c) setupts.c (91c97f41bd25cf2fbada21392acdb50722cb486e)
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

--- 184 unchanged lines hidden (view full) ---

193 PetscCall(VecZeroEntries(G));
194 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G));
195
196 // Restore vectors
197 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
198 PetscFunctionReturn(PETSC_SUCCESS);
199}
200
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

--- 184 unchanged lines hidden (view full) ---

193 PetscCall(VecZeroEntries(G));
194 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G));
195
196 // Restore vectors
197 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc));
198 PetscFunctionReturn(PETSC_SUCCESS);
199}
200
201static PetscErrorCode FormPreallocation(User user, PetscBool pbdiagonal, Mat J, CeedVector *coo_values) {
202 PetscCount ncoo;
203 PetscInt *rows_petsc, *cols_petsc;
204 CeedInt *rows_ceed, *cols_ceed;
205
206 PetscFunctionBeginUser;
207 if (pbdiagonal) {
208 PetscCallCeed(user->ceed, CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(user->op_ijacobian, &ncoo, &rows_ceed, &cols_ceed));
209 } else {
210 PetscCallCeed(user->ceed, CeedOperatorLinearAssembleSymbolic(user->op_ijacobian, &ncoo, &rows_ceed, &cols_ceed));
211 }
212 PetscCall(IntArrayC2P(ncoo, &rows_ceed, &rows_petsc));
213 PetscCall(IntArrayC2P(ncoo, &cols_ceed, &cols_petsc));
214 PetscCall(MatSetPreallocationCOOLocal(J, ncoo, rows_petsc, cols_petsc));
215 free(rows_petsc);
216 free(cols_petsc);
217 PetscCallCeed(user->ceed, CeedVectorCreate(user->ceed, ncoo, coo_values));
218 PetscFunctionReturn(PETSC_SUCCESS);
219}
220
221static PetscErrorCode FormSetValues(User user, PetscBool pbdiagonal, Mat J, CeedVector coo_values) {
222 CeedMemType mem_type = CEED_MEM_HOST;
223 const PetscScalar *values;
224 MatType mat_type;
225
226 PetscFunctionBeginUser;
227 PetscCall(MatGetType(J, &mat_type));
228 if (strstr(mat_type, "kokkos") || strstr(mat_type, "cusparse")) mem_type = CEED_MEM_DEVICE;
229 if (pbdiagonal) {
230 PetscCall(PetscLogEventBegin(FLUIDS_CeedOperatorAssemblePointBlockDiagonal, J, 0, 0, 0));
231 PetscCall(PetscLogGpuTimeBegin());
232 PetscCallCeed(user->ceed, CeedOperatorLinearAssemblePointBlockDiagonal(user->op_ijacobian, coo_values, CEED_REQUEST_IMMEDIATE));
233 PetscCall(PetscLogGpuTimeEnd());
234 PetscCall(PetscLogEventEnd(FLUIDS_CeedOperatorAssemblePointBlockDiagonal, J, 0, 0, 0));
235 } else {
236 PetscCall(PetscLogEventBegin(FLUIDS_CeedOperatorAssemble, J, 0, 0, 0));
237 PetscCall(PetscLogGpuTimeBegin());
238 PetscCallCeed(user->ceed, CeedOperatorLinearAssemble(user->op_ijacobian, coo_values));
239 PetscCall(PetscLogGpuTimeEnd());
240 PetscCall(PetscLogEventEnd(FLUIDS_CeedOperatorAssemble, J, 0, 0, 0));
241 }
242 PetscCallCeed(user->ceed, CeedVectorGetArrayRead(coo_values, mem_type, &values));
243 PetscCall(MatSetValuesCOO(J, values, INSERT_VALUES));
244 PetscCallCeed(user->ceed, CeedVectorRestoreArrayRead(coo_values, &values));
245 PetscFunctionReturn(PETSC_SUCCESS);
246}
247
248PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) {
249 User user = *(User *)user_data;
250 Ceed ceed = user->ceed;
201PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) {
202 User user = *(User *)user_data;
203 Ceed ceed = user->ceed;
251 PetscBool J_is_shell, J_is_mffd, J_pre_is_shell;
204 PetscBool J_is_matceed, J_is_mffd, J_pre_is_matceed, J_pre_is_mffd;
252
253 PetscFunctionBeginUser;
205
206 PetscFunctionBeginUser;
254 if (user->phys->ijacobian_time_shift_label)
255 PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_ijacobian, user->phys->ijacobian_time_shift_label, &shift));
256 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd));
207 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd));
257 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATSHELL, &J_is_shell));
258 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATSHELL, &J_pre_is_shell));
259 if (!user->matrices_set_up) {
260 if (J_is_shell) {
261 OperatorApplyContext op_ijacobian_ctx;
262 OperatorApplyContextCreate(user->dm, user->dm, user->ceed, user->op_ijacobian, user->q_ceed, user->g_ceed, user->Q_dot_loc, NULL,
263 &op_ijacobian_ctx);
264 PetscCall(CreateMatShell_Ceed(op_ijacobian_ctx, &J));
265 PetscCall(MatSetUp(J));
266 }
267 if (!J_pre_is_shell) {
268 PetscCall(FormPreallocation(user, user->app_ctx->pmat_pbdiagonal, J_pre, &user->coo_values_pmat));
269 }
270 if (J != J_pre && !J_is_shell && !J_is_mffd) {
271 PetscCall(FormPreallocation(user, PETSC_FALSE, J, &user->coo_values_amat));
272 }
273 user->matrices_set_up = true;
208 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATCEED, &J_is_matceed));
209 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATMFFD, &J_pre_is_mffd));
210 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATCEED, &J_pre_is_matceed));
211 if (user->phys->ijacobian_time_shift_label) {
212 CeedOperator op_ijacobian;
213
214 PetscCall(MatCeedGetCeedOperators(user->mat_ijacobian, &op_ijacobian, NULL));
215 PetscCallCeed(ceed, CeedOperatorSetContextDouble(op_ijacobian, user->phys->ijacobian_time_shift_label, &shift));
274 }
216 }
275 if (!J_pre_is_shell) {
276 PetscCall(FormSetValues(user, user->app_ctx->pmat_pbdiagonal, J_pre, user->coo_values_pmat));
277 }
278 if (user->coo_values_amat) {
279 PetscCall(FormSetValues(user, PETSC_FALSE, J, user->coo_values_amat));
280 } else if (J_is_mffd) {
217
218 if (J_is_matceed || J_is_mffd) {
281 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
282 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
219 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
220 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
221 } else PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J));
222
223 if (J_pre_is_matceed && J != J_pre) {
224 PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY));
225 PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY));
226 } else if (!J_pre_is_matceed && !J_pre_is_mffd && J != J_pre) {
227 PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J_pre));
283 }
284 PetscFunctionReturn(PETSC_SUCCESS);
285}
286
287PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) {
288 Vec Q_loc;
289 char file_path[PETSC_MAX_PATH_LEN];
290 PetscViewer viewer;

--- 123 unchanged lines hidden (view full) ---

414 PetscCall(TSSetApplicationContext(*ts, user));
415 if (phys->implicit) {
416 PetscCall(TSSetType(*ts, TSBDF));
417 if (user->op_ifunction) {
418 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user));
419 } else { // Implicit integrators can fall back to using an RHSFunction
420 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
421 }
228 }
229 PetscFunctionReturn(PETSC_SUCCESS);
230}
231
232PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) {
233 Vec Q_loc;
234 char file_path[PETSC_MAX_PATH_LEN];
235 PetscViewer viewer;

--- 123 unchanged lines hidden (view full) ---

359 PetscCall(TSSetApplicationContext(*ts, user));
360 if (phys->implicit) {
361 PetscCall(TSSetType(*ts, TSBDF));
362 if (user->op_ifunction) {
363 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user));
364 } else { // Implicit integrators can fall back to using an RHSFunction
365 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
366 }
422 if (user->op_ijacobian) {
367 if (user->mat_ijacobian) {
423 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user));
368 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user));
424 if (app_ctx->amat_type) {
425 Mat Pmat, Amat;
426 PetscCall(DMCreateMatrix(dm, &Pmat));
427 PetscCall(DMSetMatType(dm, app_ctx->amat_type));
428 PetscCall(DMCreateMatrix(dm, &Amat));
429 PetscCall(TSSetIJacobian(*ts, Amat, Pmat, NULL, NULL));
430 PetscCall(MatDestroy(&Amat));
431 PetscCall(MatDestroy(&Pmat));
432 }
433 }
434 } else {
435 PetscCheck(user->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction");
436 PetscCall(TSSetType(*ts, TSRK));
437 PetscCall(TSRKSetType(*ts, TSRK5F));
438 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
439 }
440 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second));
441 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
442 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE));
443 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second));
444 PetscCall(TSGetAdapt(*ts, &adapt));
445 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second));
446 PetscCall(TSSetFromOptions(*ts));
369 }
370 } else {
371 PetscCheck(user->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction");
372 PetscCall(TSSetType(*ts, TSRK));
373 PetscCall(TSRKSetType(*ts, TSRK5F));
374 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user));
375 }
376 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second));
377 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
378 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE));
379 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second));
380 PetscCall(TSGetAdapt(*ts, &adapt));
381 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second));
382 PetscCall(TSSetFromOptions(*ts));
383 if (user->mat_ijacobian) {
384 if (app_ctx->amat_type && !strcmp(app_ctx->amat_type, MATSHELL)) {
385 PetscBool use_matceed_pmat;
386 SNES snes;
387 KSP ksp;
388 PC pc;
389 PCType pc_type;
390 Mat Pmat;
391
392 PetscCall(TSGetSNES(*ts, &snes));
393 PetscCall(SNESGetKSP(snes, &ksp));
394 PetscCall(KSPGetPC(ksp, &pc));
395 PetscCall(PCGetType(pc, &pc_type));
396 PetscCall(PetscStrcmpAny(pc_type, &use_matceed_pmat, PCJACOBI, PCVPBJACOBI, PCPBJACOBI, ""));
397 Pmat = use_matceed_pmat ? user->mat_ijacobian : NULL;
398 PetscCall(TSSetIJacobian(*ts, user->mat_ijacobian, Pmat, NULL, NULL));
399 }
400 }
447 user->time_bc_set = -1.0; // require all BCs be updated
448 if (app_ctx->cont_steps) { // continue from previous timestep data
449 PetscInt count;
450 PetscViewer viewer;
451
452 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time
453 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer));
454 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL));

--- 75 unchanged lines hidden --- 		401 user->time_bc_set = -1.0; // require all BCs be updated
402 if (app_ctx->cont_steps) { // continue from previous timestep data
403 PetscInt count;
404 PetscViewer viewer;
405
406 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time
407 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer));
408 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL));

--- 75 unchanged lines hidden ---