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 = TSGetTime(ts,&t);CHKERRQ(ierr); 66 67 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 68 ierr = VecDuplicate(U,&X);CHKERRQ(ierr); 69 ierr = TSEvaluateStep(ts,order,X,NULL);CHKERRQ(ierr); 70 ierr = VecGetArrayRead(X,&xx);CHKERRQ(ierr); 71 for (i = 0; i < order; i++) { 72 PetscReal error = PetscAbsReal(PetscRealPart(xx[i]) - PetscPowReal(t,i+1)); 73 if (error > eps) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bad solution, error %g, %s '%s'",(double)error,type,subtype); 74 } 75 ierr = VecRestoreArrayRead(X,&xx);CHKERRQ(ierr); 76 ierr = VecDestroy(&X);CHKERRQ(ierr); 77 78 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 79 ierr = VecDuplicate(U,&Y);CHKERRQ(ierr); 80 ierr = TSEvaluateStep(ts,order-1,Y,&done);CHKERRQ(ierr); 81 ierr = VecGetArrayRead(Y,&yy);CHKERRQ(ierr); 82 for (i = 0; done && i < order-1; i++) { 83 PetscReal error = PetscAbsReal(PetscRealPart(yy[i]) - PetscPowReal(t,i+1)); 84 if (error > eps) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bad estimator, error %g, %s '%s'",(double)error,type,subtype); 85 } 86 ierr = VecRestoreArrayRead(Y,&yy);CHKERRQ(ierr); 87 ierr = VecDestroy(&Y);CHKERRQ(ierr); 88 89 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 90 ierr = TSAdaptSetCheckStage(adapt,TestCheckStage);CHKERRQ(ierr); 91 ierr = TSSetErrorIfStepFails(ts,PETSC_FALSE);CHKERRQ(ierr); 92 ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr); 93 ierr = TSSetTime(ts,0);CHKERRQ(ierr); 94 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 95 ierr = TSSolve(ts,NULL);CHKERRQ(ierr); 96 ierr = TSAdaptSetCheckStage(adapt,NULL);CHKERRQ(ierr); 97 ierr = TSSetErrorIfStepFails(ts,PETSC_TRUE);CHKERRQ(ierr); 98 ierr = TSGetStepNumber(ts,&step);CHKERRQ(ierr); 99 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 100 if (step != 2) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bad step number %D, %s '%s'",step,type,subtype); 101 if (reason != TS_DIVERGED_STEP_REJECTED) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bad reason %s, %s '%s'",TSConvergedReasons[reason],type,subtype); 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode TestTSRK(TS ts,TSRKType type) 106 { 107 PetscInt order; 108 TSAdapt adapt; 109 PetscBool rk1,rk3,rk4; 110 TSAdaptType adapttype; 111 char savetype[32]; 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 ierr = TSRKSetType(ts,type);CHKERRQ(ierr); 116 ierr = TSRKGetType(ts,&type);CHKERRQ(ierr); 117 ierr = TSRKGetOrder(ts,&order);CHKERRQ(ierr); 118 119 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 120 ierr = TSAdaptGetType(adapt,&adapttype);CHKERRQ(ierr); 121 ierr = PetscStrncpy(savetype,adapttype,sizeof(savetype));CHKERRQ(ierr); 122 ierr = PetscStrcmp(type,TSRK1FE,&rk1);CHKERRQ(ierr); 123 ierr = PetscStrcmp(type,TSRK3,&rk3);CHKERRQ(ierr); 124 ierr = PetscStrcmp(type,TSRK4,&rk4);CHKERRQ(ierr); 125 if (rk1 || rk3 || rk4) {ierr = TSAdaptSetType(adapt,TSADAPTNONE);CHKERRQ(ierr);} 126 127 ierr = TestExplicitTS(ts,order,type);CHKERRQ(ierr); 128 129 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 130 ierr = TSAdaptSetType(adapt,savetype);CHKERRQ(ierr); 131 PetscFunctionReturn(0); 132 } 133 134 int main(int argc, char *argv[]) 135 { 136 TS ts; 137 Vec X; 138 PetscInt N = 9; 139 PetscErrorCode ierr; 140 141 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 142 143 ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr); 144 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr); 145 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,NULL);CHKERRQ(ierr); 146 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&X);CHKERRQ(ierr); 147 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 148 ierr = VecDestroy(&X);CHKERRQ(ierr); 149 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 150 151 ierr = TestTSRK(ts,TSRK1FE);CHKERRQ(ierr); 152 ierr = TestTSRK(ts,TSRK2A);CHKERRQ(ierr); 153 ierr = TestTSRK(ts,TSRK3);CHKERRQ(ierr); 154 ierr = TestTSRK(ts,TSRK3BS);CHKERRQ(ierr); 155 ierr = TestTSRK(ts,TSRK4);CHKERRQ(ierr); 156 ierr = TestTSRK(ts,TSRK5F);CHKERRQ(ierr); 157 ierr = TestTSRK(ts,TSRK5DP);CHKERRQ(ierr); 158 ierr = TestTSRK(ts,TSRK5BS);CHKERRQ(ierr); 159 ierr = TestTSRK(ts,TSRK6VR);CHKERRQ(ierr); 160 ierr = TestTSRK(ts,TSRK7VR);CHKERRQ(ierr); 161 ierr = TestTSRK(ts,TSRK8VR);CHKERRQ(ierr); 162 163 ierr = TSRollBack(ts);CHKERRQ(ierr); 164 ierr = TSDestroy(&ts);CHKERRQ(ierr); 165 166 ierr = PetscFinalize(); 167 return ierr; 168 } 169 170 /*TEST 171 172 testset: 173 output_file: output/ex14.out 174 test: 175 suffix: 0 176 test: 177 suffix: 1 178 args: -ts_adapt_type none 179 test: 180 suffix: 2 181 args: -ts_adapt_type basic 182 test: 183 suffix: 3 184 args: -ts_adapt_type dsp 185 186 TEST*/ 187