1! Time-dependent advection-reaction PDE in 1d. Demonstrates IMEX methods 2! 3! u_t + a1*u_x = -k1*u + k2*v + s1 4! v_t + a2*v_x = k1*u - k2*v + s2 5! 0 < x < 1; 6! a1 = 1, k1 = 10^6, s1 = 0, 7! a2 = 0, k2 = 2*k1, s2 = 1 8! 9! Initial conditions: 10! u(x,0) = 1 + s2*x 11! v(x,0) = k0/k1*u(x,0) + s1/k1 12! 13! Upstream boundary conditions: 14! u(0,t) = 1-sin(12*t)^4 15! 16 17 module ex22f_mfmodule 18#include <petsc/finclude/petscts.h> 19 use petscts 20 PetscScalar::PETSC_SHIFT 21 TS::tscontext 22 Mat::Jmat 23 PetscReal::MFuser(6) 24 end module ex22f_mfmodule 25 26program main 27 use ex22f_mfmodule 28 use petscdm 29 implicit none 30 31 ! 32 ! Create an application context to contain data needed by the 33 ! application-provided call-back routines, FormJacobian() and 34 ! FormFunction(). We use a double precision array with six 35 ! entries, two for each problem parameter a, k, s. 36 ! 37 PetscReal user(6) 38 integer user_a,user_k,user_s 39 parameter (user_a = 0,user_k = 2,user_s = 4) 40 41 external FormRHSFunction,FormIFunction 42 external FormInitialSolution 43 external FormIJacobian 44 external MyMult,FormIJacobianMF 45 46 TS ts 47 Vec X 48 Mat J 49 PetscInt mx 50 PetscBool OptionSaveToDisk 51 PetscErrorCode ierr 52 DM da 53 PetscReal ftime,dt 54 PetscReal one,pone 55 PetscInt im11,i2 56 PetscBool flg 57 58 PetscInt xs,xe,gxs,gxe,dof,gdof 59 PetscScalar shell_shift 60 Mat A 61 62 im11 = 11 63 i2 = 2 64 one = 1.0 65 pone = one / 10 66 67 PetscCallA(PetscInitialize(ierr)) 68 69 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 70 ! Create distributed array (DMDA) to manage parallel grid and vectors 71 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 72 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER_ARRAY,da,ierr)) 73 PetscCallA(DMSetFromOptions(da,ierr)) 74 PetscCallA(DMSetUp(da,ierr)) 75 76 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 ! Extract global vectors from DMDA; 78 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 PetscCallA(DMCreateGlobalVector(da,X,ierr)) 80 81 ! Initialize user application context 82 ! Use zero-based indexing for command line parameters to match ex22.c 83 user(user_a+1) = 1.0 84 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr)) 85 user(user_a+2) = 0.0 86 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr)) 87 user(user_k+1) = 1000000.0 88 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0', user(user_k+1),flg,ierr)) 89 user(user_k+2) = 2*user(user_k+1) 90 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr)) 91 user(user_s+1) = 0.0 92 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr)) 93 user(user_s+2) = 1.0 94 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr)) 95 96 OptionSaveToDisk=.FALSE. 97 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sdisk',OptionSaveToDisk,flg,ierr)) 98 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 ! Create timestepping solver context 100 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr)) 102 tscontext=ts 103 PetscCallA(TSSetDM(ts,da,ierr)) 104 PetscCallA(TSSetType(ts,TSARKIMEX,ierr)) 105 PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr)) 106 107 ! - - - - - - - - -- - - - - 108 ! Matrix free setup 109 PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 110 dof=i2*(xe-xs+1) 111 gdof=i2*(gxe-gxs+1) 112 PetscCallA(MatCreateShell(PETSC_COMM_WORLD,dof,dof,gdof,gdof,shell_shift,A,ierr)) 113 114 PetscCallA(MatShellSetOperation(A,MATOP_MULT,MyMult,ierr)) 115 ! - - - - - - - - - - - - 116 117 PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr)) 118 PetscCallA(DMSetMatType(da,MATAIJ,ierr)) 119 PetscCallA(DMCreateMatrix(da,J,ierr)) 120 121 Jmat=J 122 123 PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)) 124 PetscCallA(TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr)) 125 126 ftime = 1.0 127 PetscCallA(TSSetMaxTime(ts,ftime,ierr)) 128 PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)) 129 130 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 ! Set initial conditions 132 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 PetscCallA(FormInitialSolution(ts,X,user,ierr)) 134 PetscCallA(TSSetSolution(ts,X,ierr)) 135 PetscCallA(VecGetSize(X,mx,ierr)) 136 ! Advective CFL, I don't know why it needs so much safety factor. 137 dt = pone * max(user(user_a+1),user(user_a+2)) / mx; 138 PetscCallA(TSSetTimeStep(ts,dt,ierr)) 139 140 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 ! Set runtime options 142 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 PetscCallA(TSSetFromOptions(ts,ierr)) 144 145 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 ! Solve nonlinear system 147 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 148 PetscCallA(TSSolve(ts,X,ierr)) 149 150 if (OptionSaveToDisk) then 151 PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 152 dof=i2*(xe-xs+1) 153 gdof=i2*(gxe-gxs+1) 154 call SaveSolutionToDisk(da,X,gdof,xs,xe) 155 end if 156 157 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 ! Free work space. 159! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 PetscCallA(MatDestroy(A,ierr)) 161 PetscCallA(MatDestroy(J,ierr)) 162 PetscCallA(VecDestroy(X,ierr)) 163 PetscCallA(TSDestroy(ts,ierr)) 164 PetscCallA(DMDestroy(da,ierr)) 165 PetscCallA(PetscFinalize(ierr)) 166end program main 167 168! Small helper to extract the layout, result uses 1-based indexing. 169 subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr) 170 use petscdm 171 implicit none 172 173 DM da 174 PetscInt mx,xs,xe,gxs,gxe 175 PetscErrorCode ierr 176 PetscInt xm,gxm 177 PetscCall(DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,ierr)) 178 PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 179 PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 180 xs = xs + 1 181 gxs = gxs + 1 182 xe = xs + xm - 1 183 gxe = gxs + gxm - 1 184end subroutine GetLayout 185 186subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr) 187 implicit none 188 PetscInt mx,xs,xe,gxs,gxe 189 PetscScalar x(2,xs:xe) 190 PetscScalar xdot(2,xs:xe) 191 PetscScalar f(2,xs:xe) 192 PetscReal a(2),k(2),s(2) 193 PetscErrorCode ierr 194 PetscInt i 195 do i = xs,xe 196 f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1) 197 f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2) 198 end do 199end subroutine FormIFunctionLocal 200 201subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr) 202 use petscdm 203 use petscts 204 implicit none 205 206 TS ts 207 PetscReal t 208 Vec X,Xdot,F 209 PetscReal user(6) 210 PetscErrorCode ierr 211 integer user_a,user_k,user_s 212 parameter (user_a = 1,user_k = 3,user_s = 5) 213 214 DM da 215 PetscInt mx,xs,xe,gxs,gxe 216 PetscScalar,pointer :: xx(:),xxdot(:),ff(:) 217 218 PetscCall(TSGetDM(ts,da,ierr)) 219 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 220 221 ! Get access to vector data 222 PetscCall(VecGetArrayRead(X,xx,ierr)) 223 PetscCall(VecGetArrayRead(Xdot,xxdot,ierr)) 224 PetscCall(VecGetArray(F,ff,ierr)) 225 226 PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx,xxdot,ff,user(user_a),user(user_k),user(user_s),ierr)) 227 228 PetscCall(VecRestoreArrayRead(X,xx,ierr)) 229 PetscCall(VecRestoreArrayRead(Xdot,xxdot,ierr)) 230 PetscCall(VecRestoreArray(F,ff,ierr)) 231end subroutine FormIFunction 232 233subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr) 234 implicit none 235 PetscInt mx,xs,xe,gxs,gxe 236 PetscReal t 237 PetscScalar x(2,gxs:gxe),f(2,xs:xe) 238 PetscReal a(2),k(2),s(2) 239 PetscErrorCode ierr 240 PetscInt i,j 241 PetscReal hx,u0t(2) 242 PetscReal one,two,three,four,six,twelve 243 PetscReal half,third,twothird,sixth 244 PetscReal twelfth 245 246 one = 1.0 247 two = 2.0 248 three = 3.0 249 four = 4.0 250 six = 6.0 251 twelve = 12.0 252 hx = one / mx 253 u0t(1) = one - sin(twelve*t)**four 254 u0t(2) = 0.0 255 half = one/two 256 third = one / three 257 twothird = two / three 258 sixth = one / six 259 twelfth = one / twelve 260 do i = xs,xe 261 do j = 1,2 262 if (i .eq. 1) then 263 f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) + sixth*x(j,i+2)) 264 else if (i .eq. 2) then 265 f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2)) 266 else if (i .eq. mx-1) then 267 f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) - half*x(j,i) -third*x(j,i+1)) 268 else if (i .eq. mx) then 269 f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1)) 270 else 271 f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2)) 272 end if 273 end do 274 end do 275 276#ifdef EXPLICIT_INTEGRATOR22 277 do i = xs,xe 278 f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1)) 279 f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2)) 280 end do 281#endif 282 283end subroutine FormRHSFunctionLocal 284 285subroutine FormRHSFunction(ts,t,X,F,user,ierr) 286 use petscts 287 use petscdm 288 implicit none 289 290 TS ts 291 PetscReal t 292 Vec X,F 293 PetscReal user(6) 294 PetscErrorCode ierr 295 integer user_a,user_k,user_s 296 parameter (user_a = 1,user_k = 3,user_s = 5) 297 DM da 298 Vec Xloc 299 PetscInt mx,xs,xe,gxs,gxe 300 PetscScalar,pointer :: xx(:), ff(:) 301 302 PetscCall(TSGetDM(ts,da,ierr)) 303 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 304 305 ! Scatter ghost points to local vector,using the 2-step process 306 ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 307 ! By placing code between these two statements, computations can be 308 ! done while messages are in transition. 309 PetscCall(DMGetLocalVector(da,Xloc,ierr)) 310 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)) 311 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)) 312 313 ! Get access to vector data 314 PetscCall(VecGetArrayRead(Xloc,xx,ierr)) 315 PetscCall(VecGetArray(F,ff,ierr)) 316 317 PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx,ff,user(user_a),user(user_k),user(user_s),ierr)) 318 319 PetscCall(VecRestoreArrayRead(Xloc,xx,ierr)) 320 PetscCall(VecRestoreArray(F,ff,ierr)) 321 PetscCall(DMRestoreLocalVector(da,Xloc,ierr)) 322end subroutine FormRHSFunction 323 324! --------------------------------------------------------------------- 325! 326! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 327! 328subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 329 use petscts 330 use petscdm 331 implicit none 332 333 TS ts 334 PetscReal t,shift 335 Vec X,Xdot 336 Mat J,Jpre 337 PetscReal user(6) 338 PetscErrorCode ierr 339 integer user_a,user_k,user_s 340 parameter (user_a = 0,user_k = 2,user_s = 4) 341 342 DM da 343 PetscInt mx,xs,xe,gxs,gxe 344 PetscInt i,i1,row,col 345 PetscReal k1,k2; 346 PetscScalar val(4) 347 348 PetscCall(TSGetDM(ts,da,ierr)) 349 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 350 351 i1 = 1 352 k1 = user(user_k+1) 353 k2 = user(user_k+2) 354 do i = xs,xe 355 row = i-gxs 356 col = i-gxs 357 val(1) = shift + k1 358 val(2) = -k2 359 val(3) = -k1 360 val(4) = shift + k2 361 PetscCall(MatSetValuesBlockedLocal(Jpre,i1,[row],i1,[col],val,INSERT_VALUES,ierr)) 362 end do 363 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 364 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 365 if (J /= Jpre) then 366 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)) 367 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)) 368 end if 369end subroutine FormIJacobian 370 371subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr) 372 implicit none 373 PetscInt mx,xs,xe,gxs,gxe 374 PetscScalar x(2,xs:xe) 375 PetscReal a(2),k(2),s(2) 376 PetscErrorCode ierr 377 378 PetscInt i 379 PetscReal one,hx,r,ik 380 one = 1.0 381 hx = one / mx 382 do i=xs,xe 383 r = i*hx 384 if (k(2) .ne. 0.0) then 385 ik = one/k(2) 386 else 387 ik = one 388 end if 389 x(1,i) = one + s(2)*r 390 x(2,i) = k(1)*ik*x(1,i) + s(2)*ik 391 end do 392end subroutine FormInitialSolutionLocal 393 394subroutine FormInitialSolution(ts,X,user,ierr) 395 use petscts 396 use petscdm 397 implicit none 398 399 TS ts 400 Vec X 401 PetscReal user(6) 402 PetscErrorCode ierr 403 integer user_a,user_k,user_s 404 parameter (user_a = 1,user_k = 3,user_s = 5) 405 406 DM da 407 PetscInt mx,xs,xe,gxs,gxe 408 PetscScalar,pointer :: xx(:) 409 410 PetscCall(TSGetDM(ts,da,ierr)) 411 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 412 413 ! Get access to vector data 414 PetscCall(VecGetArray(X,xx,ierr)) 415 416 PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx,user(user_a),user(user_k),user(user_s),ierr)) 417 418 PetscCall(VecRestoreArray(X,xx,ierr)) 419end subroutine FormInitialSolution 420 421! --------------------------------------------------------------------- 422! 423! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 424! 425subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 426 use ex22f_mfmodule 427 implicit none 428 TS ts 429 PetscReal t,shift 430 Vec X,Xdot 431 Mat J,Jpre 432 PetscReal user(6) 433 PetscErrorCode ierr 434 435 PETSC_SHIFT=shift 436 MFuser=user 437 438end subroutine FormIJacobianMF 439 440! ------------------------------------------------------------------- 441! 442! MyMult - user provided matrix multiply 443! 444! Input Parameters: 445!. X - input vector 446! 447! Output Parameter: 448!. F - function vector 449! 450subroutine MyMult(A,X,F,ierr) 451 use ex22f_mfmodule 452 implicit none 453 454 Mat A 455 Vec X,F 456 457 PetscErrorCode ierr 458 PetscScalar shift 459 460! Mat J,Jpre 461 462 PetscReal user(6) 463 464 integer user_a,user_k,user_s 465 parameter (user_a = 0,user_k = 2,user_s = 4) 466 467 DM da 468 PetscInt mx,xs,xe,gxs,gxe 469 PetscInt i,i1,row,col 470 PetscReal k1,k2; 471 PetscScalar val(4) 472 473 shift=PETSC_SHIFT 474 user=MFuser 475 476 PetscCall(TSGetDM(tscontext,da,ierr)) 477 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 478 479 i1 = 1 480 k1 = user(user_k+1) 481 k2 = user(user_k+2) 482 483 do i = xs,xe 484 row = i-gxs 485 col = i-gxs 486 val(1) = shift + k1 487 val(2) = -k2 488 val(3) = -k1 489 val(4) = shift + k2 490 PetscCall(MatSetValuesBlockedLocal(Jmat,i1,[row],i1,[col],val,INSERT_VALUES,ierr)) 491 end do 492 493! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 494! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 495! if (J /= Jpre) then 496 PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr)) 497 PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr)) 498! end if 499 500 PetscCall(MatMult(Jmat,X,F,ierr)) 501 502end subroutine MyMult 503 504! 505subroutine SaveSolutionToDisk(da,X,gdof,xs,xe) 506 use petscdm 507 implicit none 508 509 Vec X 510 DM da 511 PetscInt xs,xe,two 512 PetscInt gdof,i 513 PetscErrorCode ierr 514 PetscScalar data2(2,xs:xe),data(gdof) 515 PetscScalar,pointer :: xx(:) 516 517 PetscCall(VecGetArrayRead(X,xx,ierr)) 518 519 two = 2 520 data2=reshape(xx(gdof:gdof),(/two,xe-xs+1/)) 521 data=reshape(data2,(/gdof/)) 522 open(1020,file='solution_out_ex22f_mf.txt') 523 do i=1,gdof 524 write(1020,'(e24.16,1x)') data(i) 525 end do 526 close(1020) 527 528 PetscCall(VecRestoreArrayRead(X,xx,ierr)) 529end subroutine SaveSolutionToDisk 530 531!/*TEST 532! 533! test: 534! args: -da_grid_x 200 -ts_arkimex_type 4 535! requires: !single 536! output_file: output/ex22f_mf_1.out 537! 538!TEST*/ 539