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