1 /* 2 Formatted test for TS routines. 3 4 Solves U_t=F(t,u) 5 Where: 6 7 [2*u1+u2 8 F(t,u)= [u1+2*u2+u3 9 [ u2+2*u3 10 We can compare the solutions from euler, beuler and SUNDIALS to 11 see what is the difference. 12 13 */ 14 15 static char help[] = "Solves a linear ODE. \n\n"; 16 17 #include <petscts.h> 18 #include <petscpc.h> 19 20 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 21 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 22 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); 23 extern PetscErrorCode Initial(Vec,void*); 24 extern PetscErrorCode MyMatMult(Mat,Vec,Vec); 25 26 extern PetscReal solx(PetscReal); 27 extern PetscReal soly(PetscReal); 28 extern PetscReal solz(PetscReal); 29 30 int main(int argc,char **argv) 31 { 32 PetscErrorCode ierr; 33 PetscInt time_steps = 100,steps; 34 Vec global; 35 PetscReal dt,ftime; 36 TS ts; 37 Mat A = 0,S; 38 39 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 40 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL)); 41 42 /* set initial conditions */ 43 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&global)); 44 CHKERRQ(VecSetSizes(global,PETSC_DECIDE,3)); 45 CHKERRQ(VecSetFromOptions(global)); 46 CHKERRQ(Initial(global,NULL)); 47 48 /* make timestep context */ 49 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 50 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 51 CHKERRQ(TSMonitorSet(ts,Monitor,NULL,NULL)); 52 dt = 0.001; 53 54 /* 55 The user provides the RHS and Jacobian 56 */ 57 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunction,NULL)); 58 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 59 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3)); 60 CHKERRQ(MatSetFromOptions(A)); 61 CHKERRQ(MatSetUp(A)); 62 CHKERRQ(RHSJacobian(ts,0.0,global,A,A,NULL)); 63 CHKERRQ(TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL)); 64 65 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S)); 66 CHKERRQ(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult)); 67 CHKERRQ(TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL)); 68 69 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 70 CHKERRQ(TSSetFromOptions(ts)); 71 72 CHKERRQ(TSSetTimeStep(ts,dt)); 73 CHKERRQ(TSSetMaxSteps(ts,time_steps)); 74 CHKERRQ(TSSetMaxTime(ts,1)); 75 CHKERRQ(TSSetSolution(ts,global)); 76 77 CHKERRQ(TSSolve(ts,global)); 78 CHKERRQ(TSGetSolveTime(ts,&ftime)); 79 CHKERRQ(TSGetStepNumber(ts,&steps)); 80 81 /* free the memories */ 82 83 CHKERRQ(TSDestroy(&ts)); 84 CHKERRQ(VecDestroy(&global)); 85 CHKERRQ(MatDestroy(&A)); 86 CHKERRQ(MatDestroy(&S)); 87 88 ierr = PetscFinalize(); 89 return ierr; 90 } 91 92 PetscErrorCode MyMatMult(Mat S,Vec x,Vec y) 93 { 94 const PetscScalar *inptr; 95 PetscScalar *outptr; 96 97 PetscFunctionBegin; 98 CHKERRQ(VecGetArrayRead(x,&inptr)); 99 CHKERRQ(VecGetArrayWrite(y,&outptr)); 100 101 outptr[0] = 2.0*inptr[0]+inptr[1]; 102 outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2]; 103 outptr[2] = inptr[1]+2.0*inptr[2]; 104 105 CHKERRQ(VecRestoreArrayRead(x,&inptr)); 106 CHKERRQ(VecRestoreArrayWrite(y,&outptr)); 107 PetscFunctionReturn(0); 108 } 109 110 /* this test problem has initial values (1,1,1). */ 111 PetscErrorCode Initial(Vec global,void *ctx) 112 { 113 PetscScalar *localptr; 114 PetscInt i,mybase,myend,locsize; 115 116 /* determine starting point of each processor */ 117 CHKERRQ(VecGetOwnershipRange(global,&mybase,&myend)); 118 CHKERRQ(VecGetLocalSize(global,&locsize)); 119 120 /* Initialize the array */ 121 CHKERRQ(VecGetArrayWrite(global,&localptr)); 122 for (i=0; i<locsize; i++) localptr[i] = 1.0; 123 124 if (mybase == 0) localptr[0]=1.0; 125 126 CHKERRQ(VecRestoreArrayWrite(global,&localptr)); 127 return 0; 128 } 129 130 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx) 131 { 132 VecScatter scatter; 133 IS from,to; 134 PetscInt i,n,*idx; 135 Vec tmp_vec; 136 PetscErrorCode ierr; 137 const PetscScalar *tmp; 138 139 /* Get the size of the vector */ 140 CHKERRQ(VecGetSize(global,&n)); 141 142 /* Set the index sets */ 143 CHKERRQ(PetscMalloc1(n,&idx)); 144 for (i=0; i<n; i++) idx[i]=i; 145 146 /* Create local sequential vectors */ 147 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec)); 148 149 /* Create scatter context */ 150 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 151 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 152 CHKERRQ(VecScatterCreate(global,from,tmp_vec,to,&scatter)); 153 CHKERRQ(VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 154 CHKERRQ(VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 155 156 CHKERRQ(VecGetArrayRead(tmp_vec,&tmp)); 157 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n", 158 (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));CHKERRQ(ierr); 159 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n", 160 (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));CHKERRQ(ierr); 161 CHKERRQ(VecRestoreArrayRead(tmp_vec,&tmp)); 162 CHKERRQ(VecScatterDestroy(&scatter)); 163 CHKERRQ(ISDestroy(&from)); 164 CHKERRQ(ISDestroy(&to)); 165 CHKERRQ(PetscFree(idx)); 166 CHKERRQ(VecDestroy(&tmp_vec)); 167 return 0; 168 } 169 170 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 171 { 172 PetscScalar *outptr; 173 const PetscScalar *inptr; 174 PetscInt i,n,*idx; 175 IS from,to; 176 VecScatter scatter; 177 Vec tmp_in,tmp_out; 178 179 /* Get the length of parallel vector */ 180 CHKERRQ(VecGetSize(globalin,&n)); 181 182 /* Set the index sets */ 183 CHKERRQ(PetscMalloc1(n,&idx)); 184 for (i=0; i<n; i++) idx[i]=i; 185 186 /* Create local sequential vectors */ 187 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in)); 188 CHKERRQ(VecDuplicate(tmp_in,&tmp_out)); 189 190 /* Create scatter context */ 191 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 192 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 193 CHKERRQ(VecScatterCreate(globalin,from,tmp_in,to,&scatter)); 194 CHKERRQ(VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 195 CHKERRQ(VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 196 CHKERRQ(VecScatterDestroy(&scatter)); 197 198 /*Extract income array */ 199 CHKERRQ(VecGetArrayRead(tmp_in,&inptr)); 200 201 /* Extract outcome array*/ 202 CHKERRQ(VecGetArrayWrite(tmp_out,&outptr)); 203 204 outptr[0] = 2.0*inptr[0]+inptr[1]; 205 outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2]; 206 outptr[2] = inptr[1]+2.0*inptr[2]; 207 208 CHKERRQ(VecRestoreArrayRead(tmp_in,&inptr)); 209 CHKERRQ(VecRestoreArrayWrite(tmp_out,&outptr)); 210 211 CHKERRQ(VecScatterCreate(tmp_out,from,globalout,to,&scatter)); 212 CHKERRQ(VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 213 CHKERRQ(VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 214 215 /* Destroy idx aand scatter */ 216 CHKERRQ(ISDestroy(&from)); 217 CHKERRQ(ISDestroy(&to)); 218 CHKERRQ(VecScatterDestroy(&scatter)); 219 CHKERRQ(VecDestroy(&tmp_in)); 220 CHKERRQ(VecDestroy(&tmp_out)); 221 CHKERRQ(PetscFree(idx)); 222 return 0; 223 } 224 225 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx) 226 { 227 PetscScalar v[3]; 228 const PetscScalar *tmp; 229 PetscInt idx[3],i; 230 231 idx[0]=0; idx[1]=1; idx[2]=2; 232 CHKERRQ(VecGetArrayRead(x,&tmp)); 233 234 i = 0; 235 v[0] = 2.0; v[1] = 1.0; v[2] = 0.0; 236 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 237 238 i = 1; 239 v[0] = 1.0; v[1] = 2.0; v[2] = 1.0; 240 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 241 242 i = 2; 243 v[0] = 0.0; v[1] = 1.0; v[2] = 2.0; 244 CHKERRQ(MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES)); 245 246 CHKERRQ(MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY)); 247 CHKERRQ(MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY)); 248 249 if (A != BB) { 250 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 251 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 252 } 253 CHKERRQ(VecRestoreArrayRead(x,&tmp)); 254 255 return 0; 256 } 257 258 /* 259 The exact solutions 260 */ 261 PetscReal solx(PetscReal t) 262 { 263 return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) + 264 PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)); 265 } 266 267 PetscReal soly(PetscReal t) 268 { 269 return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) + 270 PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0); 271 } 272 273 PetscReal solz(PetscReal t) 274 { 275 return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) + 276 PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)); 277 } 278 279 /*TEST 280 281 test: 282 suffix: euler 283 args: -ts_type euler 284 requires: !single 285 286 test: 287 suffix: beuler 288 args: -ts_type beuler 289 requires: !single 290 291 TEST*/ 292