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