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