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