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#include <petsc/finclude/petscts.h> 18#include <petsc/finclude/petscdmda.h> 19 20module ex22f_modctx 21 use petscts 22 type AppCtx 23 PetscReal a(2), k(2), s(2) 24 end type AppCtx 25 26end module ex22f_modctx 27program main 28 use ex22f_modctx 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 38 TS ts 39 SNES snes 40 SNESLineSearch linesearch 41 Vec X 42 Mat J 43 PetscInt mx 44 PetscErrorCode ierr 45 DM da 46 PetscReal ftime, dt 47 PetscReal one, pone 48 PetscInt im11, i2 49 PetscBool flg 50 type(AppCtx) ctx 51 52 im11 = 11 53 i2 = 2 54 one = 1.0 55 pone = one/10 56 57 PetscCallA(PetscInitialize(ierr)) 58 59! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60! Create distributed array (DMDA) to manage parallel grid and vectors 61! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, im11, i2, i2, PETSC_NULL_INTEGER, da, ierr)) 63 PetscCallA(DMSetFromOptions(da, ierr)) 64 PetscCallA(DMSetUp(da, ierr)) 65 66! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67! Extract global vectors from DMDA 68! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 69 PetscCallA(DMCreateGlobalVector(da, X, ierr)) 70 71! Initialize user application context 72! Use zero-based indexing for command line parameters to match ex22.c 73 ctx%a(1) = 1.0 74 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a0', ctx%a(1), flg, ierr)) 75 ctx%a(2) = 0.0 76 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-a1', ctx%a(2), flg, ierr)) 77 ctx%k(1) = 1000000.0 78 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k0', ctx%k(1), flg, ierr)) 79 ctx%k(2) = 2*ctx%k(1) 80 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-k1', ctx%k(2), flg, ierr)) 81 ctx%s(1) = 0.0 82 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s0', ctx%s(1), flg, ierr)) 83 ctx%s(2) = 1.0 84 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-s1', ctx%s(2), flg, ierr)) 85 86! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 87! Create timestepping solver context 88! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 89 PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr)) 90 PetscCallA(TSSetDM(ts, da, ierr)) 91 PetscCallA(TSSetType(ts, TSARKIMEX, ierr)) 92 PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormRHSFunction, ctx, ierr)) 93 PetscCallA(TSSetIFunction(ts, PETSC_NULL_VEC, FormIFunction, ctx, ierr)) 94 PetscCallA(DMSetMatType(da, MATAIJ, ierr)) 95 PetscCallA(DMCreateMatrix(da, J, ierr)) 96 PetscCallA(TSSetIJacobian(ts, J, J, FormIJacobian, ctx, ierr)) 97 98 PetscCallA(TSGetSNES(ts, snes, ierr)) 99 PetscCallA(SNESGetLineSearch(snes, linesearch, ierr)) 100 PetscCallA(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC, ierr)) 101 102 ftime = 1.0 103 PetscCallA(TSSetMaxTime(ts, ftime, ierr)) 104 PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr)) 105 106! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 107! Set initial conditions 108! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 PetscCallA(FormInitialSolution(ts, X, ctx, ierr)) 110 PetscCallA(TSSetSolution(ts, X, ierr)) 111 PetscCallA(VecGetSize(X, mx, ierr)) 112! Advective CFL, I don't know why it needs so much safety factor. 113 dt = pone*max(ctx%a(1), ctx%a(2))/mx 114 PetscCallA(TSSetTimeStep(ts, dt, ierr)) 115 116! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 117! Set runtime options 118! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 PetscCallA(TSSetFromOptions(ts, ierr)) 120 121! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122! Solve nonlinear system 123! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124 PetscCallA(TSSolve(ts, X, ierr)) 125 126! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127! Free work space. 128! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 PetscCallA(MatDestroy(J, ierr)) 130 PetscCallA(VecDestroy(X, ierr)) 131 PetscCallA(TSDestroy(ts, ierr)) 132 PetscCallA(DMDestroy(da, ierr)) 133 PetscCallA(PetscFinalize(ierr)) 134contains 135 136! Small helper to extract the layout, result uses 1-based indexing. 137 subroutine GetLayout(da, mx, xs, xe, gxs, gxe, ierr) 138 use petscdm 139 implicit none 140 141 DM da 142 PetscInt mx, xs, xe, gxs, gxe 143 PetscErrorCode ierr 144 PetscInt xm, gxm 145 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)) 146 PetscCall(DMDAGetCorners(da, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)) 147 PetscCall(DMDAGetGhostCorners(da, gxs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, gxm, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)) 148 xs = xs + 1 149 gxs = gxs + 1 150 xe = xs + xm - 1 151 gxe = gxs + gxm - 1 152 end subroutine 153 154 subroutine FormIFunctionLocal(mx, xs, xe, gxs, gxe, x, xdot, f, a, k, s, ierr) 155 implicit none 156 PetscInt mx, xs, xe, gxs, gxe 157 PetscScalar x(2, xs:xe) 158 PetscScalar xdot(2, xs:xe) 159 PetscScalar f(2, xs:xe) 160 PetscReal a(2), k(2), s(2) 161 PetscErrorCode ierr 162 PetscInt i 163 do i = xs, xe 164 f(1, i) = xdot(1, i) + k(1)*x(1, i) - k(2)*x(2, i) - s(1) 165 f(2, i) = xdot(2, i) - k(1)*x(1, i) + k(2)*x(2, i) - s(2) 166 end do 167 end subroutine 168 169 subroutine FormIFunction(ts, t, X, Xdot, F, ctx, ierr) 170 use ex22f_modctx 171 implicit none 172 173 TS ts 174 type(AppCtx) ctx 175 PetscReal t 176 Vec X, Xdot, F 177 PetscErrorCode ierr 178 179 DM da 180 PetscInt mx, xs, xe, gxs, gxe 181 PetscScalar, pointer :: xx(:), xxdot(:), ff(:) 182 183 PetscCall(TSGetDM(ts, da, ierr)) 184 PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr)) 185 186! Get access to vector data 187 PetscCall(VecGetArrayRead(X, xx, ierr)) 188 PetscCall(VecGetArrayRead(Xdot, xxdot, ierr)) 189 PetscCall(VecGetArray(F, ff, ierr)) 190 191 PetscCall(FormIFunctionLocal(mx, xs, xe, gxs, gxe, xx, xxdot, ff, ctx%a, ctx%k, ctx%s, ierr)) 192 193 PetscCall(VecRestoreArrayRead(X, xx, ierr)) 194 PetscCall(VecRestoreArrayRead(Xdot, xxdot, ierr)) 195 PetscCall(VecRestoreArray(F, ff, ierr)) 196 end subroutine 197 198 subroutine FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, x, f, a, k, s, ierr) 199 implicit none 200 PetscInt mx, xs, xe, gxs, gxe 201 PetscReal t 202 PetscScalar x(2, gxs:gxe), f(2, xs:xe) 203 PetscReal a(2), k(2), s(2) 204 PetscErrorCode ierr 205 PetscInt i, j 206 PetscReal hx, u0t(2) 207 PetscReal one, two, three, four, six, twelve 208 PetscReal half, third, twothird, sixth 209 PetscReal twelfth 210 211 one = 1.0 212 two = 2.0 213 three = 3.0 214 four = 4.0 215 six = 6.0 216 twelve = 12.0 217 hx = one/mx 218! The Fortran standard only allows positive base for power functions; Nag compiler fails on this 219 u0t(1) = one - abs(sin(twelve*t))**four 220 u0t(2) = 0.0 221 half = one/two 222 third = one/three 223 twothird = two/three 224 sixth = one/six 225 twelfth = one/twelve 226 do i = xs, xe 227 do j = 1, 2 228 if (i == 1) then 229 f(j, i) = a(j)/hx*(third*u0t(j) + half*x(j, i) - x(j, i + 1) & 230& + sixth*x(j, i + 2)) 231 else if (i == 2) then 232 f(j, i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j, i - 1) & 233& - twothird*x(j, i + 1) + twelfth*x(j, i + 2)) 234 else if (i == mx - 1) then 235 f(j, i) = a(j)/hx*(-sixth*x(j, i - 2) + x(j, i - 1) & 236& - half*x(j, i) - third*x(j, i + 1)) 237 else if (i == mx) then 238 f(j, i) = a(j)/hx*(-x(j, i) + x(j, i - 1)) 239 else 240 f(j, i) = a(j)/hx*(-twelfth*x(j, i - 2) & 241& + twothird*x(j, i - 1) & 242& - twothird*x(j, i + 1) + twelfth*x(j, i + 2)) 243 end if 244 end do 245 end do 246 end subroutine 247 248 subroutine FormRHSFunction(ts, t, X, F, ctx, ierr) 249 use ex22f_modctx 250 implicit none 251 252 type(AppCtx) ctx 253 TS ts 254 PetscReal t 255 Vec X, F 256 PetscErrorCode ierr 257 DM da 258 Vec Xloc 259 PetscInt mx, xs, xe, gxs, gxe 260 PetscScalar, pointer :: xx(:), ff(:) 261 262 PetscCall(TSGetDM(ts, da, ierr)) 263 PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr)) 264 265! Scatter ghost points to local vector,using the 2-step process 266! DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 267! By placing code between these two statements, computations can be 268! done while messages are in transition. 269 PetscCall(DMGetLocalVector(da, Xloc, ierr)) 270 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc, ierr)) 271 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc, ierr)) 272 273! Get access to vector data 274 PetscCall(VecGetArrayRead(Xloc, xx, ierr)) 275 PetscCall(VecGetArray(F, ff, ierr)) 276 277 PetscCall(FormRHSFunctionLocal(mx, xs, xe, gxs, gxe, t, xx, ff, ctx%a, ctx%k, ctx%s, ierr)) 278 279 PetscCall(VecRestoreArrayRead(Xloc, xx, ierr)) 280 PetscCall(VecRestoreArray(F, ff, ierr)) 281 PetscCall(DMRestoreLocalVector(da, Xloc, ierr)) 282 end subroutine 283 284! --------------------------------------------------------------------- 285! 286! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 287! 288 subroutine FormIJacobian(ts, t, X, Xdot, shift, J, Jpre, ctx, ierr) 289 use ex22f_modctx 290 implicit none 291 292 type(AppCtx) ctx 293 TS ts 294 PetscReal t, shift 295 Vec X, Xdot 296 Mat J, Jpre 297 PetscErrorCode ierr 298 299 DM da 300 PetscInt mx, xs, xe, gxs, gxe 301 PetscInt i, i1, row, col 302 PetscReal k1, k2 303 PetscScalar val(4) 304 305 PetscCall(TSGetDM(ts, da, ierr)) 306 PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr)) 307 308 i1 = 1 309 k1 = ctx%k(1) 310 k2 = ctx%k(2) 311 do i = xs, xe 312 row = i - gxs 313 col = i - gxs 314 val(1) = shift + k1 315 val(2) = -k2 316 val(3) = -k1 317 val(4) = shift + k2 318 PetscCall(MatSetValuesBlockedLocal(Jpre, i1, [row], i1, [col], val, INSERT_VALUES, ierr)) 319 end do 320 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY, ierr)) 321 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY, ierr)) 322 if (J /= Jpre) then 323 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr)) 324 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr)) 325 end if 326 end subroutine 327 328 subroutine FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, x, a, k, s, ierr) 329 implicit none 330 PetscInt mx, xs, xe, gxs, gxe 331 PetscScalar x(2, xs:xe) 332 PetscReal a(2), k(2), s(2) 333 PetscErrorCode ierr 334 335 PetscInt i 336 PetscReal one, hx, r, ik 337 one = 1.0 338 hx = one/mx 339 do i = xs, xe 340 r = i*hx 341 if (k(2) /= 0.0) then 342 ik = one/k(2) 343 else 344 ik = one 345 end if 346 x(1, i) = one + s(2)*r 347 x(2, i) = k(1)*ik*x(1, i) + s(2)*ik 348 end do 349 end subroutine 350 351 subroutine FormInitialSolution(ts, X, ctx, ierr) 352 use ex22f_modctx 353 implicit none 354 355 type(AppCtx) ctx 356 TS ts 357 Vec X 358 PetscErrorCode ierr 359 360 DM da 361 PetscInt mx, xs, xe, gxs, gxe 362 PetscScalar, pointer :: xx(:) 363 364 PetscCall(TSGetDM(ts, da, ierr)) 365 PetscCall(GetLayout(da, mx, xs, xe, gxs, gxe, ierr)) 366 367! Get access to vector data 368 PetscCall(VecGetArray(X, xx, ierr)) 369 370 PetscCall(FormInitialSolutionLocal(mx, xs, xe, gxs, gxe, xx, ctx%a, ctx%k, ctx%s, ierr)) 371 372 PetscCall(VecRestoreArray(X, xx, ierr)) 373 end subroutine 374end program 375!/*TEST 376! 377! test: 378! args: -da_grid_x 200 -ts_arkimex_type 4 379! requires: !single 380! output_file: output/empty.out 381! 382!TEST*/ 383