1 static char help[] = "Small ODE to test TS accuracy.\n"; 2 3 /* 4 The ODE 5 u1_t = cos(t), 6 u2_t = sin(u2) 7 with analytical solution 8 u1(t) = sin(t), 9 u2(t) = 2 * atan(exp(t) * tan(0.5)) 10 is used to test the accuracy of TS schemes. 11 */ 12 13 #include <petscts.h> 14 15 /* 16 Defines the ODE passed to the ODE solver in explicit form: U_t = F(U) 17 */ 18 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s) { 19 PetscScalar *f; 20 const PetscScalar *u; 21 22 PetscFunctionBeginUser; 23 PetscCall(VecGetArrayRead(U, &u)); 24 PetscCall(VecGetArray(F, &f)); 25 26 f[0] = PetscCosReal(t); 27 f[1] = PetscSinReal(u[1]); 28 29 PetscCall(VecRestoreArrayRead(U, &u)); 30 PetscCall(VecRestoreArray(F, &f)); 31 PetscFunctionReturn(0); 32 } 33 34 /* 35 Defines the exact solution. 36 */ 37 static PetscErrorCode ExactSolution(PetscReal t, Vec U) { 38 PetscScalar *u; 39 40 PetscFunctionBeginUser; 41 PetscCall(VecGetArray(U, &u)); 42 43 u[0] = PetscSinReal(t); 44 u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5)); 45 46 PetscCall(VecRestoreArray(U, &u)); 47 PetscFunctionReturn(0); 48 } 49 50 int main(int argc, char **argv) { 51 TS ts; /* ODE integrator */ 52 Vec U; /* numerical solution will be stored here */ 53 Vec Uex; /* analytical (exact) solution will be stored here */ 54 PetscMPIInt size; 55 PetscInt n = 2; 56 PetscScalar *u; 57 PetscReal t, final_time = 1.0, dt = 0.25; 58 PetscReal error; 59 TSAdapt adapt; 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Initialize program 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 PetscFunctionBeginUser; 65 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 66 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 67 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 68 69 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 70 Create timestepping solver context 71 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 72 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 73 PetscCall(TSSetType(ts, TSROSW)); 74 75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 Set ODE routines 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 79 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Set initial conditions 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 PetscCall(VecCreate(PETSC_COMM_WORLD, &U)); 85 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE)); 86 PetscCall(VecSetUp(U)); 87 PetscCall(VecGetArray(U, &u)); 88 u[0] = 0.0; 89 u[1] = 1.0; 90 PetscCall(VecRestoreArray(U, &u)); 91 PetscCall(TSSetSolution(ts, U)); 92 93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94 Set solver options 95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 96 PetscCall(TSSetSaveTrajectory(ts)); 97 PetscCall(TSSetMaxTime(ts, final_time)); 98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 99 PetscCall(TSSetTimeStep(ts, dt)); 100 /* The adaptive time step controller is forced to take constant time steps. */ 101 PetscCall(TSGetAdapt(ts, &adapt)); 102 PetscCall(TSAdaptSetType(adapt, TSADAPTNONE)); 103 104 PetscCall(TSSetFromOptions(ts)); 105 106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 107 Run timestepping solver and compute error 108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 109 PetscCall(TSSolve(ts, U)); 110 PetscCall(TSGetTime(ts, &t)); 111 112 if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) { 113 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference of %g between the prescribed final time %g and the actual final time.\n", (double)(final_time - t), (double)final_time)); 114 } 115 PetscCall(VecDuplicate(U, &Uex)); 116 PetscCall(ExactSolution(t, Uex)); 117 118 PetscCall(VecAYPX(Uex, -1.0, U)); 119 PetscCall(VecNorm(Uex, NORM_2, &error)); 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error)); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Free work space. All PETSc objects should be destroyed when they are no longer needed. 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 PetscCall(VecDestroy(&U)); 126 PetscCall(VecDestroy(&Uex)); 127 PetscCall(TSDestroy(&ts)); 128 129 PetscCall(PetscFinalize()); 130 return 0; 131 } 132 133 /*TEST 134 135 test: 136 suffix: 3bs 137 args: -ts_type rk -ts_rk_type 3bs 138 requires: !single 139 140 test: 141 suffix: 5bs 142 args: -ts_type rk -ts_rk_type 5bs 143 requires: !single 144 145 test: 146 suffix: 5dp 147 args: -ts_type rk -ts_rk_type 5dp 148 requires: !single 149 150 test: 151 suffix: 6vr 152 args: -ts_type rk -ts_rk_type 6vr 153 requires: !single 154 155 test: 156 suffix: 7vr 157 args: -ts_type rk -ts_rk_type 7vr 158 requires: !single 159 160 test: 161 suffix: 8vr 162 args: -ts_type rk -ts_rk_type 8vr 163 requires: !single 164 165 TEST*/ 166