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