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 CHKERRQ(VecGetLocalSize(X,&n)); 13 CHKERRQ(VecGetArrayRead(X,&xx)); 14 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(X,&xx)); 22 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TSGetType(ts,&type)); 54 CHKERRQ(TSGetSolution(ts,&U)); 55 CHKERRQ(VecZeroEntries(U)); 56 CHKERRQ(TSSetStepNumber(ts,0)); 57 CHKERRQ(TSSetTime(ts,0)); 58 CHKERRQ(TSSetTimeStep(ts,dt)); 59 CHKERRQ(TSSetMaxTime(ts,Tf)); 60 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 61 CHKERRQ(TSSolve(ts,NULL)); 62 CHKERRQ(TSRollBack(ts)); 63 CHKERRQ(TSSolve(ts,NULL)); 64 CHKERRQ(TSGetTime(ts,&t)); 65 66 CHKERRQ(TSGetSolution(ts,&U)); 67 CHKERRQ(VecDuplicate(U,&X)); 68 CHKERRQ(TSEvaluateStep(ts,order,X,NULL)); 69 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(X,&xx)); 75 CHKERRQ(VecDestroy(&X)); 76 77 CHKERRQ(TSGetSolution(ts,&U)); 78 CHKERRQ(VecDuplicate(U,&Y)); 79 CHKERRQ(TSEvaluateStep(ts,order-1,Y,&done)); 80 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(Y,&yy)); 86 CHKERRQ(VecDestroy(&Y)); 87 88 CHKERRQ(TSGetAdapt(ts,&adapt)); 89 CHKERRQ(TSAdaptSetCheckStage(adapt,TestCheckStage)); 90 CHKERRQ(TSSetErrorIfStepFails(ts,PETSC_FALSE)); 91 CHKERRQ(TSSetStepNumber(ts,0)); 92 CHKERRQ(TSSetTime(ts,0)); 93 CHKERRQ(TSSetTimeStep(ts,dt)); 94 CHKERRQ(TSSolve(ts,NULL)); 95 CHKERRQ(TSAdaptSetCheckStage(adapt,NULL)); 96 CHKERRQ(TSSetErrorIfStepFails(ts,PETSC_TRUE)); 97 CHKERRQ(TSGetStepNumber(ts,&step)); 98 CHKERRQ(TSGetConvergedReason(ts,&reason)); 99 PetscCheck(step == 2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bad step number %D, %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 CHKERRQ(TSRKSetType(ts,type)); 114 CHKERRQ(TSRKGetType(ts,&type)); 115 CHKERRQ(TSRKGetOrder(ts,&order)); 116 117 CHKERRQ(TSGetAdapt(ts,&adapt)); 118 CHKERRQ(TSAdaptGetType(adapt,&adapttype)); 119 CHKERRQ(PetscStrncpy(savetype,adapttype,sizeof(savetype))); 120 CHKERRQ(PetscStrcmp(type,TSRK1FE,&rk1)); 121 CHKERRQ(PetscStrcmp(type,TSRK3,&rk3)); 122 CHKERRQ(PetscStrcmp(type,TSRK4,&rk4)); 123 if (rk1 || rk3 || rk4) CHKERRQ(TSAdaptSetType(adapt,TSADAPTNONE)); 124 125 CHKERRQ(TestExplicitTS(ts,order,type)); 126 127 CHKERRQ(TSGetAdapt(ts,&adapt)); 128 CHKERRQ(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 PetscErrorCode ierr; 138 139 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 140 141 CHKERRQ(TSCreate(PETSC_COMM_SELF,&ts)); 142 CHKERRQ(TSSetType(ts,TSRK)); 143 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 144 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,N,&X)); 145 CHKERRQ(TSSetSolution(ts,X)); 146 CHKERRQ(VecDestroy(&X)); 147 CHKERRQ(TSSetFromOptions(ts)); 148 149 CHKERRQ(TestTSRK(ts,TSRK1FE)); 150 CHKERRQ(TestTSRK(ts,TSRK2A)); 151 CHKERRQ(TestTSRK(ts,TSRK3)); 152 CHKERRQ(TestTSRK(ts,TSRK3BS)); 153 CHKERRQ(TestTSRK(ts,TSRK4)); 154 CHKERRQ(TestTSRK(ts,TSRK5F)); 155 CHKERRQ(TestTSRK(ts,TSRK5DP)); 156 CHKERRQ(TestTSRK(ts,TSRK5BS)); 157 CHKERRQ(TestTSRK(ts,TSRK6VR)); 158 CHKERRQ(TestTSRK(ts,TSRK7VR)); 159 CHKERRQ(TestTSRK(ts,TSRK8VR)); 160 161 CHKERRQ(TSRollBack(ts)); 162 CHKERRQ(TSDestroy(&ts)); 163 164 ierr = PetscFinalize(); 165 return ierr; 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