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