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