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 PetscCallA(PetscInitialize(ierr)) 55 56! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 57! Create distributed array (DMDA) to manage parallel grid and vectors 58! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr)) 60 PetscCallA(DMSetFromOptions(da,ierr)) 61 PetscCallA(DMSetUp(da,ierr)) 62 63! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 64! Extract global vectors from DMDA; 65! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66 PetscCallA(DMCreateGlobalVector(da,X,ierr)) 67 68! Initialize user application context 69! Use zero-based indexing for command line parameters to match ex22.c 70 user(user_a+1) = 1.0 71 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr)) 72 user(user_a+2) = 0.0 73 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr)) 74 user(user_k+1) = 1000000.0 75 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0',user(user_k+1),flg,ierr)) 76 user(user_k+2) = 2*user(user_k+1) 77 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr)) 78 user(user_s+1) = 0.0 79 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr)) 80 user(user_s+2) = 1.0 81 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr)) 82 83! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84! Create timestepping solver context 85! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr)) 87 PetscCallA(TSSetDM(ts,da,ierr)) 88 PetscCallA(TSSetType(ts,TSARKIMEX,ierr)) 89 PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr)) 90 PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr)) 91 PetscCallA(DMSetMatType(da,MATAIJ,ierr)) 92 PetscCallA(DMCreateMatrix(da,J,ierr)) 93 PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)) 94 95 PetscCallA(TSGetSNES(ts,snes,ierr)) 96 PetscCallA(SNESGetLineSearch(snes,linesearch,ierr)) 97 PetscCallA(SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC,ierr)) 98 99 ftime = 1.0 100 PetscCallA(TSSetMaxTime(ts,ftime,ierr)) 101 PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)) 102 103! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 104! Set initial conditions 105! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 106 PetscCallA(FormInitialSolution(ts,X,user,ierr)) 107 PetscCallA(TSSetSolution(ts,X,ierr)) 108 PetscCallA(VecGetSize(X,mx,ierr)) 109! Advective CFL, I don't know why it needs so much safety factor. 110 dt = pone * max(user(user_a+1),user(user_a+2)) / mx; 111 PetscCallA(TSSetTimeStep(ts,dt,ierr)) 112 113! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114! Set runtime options 115! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116 PetscCallA(TSSetFromOptions(ts,ierr)) 117 118! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119! Solve nonlinear system 120! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 PetscCallA(TSSolve(ts,X,ierr)) 122 123! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124! Free work space. 125! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126 PetscCallA(MatDestroy(J,ierr)) 127 PetscCallA(VecDestroy(X,ierr)) 128 PetscCallA(TSDestroy(ts,ierr)) 129 PetscCallA(DMDestroy(da,ierr)) 130 PetscCallA(PetscFinalize(ierr)) 131 end program 132 133! Small helper to extract the layout, result uses 1-based indexing. 134 subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr) 135 use petscdmda 136 implicit none 137 138 DM da 139 PetscInt mx,xs,xe,gxs,gxe 140 PetscErrorCode ierr 141 PetscInt xm,gxm 142 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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 143 PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 144 PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 145 xs = xs + 1 146 gxs = gxs + 1 147 xe = xs + xm - 1 148 gxe = gxs + gxm - 1 149 end subroutine 150 151 subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr) 152 implicit none 153 PetscInt mx,xs,xe,gxs,gxe 154 PetscScalar x(2,xs:xe) 155 PetscScalar xdot(2,xs:xe) 156 PetscScalar f(2,xs:xe) 157 PetscReal a(2),k(2),s(2) 158 PetscErrorCode ierr 159 PetscInt i 160 do 10 i = xs,xe 161 f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1) 162 f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2) 163 10 continue 164 end subroutine 165 166 subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr) 167 use petscts 168 implicit none 169 170 TS ts 171 PetscReal t 172 Vec X,Xdot,F 173 PetscReal user(6) 174 PetscErrorCode ierr 175 integer user_a,user_k,user_s 176 parameter (user_a = 1,user_k = 3,user_s = 5) 177 178 DM da 179 PetscInt mx,xs,xe,gxs,gxe 180 PetscOffset ixx,ixxdot,iff 181 PetscScalar xx(0:1),xxdot(0:1),ff(0:1) 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,ixx,ierr)) 188 PetscCall(VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)) 189 PetscCall(VecGetArray(F,ff,iff,ierr)) 190 191 PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr)) 192 193 PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr)) 194 PetscCall(VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)) 195 PetscCall(VecRestoreArray(F,ff,iff,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 20 i = xs,xe 227 do 10 j = 1,2 228 if (i .eq. 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 .eq. 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 .eq. 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 .eq. 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 10 continue 245 20 continue 246 end subroutine 247 248 subroutine FormRHSFunction(ts,t,X,F,user,ierr) 249 use petscts 250 implicit none 251 252 TS ts 253 PetscReal t 254 Vec X,F 255 PetscReal user(6) 256 PetscErrorCode ierr 257 integer user_a,user_k,user_s 258 parameter (user_a = 1,user_k = 3,user_s = 5) 259 DM da 260 Vec Xloc 261 PetscInt mx,xs,xe,gxs,gxe 262 PetscOffset ixx,iff 263 PetscScalar xx(0:1),ff(0:1) 264 265 PetscCall(TSGetDM(ts,da,ierr)) 266 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 267 268! Scatter ghost points to local vector,using the 2-step process 269! DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 270! By placing code between these two statements, computations can be 271! done while messages are in transition. 272 PetscCall(DMGetLocalVector(da,Xloc,ierr)) 273 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)) 274 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)) 275 276! Get access to vector data 277 PetscCall(VecGetArrayRead(Xloc,xx,ixx,ierr)) 278 PetscCall(VecGetArray(F,ff,iff,ierr)) 279 280 PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx(ixx),ff(iff),user(user_a),user(user_k),user(user_s),ierr)) 281 282 PetscCall(VecRestoreArrayRead(Xloc,xx,ixx,ierr)) 283 PetscCall(VecRestoreArray(F,ff,iff,ierr)) 284 PetscCall(DMRestoreLocalVector(da,Xloc,ierr)) 285 end subroutine 286 287! --------------------------------------------------------------------- 288! 289! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 290! 291 subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 292 use petscts 293 implicit none 294 295 TS ts 296 PetscReal t,shift 297 Vec X,Xdot 298 Mat J,Jpre 299 PetscReal user(6) 300 PetscErrorCode ierr 301 integer user_a,user_k,user_s 302 parameter (user_a = 0,user_k = 2,user_s = 4) 303 304 DM da 305 PetscInt mx,xs,xe,gxs,gxe 306 PetscInt i,i1,row,col 307 PetscReal k1,k2; 308 PetscScalar val(4) 309 310 PetscCall(TSGetDM(ts,da,ierr)) 311 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 312 313 i1 = 1 314 k1 = user(user_k+1) 315 k2 = user(user_k+2) 316 do 10 i = xs,xe 317 row = i-gxs 318 col = i-gxs 319 val(1) = shift + k1 320 val(2) = -k2 321 val(3) = -k1 322 val(4) = shift + k2 323 PetscCall(MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,INSERT_VALUES,ierr)) 324 10 continue 325 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 326 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 327 if (J /= Jpre) then 328 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)) 329 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)) 330 end if 331 end subroutine 332 333 subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr) 334 implicit none 335 PetscInt mx,xs,xe,gxs,gxe 336 PetscScalar x(2,xs:xe) 337 PetscReal a(2),k(2),s(2) 338 PetscErrorCode ierr 339 340 PetscInt i 341 PetscReal one,hx,r,ik 342 one = 1.0 343 hx = one / mx 344 do 10 i=xs,xe 345 r = i*hx 346 if (k(2) .ne. 0.0) then 347 ik = one/k(2) 348 else 349 ik = one 350 end if 351 x(1,i) = one + s(2)*r 352 x(2,i) = k(1)*ik*x(1,i) + s(2)*ik 353 10 continue 354 end subroutine 355 356 subroutine FormInitialSolution(ts,X,user,ierr) 357 use petscts 358 implicit none 359 360 TS ts 361 Vec X 362 PetscReal user(6) 363 PetscErrorCode ierr 364 integer user_a,user_k,user_s 365 parameter (user_a = 1,user_k = 3,user_s = 5) 366 367 DM da 368 PetscInt mx,xs,xe,gxs,gxe 369 PetscOffset ixx 370 PetscScalar xx(0:1) 371 372 PetscCall(TSGetDM(ts,da,ierr)) 373 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 374 375! Get access to vector data 376 PetscCall(VecGetArray(X,xx,ixx,ierr)) 377 378 PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx(ixx),user(user_a),user(user_k),user(user_s),ierr)) 379 380 PetscCall(VecRestoreArray(X,xx,ixx,ierr)) 381 end subroutine 382 383!/*TEST 384! 385! test: 386! args: -da_grid_x 200 -ts_arkimex_type 4 387! requires: !single 388! 389!TEST*/ 390