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