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