1 static char help[] = "Tests all TSRK types \n\n"; 2 3 #include <petscts.h> 4 5 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ctx) { 6 PetscInt i, n; 7 const PetscScalar *xx; 8 /* */ PetscScalar *ff; 9 10 PetscFunctionBeginUser; 11 PetscCall(VecGetLocalSize(X, &n)); 12 PetscCall(VecGetArrayRead(X, &xx)); 13 PetscCall(VecGetArray(F, &ff)); 14 15 if (n >= 1) ff[0] = 1; 16 for (i = 1; i < n; i++) ff[i] = (i + 1) * (xx[i - 1] + PetscPowReal(t, i)) / 2; 17 18 PetscCall(VecRestoreArrayRead(X, &xx)); 19 PetscCall(VecRestoreArray(F, &ff)); 20 PetscFunctionReturn(0); 21 } 22 23 PetscErrorCode TestCheckStage(TSAdapt adapt, TS ts, PetscReal t, Vec X, PetscBool *accept) { 24 PetscInt step; 25 26 PetscFunctionBeginUser; 27 PetscCall(TSGetStepNumber(ts, &step)); 28 *accept = (step >= 2) ? PETSC_FALSE : PETSC_TRUE; 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode TestExplicitTS(TS ts, PetscInt order, const char subtype[]) { 33 PetscInt i; 34 PetscReal t; 35 Vec U, X, Y; 36 TSType type; 37 PetscBool done; 38 const PetscScalar *xx; 39 const PetscScalar *yy; 40 const PetscReal Tf = 1; 41 const PetscReal dt = Tf / 8; 42 const PetscReal eps = 100 * PETSC_MACHINE_EPSILON; 43 TSAdapt adapt; 44 PetscInt step; 45 TSConvergedReason reason; 46 47 PetscFunctionBeginUser; 48 PetscCall(TSGetType(ts, &type)); 49 PetscCall(TSGetSolution(ts, &U)); 50 PetscCall(VecZeroEntries(U)); 51 PetscCall(TSSetStepNumber(ts, 0)); 52 PetscCall(TSSetTime(ts, 0)); 53 PetscCall(TSSetTimeStep(ts, dt)); 54 PetscCall(TSSetMaxTime(ts, Tf)); 55 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 56 PetscCall(TSSolve(ts, NULL)); 57 PetscCall(TSRollBack(ts)); 58 PetscCall(TSSolve(ts, NULL)); 59 PetscCall(TSGetTime(ts, &t)); 60 61 PetscCall(TSGetSolution(ts, &U)); 62 PetscCall(VecDuplicate(U, &X)); 63 PetscCall(TSEvaluateStep(ts, order, X, NULL)); 64 PetscCall(VecGetArrayRead(X, &xx)); 65 for (i = 0; i < order; i++) { 66 PetscReal error = PetscAbsReal(PetscRealPart(xx[i]) - PetscPowReal(t, i + 1)); 67 PetscCheck(error <= eps, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad solution, error %g, %s '%s'", (double)error, type, subtype); 68 } 69 PetscCall(VecRestoreArrayRead(X, &xx)); 70 PetscCall(VecDestroy(&X)); 71 72 PetscCall(TSGetSolution(ts, &U)); 73 PetscCall(VecDuplicate(U, &Y)); 74 PetscCall(TSEvaluateStep(ts, order - 1, Y, &done)); 75 PetscCall(VecGetArrayRead(Y, &yy)); 76 for (i = 0; done && i < order - 1; i++) { 77 PetscReal error = PetscAbsReal(PetscRealPart(yy[i]) - PetscPowReal(t, i + 1)); 78 PetscCheck(error <= eps, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad estimator, error %g, %s '%s'", (double)error, type, subtype); 79 } 80 PetscCall(VecRestoreArrayRead(Y, &yy)); 81 PetscCall(VecDestroy(&Y)); 82 83 PetscCall(TSGetAdapt(ts, &adapt)); 84 PetscCall(TSAdaptSetCheckStage(adapt, TestCheckStage)); 85 PetscCall(TSSetErrorIfStepFails(ts, PETSC_FALSE)); 86 PetscCall(TSSetStepNumber(ts, 0)); 87 PetscCall(TSSetTime(ts, 0)); 88 PetscCall(TSSetTimeStep(ts, dt)); 89 PetscCall(TSSolve(ts, NULL)); 90 PetscCall(TSAdaptSetCheckStage(adapt, NULL)); 91 PetscCall(TSSetErrorIfStepFails(ts, PETSC_TRUE)); 92 PetscCall(TSGetStepNumber(ts, &step)); 93 PetscCall(TSGetConvergedReason(ts, &reason)); 94 PetscCheck(step == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad step number %" PetscInt_FMT ", %s '%s'", step, type, subtype); 95 PetscCheck(reason == TS_DIVERGED_STEP_REJECTED, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad reason %s, %s '%s'", TSConvergedReasons[reason], type, subtype); 96 PetscFunctionReturn(0); 97 } 98 99 static PetscErrorCode TestTSRK(TS ts, TSRKType type) { 100 PetscInt order; 101 TSAdapt adapt; 102 PetscBool rk1, rk3, rk4; 103 TSAdaptType adapttype; 104 char savetype[32]; 105 106 PetscFunctionBeginUser; 107 PetscCall(TSRKSetType(ts, type)); 108 PetscCall(TSRKGetType(ts, &type)); 109 PetscCall(TSRKGetOrder(ts, &order)); 110 111 PetscCall(TSGetAdapt(ts, &adapt)); 112 PetscCall(TSAdaptGetType(adapt, &adapttype)); 113 PetscCall(PetscStrncpy(savetype, adapttype, sizeof(savetype))); 114 PetscCall(PetscStrcmp(type, TSRK1FE, &rk1)); 115 PetscCall(PetscStrcmp(type, TSRK3, &rk3)); 116 PetscCall(PetscStrcmp(type, TSRK4, &rk4)); 117 if (rk1 || rk3 || rk4) PetscCall(TSAdaptSetType(adapt, TSADAPTNONE)); 118 119 PetscCall(TestExplicitTS(ts, order, type)); 120 121 PetscCall(TSGetAdapt(ts, &adapt)); 122 PetscCall(TSAdaptSetType(adapt, savetype)); 123 PetscFunctionReturn(0); 124 } 125 126 int main(int argc, char *argv[]) { 127 TS ts; 128 Vec X; 129 PetscInt N = 9; 130 131 PetscFunctionBeginUser; 132 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 133 134 PetscCall(TSCreate(PETSC_COMM_SELF, &ts)); 135 PetscCall(TSSetType(ts, TSRK)); 136 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL)); 137 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &X)); 138 PetscCall(TSSetSolution(ts, X)); 139 PetscCall(VecDestroy(&X)); 140 PetscCall(TSSetFromOptions(ts)); 141 142 PetscCall(TestTSRK(ts, TSRK1FE)); 143 PetscCall(TestTSRK(ts, TSRK2A)); 144 PetscCall(TestTSRK(ts, TSRK3)); 145 PetscCall(TestTSRK(ts, TSRK3BS)); 146 PetscCall(TestTSRK(ts, TSRK4)); 147 PetscCall(TestTSRK(ts, TSRK5F)); 148 PetscCall(TestTSRK(ts, TSRK5DP)); 149 PetscCall(TestTSRK(ts, TSRK5BS)); 150 PetscCall(TestTSRK(ts, TSRK6VR)); 151 PetscCall(TestTSRK(ts, TSRK7VR)); 152 PetscCall(TestTSRK(ts, TSRK8VR)); 153 154 PetscCall(TSRollBack(ts)); 155 PetscCall(TSDestroy(&ts)); 156 157 PetscCall(PetscFinalize()); 158 return 0; 159 } 160 161 /*TEST 162 163 testset: 164 output_file: output/ex14.out 165 test: 166 suffix: 0 167 test: 168 suffix: 1 169 args: -ts_adapt_type none 170 test: 171 suffix: 2 172 args: -ts_adapt_type basic 173 test: 174 suffix: 3 175 args: -ts_adapt_type dsp 176 177 TEST*/ 178