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