1 static char help[] = "Solves the trivial ODE 2 du/dt = 1, u(0) = 0. \n\n"; 2 3 #include <petscts.h> 4 #include <petscpc.h> 5 6 PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *); 7 PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 8 PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 9 10 int main(int argc, char **argv) 11 { 12 TS ts; 13 Vec x; 14 Mat A; 15 PetscBool flg = PETSC_FALSE, usingimex; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 19 20 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 21 PetscCall(PetscOptionsGetBool(NULL, NULL, "-set_implicit", &flg, NULL)); 22 if (flg) PetscCall(TSSetEquationType(ts, TS_EQ_ODE_IMPLICIT)); 23 PetscCall(TSSetIFunction(ts, NULL, IFunction, NULL)); 24 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &usingimex)); 25 26 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 27 PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE)); 28 PetscCall(MatSetFromOptions(A)); 29 PetscCall(MatSetUp(A)); 30 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 31 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 32 PetscCall(TSSetIJacobian(ts, A, A, IJacobian, NULL)); 33 PetscCall(MatDestroy(&A)); 34 35 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 36 PetscCall(VecSetSizes(x, 1, PETSC_DECIDE)); 37 PetscCall(VecSetFromOptions(x)); 38 PetscCall(VecSetUp(x)); 39 PetscCall(TSSetSolution(ts, x)); 40 PetscCall(VecDestroy(&x)); 41 PetscCall(TSSetFromOptions(ts)); 42 43 /* Need to know if we are using an IMEX scheme to decide on the form 44 of the RHS function */ 45 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSARKIMEX, &usingimex)); 46 if (usingimex) { 47 PetscCall(TSARKIMEXGetFullyImplicit(ts, &flg)); 48 if (flg) usingimex = PETSC_FALSE; 49 } 50 PetscCall(TSSetStepNumber(ts, 0)); 51 PetscCall(TSSetTimeStep(ts, 1)); 52 PetscCall(TSSetTime(ts, 0)); 53 PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL)); 54 PetscCall(TSSetMaxSteps(ts, 3)); 55 56 PetscCall(TSSolve(ts, NULL)); 57 58 PetscCall(TSDestroy(&ts)); 59 PetscCall(PetscFinalize()); 60 return 0; 61 } 62 63 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx) 64 { 65 PetscBool usingimex = *(PetscBool *)ctx; 66 67 PetscFunctionBeginUser; 68 PetscCall(VecSet(f, usingimex ? 0.5 : 1)); 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 72 PetscErrorCode IFunction(TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx) 73 { 74 PetscFunctionBeginUser; 75 PetscCall(VecCopy(xdot, f)); 76 PetscCall(VecScale(f, 2.0)); 77 PetscFunctionReturn(PETSC_SUCCESS); 78 } 79 80 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec x, Vec xdot, PetscReal shift, Mat A, Mat B, void *ctx) 81 { 82 PetscScalar j; 83 84 PetscFunctionBeginUser; 85 j = shift * 2.0; 86 PetscCall(MatSetValue(B, 0, 0, j, INSERT_VALUES)); 87 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 88 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 92 /*TEST 93 94 test: 95 suffix: arkimex_explicit_stage 96 requires: !defined(PETSCTEST_VALGRIND) defined(PETSC_USE_DEBUG) 97 args: -ts_type arkimex -petsc_ci_portable_error_output -error_output_stdout -set_implicit 98 filter: grep -E -v "(memory block|leaked context|not freed before MPI_Finalize|Could be the program crashed)" 99 100 test: 101 suffix: arkimex_implicit_stage 102 args: -ts_type arkimex -ts_arkimex_type {{3 l2}} -ts_monitor_solution -ts_monitor 103 104 TEST*/ 105