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) PetscCall(TSSetPostStep(ts,PostStep)); 154 155 PetscCall(PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL)); 156 for (iout=1; iout<=NOUT; iout++) { 157 PetscCall(TSSetMaxSteps(ts,time_steps)); 158 PetscCall(TSSetMaxTime(ts,iout*ftime_original/NOUT)); 159 PetscCall(TSSolve(ts,global)); 160 PetscCall(TSGetSolveTime(ts,&ftime)); 161 PetscCall(TSSetTime(ts,ftime)); 162 PetscCall(TSSetTimeStep(ts,dt)); 163 } 164 /* Interpolate solution at tfinal */ 165 PetscCall(TSGetSolution(ts,&global)); 166 PetscCall(TSInterpolate(ts,ftime_original,global)); 167 168 PetscCall(PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg)); 169 if (flg) { /* print solution into a MATLAB file */ 170 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile)); 171 PetscCall(PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB)); 172 PetscCall(VecView(global,viewfile)); 173 PetscCall(PetscViewerPopFormat(viewfile)); 174 PetscCall(PetscViewerDestroy(&viewfile)); 175 } 176 177 /* free the memories */ 178 PetscCall(TSDestroy(&ts)); 179 PetscCall(VecDestroy(&global)); 180 PetscCall(VecDestroy(&x)); 181 PetscCall(MatDestroy(&J)); 182 if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 183 PetscCall(PetscFinalize()); 184 return 0; 185 } 186 187 /* -------------------------------------------------------------------*/ 188 /* the initial function */ 189 PetscReal f_ini(PetscReal x,PetscReal y) 190 { 191 PetscReal f; 192 193 f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2))); 194 return f; 195 } 196 197 PetscErrorCode Initial(Vec global,void *ctx) 198 { 199 Data *data = (Data*)ctx; 200 PetscInt m,row,col; 201 PetscReal x,y,dx,dy; 202 PetscScalar *localptr; 203 PetscInt i,mybase,myend,locsize; 204 205 PetscFunctionBeginUser; 206 /* make the local copies of parameters */ 207 m = data->m; 208 dx = data->dx; 209 dy = data->dy; 210 211 /* determine starting point of each processor */ 212 PetscCall(VecGetOwnershipRange(global,&mybase,&myend)); 213 PetscCall(VecGetLocalSize(global,&locsize)); 214 215 /* Initialize the array */ 216 PetscCall(VecGetArrayWrite(global,&localptr)); 217 218 for (i=0; i<locsize; i++) { 219 row = 1+(mybase+i)-((mybase+i)/m)*m; 220 col = (mybase+i)/m+1; 221 x = dx*row; 222 y = dy*col; 223 localptr[i] = f_ini(x,y); 224 } 225 226 PetscCall(VecRestoreArrayWrite(global,&localptr)); 227 PetscFunctionReturn(0); 228 } 229 230 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx) 231 { 232 VecScatter scatter; 233 IS from,to; 234 PetscInt i,n,*idx,nsteps,maxsteps; 235 Vec tmp_vec; 236 const PetscScalar *tmp; 237 238 PetscFunctionBeginUser; 239 PetscCall(TSGetStepNumber(ts,&nsteps)); 240 /* display output at selected time steps */ 241 PetscCall(TSGetMaxSteps(ts, &maxsteps)); 242 if (nsteps % 10 != 0) PetscFunctionReturn(0); 243 244 /* Get the size of the vector */ 245 PetscCall(VecGetSize(global,&n)); 246 247 /* Set the index sets */ 248 PetscCall(PetscMalloc1(n,&idx)); 249 for (i=0; i<n; i++) idx[i]=i; 250 251 /* Create local sequential vectors */ 252 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec)); 253 254 /* Create scatter context */ 255 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from)); 256 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to)); 257 PetscCall(VecScatterCreate(global,from,tmp_vec,to,&scatter)); 258 PetscCall(VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 259 PetscCall(VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD)); 260 261 PetscCall(VecGetArrayRead(tmp_vec,&tmp)); 262 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]))); 263 PetscCall(VecRestoreArrayRead(tmp_vec,&tmp)); 264 265 PetscCall(PetscFree(idx)); 266 PetscCall(ISDestroy(&from)); 267 PetscCall(ISDestroy(&to)); 268 PetscCall(VecScatterDestroy(&scatter)); 269 PetscCall(VecDestroy(&tmp_vec)); 270 PetscFunctionReturn(0); 271 } 272 273 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr) 274 { 275 Data *data = (Data*)ptr; 276 PetscScalar v[5]; 277 PetscInt idx[5],i,j,row; 278 PetscInt m,n,mn; 279 PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr; 280 281 PetscFunctionBeginUser; 282 m = data->m; 283 n = data->n; 284 mn = m*n; 285 dx = data->dx; 286 dy = data->dy; 287 a = data->a; 288 epsilon = data->epsilon; 289 290 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy)); 291 xl = 0.5*a/dx+epsilon/(dx*dx); 292 xr = -0.5*a/dx+epsilon/(dx*dx); 293 yl = 0.5*a/dy+epsilon/(dy*dy); 294 yr = -0.5*a/dy+epsilon/(dy*dy); 295 296 row = 0; 297 v[0] = xc; v[1] = xr; v[2] = yr; 298 idx[0] = 0; idx[1] = 2; idx[2] = m; 299 PetscCall(MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES)); 300 301 row = m-1; 302 v[0] = 2.0*xl; v[1] = xc; v[2] = yr; 303 idx[0] = m-2; idx[1] = m-1; idx[2] = m-1+m; 304 PetscCall(MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES)); 305 306 for (i=1; i<m-1; i++) { 307 row = i; 308 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = yr; 309 idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m; 310 PetscCall(MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES)); 311 } 312 313 for (j=1; j<n-1; j++) { 314 row = j*m; 315 v[0] = xc; v[1] = xr; v[2] = yl; v[3] = yr; 316 idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m; 317 PetscCall(MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES)); 318 319 row = j*m+m-1; 320 v[0] = xc; v[1] = 2.0*xl; v[2] = yl; v[3] = yr; 321 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; 322 PetscCall(MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES)); 323 324 for (i=1; i<m-1; i++) { 325 row = j*m+i; 326 v[0] = xc; v[1] = xl; v[2] = xr; v[3] = yl; v[4]=yr; 327 idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m; 328 idx[4] = j*m+i+m; 329 PetscCall(MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES)); 330 } 331 } 332 333 row = mn-m; 334 v[0] = xc; v[1] = xr; v[2] = 2.0*yl; 335 idx[0] = mn-m; idx[1] = mn-m+1; idx[2] = mn-m-m; 336 PetscCall(MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES)); 337 338 row = mn-1; 339 v[0] = xc; v[1] = 2.0*xl; v[2] = 2.0*yl; 340 idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m; 341 PetscCall(MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES)); 342 343 for (i=1; i<m-1; i++) { 344 row = mn-m+i; 345 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = 2.0*yl; 346 idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m; 347 PetscCall(MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES)); 348 } 349 350 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 351 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 352 353 PetscFunctionReturn(0); 354 } 355 356 /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */ 357 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 358 { 359 Data *data = (Data*)ctx; 360 PetscInt m,n,mn; 361 PetscReal dx,dy; 362 PetscReal xc,xl,xr,yl,yr; 363 PetscReal a,epsilon; 364 PetscScalar *outptr; 365 const PetscScalar *inptr; 366 PetscInt i,j,len; 367 IS from,to; 368 PetscInt *idx; 369 VecScatter scatter; 370 Vec tmp_in,tmp_out; 371 372 PetscFunctionBeginUser; 373 m = data->m; 374 n = data->n; 375 mn = m*n; 376 dx = data->dx; 377 dy = data->dy; 378 a = data->a; 379 epsilon = data->epsilon; 380 381 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy)); 382 xl = 0.5*a/dx+epsilon/(dx*dx); 383 xr = -0.5*a/dx+epsilon/(dx*dx); 384 yl = 0.5*a/dy+epsilon/(dy*dy); 385 yr = -0.5*a/dy+epsilon/(dy*dy); 386 387 /* Get the length of parallel vector */ 388 PetscCall(VecGetSize(globalin,&len)); 389 390 /* Set the index sets */ 391 PetscCall(PetscMalloc1(len,&idx)); 392 for (i=0; i<len; i++) idx[i]=i; 393 394 /* Create local sequential vectors */ 395 PetscCall(VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in)); 396 PetscCall(VecDuplicate(tmp_in,&tmp_out)); 397 398 /* Create scatter context */ 399 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from)); 400 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to)); 401 PetscCall(VecScatterCreate(globalin,from,tmp_in,to,&scatter)); 402 PetscCall(VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 403 PetscCall(VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD)); 404 PetscCall(VecScatterDestroy(&scatter)); 405 406 /*Extract income array - include ghost points */ 407 PetscCall(VecGetArrayRead(tmp_in,&inptr)); 408 409 /* Extract outcome array*/ 410 PetscCall(VecGetArrayWrite(tmp_out,&outptr)); 411 412 outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m]; 413 outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m]; 414 for (i=1; i<m-1; i++) { 415 outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m]; 416 } 417 418 for (j=1; j<n-1; j++) { 419 outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m]; 420 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]; 421 for (i=1; i<m-1; i++) { 422 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]; 423 } 424 } 425 426 outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m]; 427 outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m]; 428 for (i=1; i<m-1; i++) { 429 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]; 430 } 431 432 PetscCall(VecRestoreArrayRead(tmp_in,&inptr)); 433 PetscCall(VecRestoreArrayWrite(tmp_out,&outptr)); 434 435 PetscCall(VecScatterCreate(tmp_out,from,globalout,to,&scatter)); 436 PetscCall(VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 437 PetscCall(VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD)); 438 439 /* Destroy idx aand scatter */ 440 PetscCall(VecDestroy(&tmp_in)); 441 PetscCall(VecDestroy(&tmp_out)); 442 PetscCall(ISDestroy(&from)); 443 PetscCall(ISDestroy(&to)); 444 PetscCall(VecScatterDestroy(&scatter)); 445 446 PetscCall(PetscFree(idx)); 447 PetscFunctionReturn(0); 448 } 449 450 PetscErrorCode PostStep(TS ts) 451 { 452 PetscReal t; 453 454 PetscFunctionBeginUser; 455 PetscCall(TSGetTime(ts,&t)); 456 PetscCall(PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t)); 457 PetscFunctionReturn(0); 458 } 459 460 /*TEST 461 462 test: 463 args: -ts_fd -ts_type beuler 464 output_file: output/ex4.out 465 466 test: 467 suffix: 2 468 args: -ts_fd -ts_type beuler 469 nsize: 2 470 output_file: output/ex4.out 471 472 test: 473 suffix: 3 474 args: -ts_fd -ts_type cn 475 476 test: 477 suffix: 4 478 args: -ts_fd -ts_type cn 479 output_file: output/ex4_3.out 480 nsize: 2 481 482 test: 483 suffix: 5 484 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl 485 output_file: output/ex4.out 486 487 test: 488 suffix: 6 489 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl 490 output_file: output/ex4.out 491 nsize: 2 492 493 test: 494 suffix: 7 495 requires: !single 496 args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1 497 498 test: 499 suffix: 8 500 requires: !single 501 args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view 502 503 TEST*/ 504