1 /* 2 The Problem: 3 Solve the convection-diffusion equation: 4 5 u_t+a*(u_x+u_y)=epsilon*(u_xx+u_yy) 6 u=0 at x=0, y=0 7 u_x=0 at x=1 8 u_y=0 at y=1 9 u = exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0))) at t=0 10 11 This program tests the routine of computing the Jacobian by the 12 finite difference method as well as PETSc. 13 14 */ 15 16 static char help[] = "Solve the convection-diffusion equation. \n\n"; 17 18 #include <petscts.h> 19 20 typedef struct 21 { 22 PetscInt m; /* the number of mesh points in x-direction */ 23 PetscInt n; /* the number of mesh points in y-direction */ 24 PetscReal dx; /* the grid space in x-direction */ 25 PetscReal dy; /* the grid space in y-direction */ 26 PetscReal a; /* the convection coefficient */ 27 PetscReal epsilon; /* the diffusion coefficient */ 28 PetscReal tfinal; 29 } Data; 30 31 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); 32 extern PetscErrorCode Initial(Vec,void*); 33 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 34 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 35 extern PetscErrorCode PostStep(TS); 36 37 int main(int argc,char **argv) 38 { 39 PetscErrorCode ierr; 40 PetscInt time_steps=100,iout,NOUT=1; 41 Vec global; 42 PetscReal dt,ftime,ftime_original; 43 TS ts; 44 PetscViewer viewfile; 45 Mat J = 0; 46 Vec x; 47 Data data; 48 PetscInt mn; 49 PetscBool flg; 50 MatColoring mc; 51 ISColoring iscoloring; 52 MatFDColoring matfdcoloring = 0; 53 PetscBool fd_jacobian_coloring = PETSC_FALSE; 54 SNES snes; 55 KSP ksp; 56 PC pc; 57 58 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 59 60 /* set data */ 61 data.m = 9; 62 data.n = 9; 63 data.a = 1.0; 64 data.epsilon = 0.1; 65 data.dx = 1.0/(data.m+1.0); 66 data.dy = 1.0/(data.n+1.0); 67 mn = (data.m)*(data.n); 68 ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr); 69 70 /* set initial conditions */ 71 ierr = VecCreate(PETSC_COMM_WORLD,&global);CHKERRQ(ierr); 72 ierr = VecSetSizes(global,PETSC_DECIDE,mn);CHKERRQ(ierr); 73 ierr = VecSetFromOptions(global);CHKERRQ(ierr); 74 ierr = Initial(global,&data);CHKERRQ(ierr); 75 ierr = VecDuplicate(global,&x);CHKERRQ(ierr); 76 77 /* create timestep context */ 78 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 79 ierr = TSMonitorSet(ts,Monitor,&data,NULL);CHKERRQ(ierr); 80 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 81 dt = 0.1; 82 ftime_original = data.tfinal = 1.0; 83 84 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 85 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr); 86 ierr = TSSetMaxTime(ts,ftime_original);CHKERRQ(ierr); 87 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 88 ierr = TSSetSolution(ts,global);CHKERRQ(ierr); 89 90 /* set user provided RHSFunction and RHSJacobian */ 91 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&data);CHKERRQ(ierr); 92 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 93 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);CHKERRQ(ierr); 94 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 95 ierr = MatSeqAIJSetPreallocation(J,5,NULL);CHKERRQ(ierr); 96 ierr = MatMPIAIJSetPreallocation(J,5,NULL,5,NULL);CHKERRQ(ierr); 97 98 ierr = PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);CHKERRQ(ierr); 99 if (!flg) { 100 ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);CHKERRQ(ierr); 101 } else { 102 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 103 ierr = PetscOptionsHasName(NULL,NULL,"-fd_color",&fd_jacobian_coloring);CHKERRQ(ierr); 104 if (fd_jacobian_coloring) { /* Use finite differences with coloring */ 105 /* Get data structure of J */ 106 PetscBool pc_diagonal; 107 ierr = PetscOptionsHasName(NULL,NULL,"-pc_diagonal",&pc_diagonal);CHKERRQ(ierr); 108 if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */ 109 PetscInt rstart,rend,i; 110 PetscScalar zero=0.0; 111 ierr = MatGetOwnershipRange(J,&rstart,&rend);CHKERRQ(ierr); 112 for (i=rstart; i<rend; i++) { 113 ierr = MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);CHKERRQ(ierr); 114 } 115 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 117 } else { 118 /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */ 119 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); 120 ierr = TSSetUp(ts);CHKERRQ(ierr); 121 ierr = SNESComputeJacobianDefault(snes,x,J,J,ts);CHKERRQ(ierr); 122 } 123 124 /* create coloring context */ 125 ierr = MatColoringCreate(J,&mc);CHKERRQ(ierr); 126 ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr); 127 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 128 ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr); 129 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 130 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); 131 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);CHKERRQ(ierr); 132 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 133 ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr); 134 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr); 135 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 136 } else { /* Use finite differences (slow) */ 137 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr); 138 } 139 } 140 141 /* Pick up a Petsc preconditioner */ 142 /* one can always set method or preconditioner during the run time */ 143 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 144 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 145 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 146 ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr); 147 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 148 149 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 150 ierr = TSSetUp(ts);CHKERRQ(ierr); 151 152 /* Test TSSetPostStep() */ 153 ierr = PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);CHKERRQ(ierr); 154 if (flg) { 155 ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr); 156 } 157 158 ierr = PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);CHKERRQ(ierr); 159 for (iout=1; iout<=NOUT; iout++) { 160 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr); 161 ierr = TSSetMaxTime(ts,iout*ftime_original/NOUT);CHKERRQ(ierr); 162 ierr = TSSolve(ts,global);CHKERRQ(ierr); 163 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 164 ierr = TSSetTime(ts,ftime);CHKERRQ(ierr); 165 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 166 } 167 /* Interpolate solution at tfinal */ 168 ierr = TSGetSolution(ts,&global);CHKERRQ(ierr); 169 ierr = TSInterpolate(ts,ftime_original,global);CHKERRQ(ierr); 170 171 ierr = PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);CHKERRQ(ierr); 172 if (flg) { /* print solution into a MATLAB file */ 173 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);CHKERRQ(ierr); 174 ierr = PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 175 ierr = VecView(global,viewfile);CHKERRQ(ierr); 176 ierr = PetscViewerPopFormat(viewfile);CHKERRQ(ierr); 177 ierr = PetscViewerDestroy(&viewfile);CHKERRQ(ierr); 178 } 179 180 /* free the memories */ 181 ierr = TSDestroy(&ts);CHKERRQ(ierr); 182 ierr = VecDestroy(&global);CHKERRQ(ierr); 183 ierr = VecDestroy(&x);CHKERRQ(ierr); 184 ierr = MatDestroy(&J);CHKERRQ(ierr); 185 if (fd_jacobian_coloring) {ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);} 186 ierr = PetscFinalize(); 187 return ierr; 188 } 189 190 /* -------------------------------------------------------------------*/ 191 /* the initial function */ 192 PetscReal f_ini(PetscReal x,PetscReal y) 193 { 194 PetscReal f; 195 196 f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2))); 197 return f; 198 } 199 200 PetscErrorCode Initial(Vec global,void *ctx) 201 { 202 Data *data = (Data*)ctx; 203 PetscInt m,row,col; 204 PetscReal x,y,dx,dy; 205 PetscScalar *localptr; 206 PetscInt i,mybase,myend,locsize; 207 PetscErrorCode ierr; 208 209 PetscFunctionBeginUser; 210 /* make the local copies of parameters */ 211 m = data->m; 212 dx = data->dx; 213 dy = data->dy; 214 215 /* determine starting point of each processor */ 216 ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr); 217 ierr = VecGetLocalSize(global,&locsize);CHKERRQ(ierr); 218 219 /* Initialize the array */ 220 ierr = VecGetArrayWrite(global,&localptr);CHKERRQ(ierr); 221 222 for (i=0; i<locsize; i++) { 223 row = 1+(mybase+i)-((mybase+i)/m)*m; 224 col = (mybase+i)/m+1; 225 x = dx*row; 226 y = dy*col; 227 localptr[i] = f_ini(x,y); 228 } 229 230 ierr = VecRestoreArrayWrite(global,&localptr);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx) 235 { 236 VecScatter scatter; 237 IS from,to; 238 PetscInt i,n,*idx,nsteps,maxsteps; 239 Vec tmp_vec; 240 PetscErrorCode ierr; 241 const PetscScalar *tmp; 242 243 PetscFunctionBeginUser; 244 ierr = TSGetStepNumber(ts,&nsteps);CHKERRQ(ierr); 245 /* display output at selected time steps */ 246 ierr = TSGetMaxSteps(ts, &maxsteps);CHKERRQ(ierr); 247 if (nsteps % 10 != 0) PetscFunctionReturn(0); 248 249 /* Get the size of the vector */ 250 ierr = VecGetSize(global,&n);CHKERRQ(ierr); 251 252 /* Set the index sets */ 253 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 254 for (i=0; i<n; i++) idx[i]=i; 255 256 /* Create local sequential vectors */ 257 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);CHKERRQ(ierr); 258 259 /* Create scatter context */ 260 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 261 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 262 ierr = VecScatterCreate(global,from,tmp_vec,to,&scatter);CHKERRQ(ierr); 263 ierr = VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 264 ierr = VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 265 266 ierr = VecGetArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 267 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));CHKERRQ(ierr); 268 ierr = VecRestoreArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 269 270 ierr = PetscFree(idx);CHKERRQ(ierr); 271 ierr = ISDestroy(&from);CHKERRQ(ierr); 272 ierr = ISDestroy(&to);CHKERRQ(ierr); 273 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 274 ierr = VecDestroy(&tmp_vec);CHKERRQ(ierr); 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr) 279 { 280 Data *data = (Data*)ptr; 281 PetscScalar v[5]; 282 PetscInt idx[5],i,j,row; 283 PetscErrorCode ierr; 284 PetscInt m,n,mn; 285 PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr; 286 287 PetscFunctionBeginUser; 288 m = data->m; 289 n = data->n; 290 mn = m*n; 291 dx = data->dx; 292 dy = data->dy; 293 a = data->a; 294 epsilon = data->epsilon; 295 296 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy)); 297 xl = 0.5*a/dx+epsilon/(dx*dx); 298 xr = -0.5*a/dx+epsilon/(dx*dx); 299 yl = 0.5*a/dy+epsilon/(dy*dy); 300 yr = -0.5*a/dy+epsilon/(dy*dy); 301 302 row = 0; 303 v[0] = xc; v[1] = xr; v[2] = yr; 304 idx[0] = 0; idx[1] = 2; idx[2] = m; 305 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 306 307 row = m-1; 308 v[0] = 2.0*xl; v[1] = xc; v[2] = yr; 309 idx[0] = m-2; idx[1] = m-1; idx[2] = m-1+m; 310 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 311 312 for (i=1; i<m-1; i++) { 313 row = i; 314 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = yr; 315 idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m; 316 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr); 317 } 318 319 for (j=1; j<n-1; j++) { 320 row = j*m; 321 v[0] = xc; v[1] = xr; v[2] = yl; v[3] = yr; 322 idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m; 323 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr); 324 325 row = j*m+m-1; 326 v[0] = xc; v[1] = 2.0*xl; v[2] = yl; v[3] = yr; 327 idx[0] = j*m+m-1; idx[1] = j*m+m-1-1; idx[2] = j*m+m-1-m; idx[3] = j*m+m-1+m; 328 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr); 329 330 for (i=1; i<m-1; i++) { 331 row = j*m+i; 332 v[0] = xc; v[1] = xl; v[2] = xr; v[3] = yl; v[4]=yr; 333 idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m; 334 idx[4] = j*m+i+m; 335 ierr = MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);CHKERRQ(ierr); 336 } 337 } 338 339 row = mn-m; 340 v[0] = xc; v[1] = xr; v[2] = 2.0*yl; 341 idx[0] = mn-m; idx[1] = mn-m+1; idx[2] = mn-m-m; 342 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 343 344 row = mn-1; 345 v[0] = xc; v[1] = 2.0*xl; v[2] = 2.0*yl; 346 idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m; 347 ierr = MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 348 349 for (i=1; i<m-1; i++) { 350 row = mn-m+i; 351 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = 2.0*yl; 352 idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m; 353 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr); 354 } 355 356 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 357 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 358 359 PetscFunctionReturn(0); 360 } 361 362 /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */ 363 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 364 { 365 Data *data = (Data*)ctx; 366 PetscInt m,n,mn; 367 PetscReal dx,dy; 368 PetscReal xc,xl,xr,yl,yr; 369 PetscReal a,epsilon; 370 PetscScalar *outptr; 371 const PetscScalar *inptr; 372 PetscInt i,j,len; 373 PetscErrorCode ierr; 374 IS from,to; 375 PetscInt *idx; 376 VecScatter scatter; 377 Vec tmp_in,tmp_out; 378 379 PetscFunctionBeginUser; 380 m = data->m; 381 n = data->n; 382 mn = m*n; 383 dx = data->dx; 384 dy = data->dy; 385 a = data->a; 386 epsilon = data->epsilon; 387 388 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy)); 389 xl = 0.5*a/dx+epsilon/(dx*dx); 390 xr = -0.5*a/dx+epsilon/(dx*dx); 391 yl = 0.5*a/dy+epsilon/(dy*dy); 392 yr = -0.5*a/dy+epsilon/(dy*dy); 393 394 /* Get the length of parallel vector */ 395 ierr = VecGetSize(globalin,&len);CHKERRQ(ierr); 396 397 /* Set the index sets */ 398 ierr = PetscMalloc1(len,&idx);CHKERRQ(ierr); 399 for (i=0; i<len; i++) idx[i]=i; 400 401 /* Create local sequential vectors */ 402 ierr = VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);CHKERRQ(ierr); 403 ierr = VecDuplicate(tmp_in,&tmp_out);CHKERRQ(ierr); 404 405 /* Create scatter context */ 406 ierr = ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 407 ierr = ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 408 ierr = VecScatterCreate(globalin,from,tmp_in,to,&scatter);CHKERRQ(ierr); 409 ierr = VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 410 ierr = VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 411 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 412 413 /*Extract income array - include ghost points */ 414 ierr = VecGetArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 415 416 /* Extract outcome array*/ 417 ierr = VecGetArrayWrite(tmp_out,&outptr);CHKERRQ(ierr); 418 419 outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m]; 420 outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m]; 421 for (i=1; i<m-1; i++) { 422 outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m]; 423 } 424 425 for (j=1; j<n-1; j++) { 426 outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m]; 427 outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+ yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m]; 428 for (i=1; i<m-1; i++) { 429 outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]+yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m]; 430 } 431 } 432 433 outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m]; 434 outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m]; 435 for (i=1; i<m-1; i++) { 436 outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]+2*yl*inptr[mn-m+i-m]; 437 } 438 439 ierr = VecRestoreArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 440 ierr = VecRestoreArrayWrite(tmp_out,&outptr);CHKERRQ(ierr); 441 442 ierr = VecScatterCreate(tmp_out,from,globalout,to,&scatter);CHKERRQ(ierr); 443 ierr = VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 444 ierr = VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 445 446 /* Destroy idx aand scatter */ 447 ierr = VecDestroy(&tmp_in);CHKERRQ(ierr); 448 ierr = VecDestroy(&tmp_out);CHKERRQ(ierr); 449 ierr = ISDestroy(&from);CHKERRQ(ierr); 450 ierr = ISDestroy(&to);CHKERRQ(ierr); 451 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 452 453 ierr = PetscFree(idx);CHKERRQ(ierr); 454 PetscFunctionReturn(0); 455 } 456 457 PetscErrorCode PostStep(TS ts) 458 { 459 PetscErrorCode ierr; 460 PetscReal t; 461 462 PetscFunctionBeginUser; 463 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 464 ierr = PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t);CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 /*TEST 469 470 test: 471 args: -ts_fd -ts_type beuler 472 output_file: output/ex4.out 473 474 test: 475 suffix: 2 476 args: -ts_fd -ts_type beuler 477 nsize: 2 478 output_file: output/ex4.out 479 480 test: 481 suffix: 3 482 args: -ts_fd -ts_type cn 483 484 test: 485 suffix: 4 486 args: -ts_fd -ts_type cn 487 output_file: output/ex4_3.out 488 nsize: 2 489 490 test: 491 suffix: 5 492 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl 493 output_file: output/ex4.out 494 495 test: 496 suffix: 6 497 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl 498 output_file: output/ex4.out 499 nsize: 2 500 501 test: 502 suffix: 7 503 requires: !single 504 args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1 505 506 test: 507 suffix: 8 508 requires: !single 509 args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view 510 511 TEST*/ 512 513