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 petscdm 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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,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 PetscScalar,pointer :: xx(:),xxdot(:),ff(:) 181 182 PetscCall(TSGetDM(ts,da,ierr)) 183 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 184 185! Get access to vector data 186 PetscCall(VecGetArrayRead(X,xx,ierr)) 187 PetscCall(VecGetArrayRead(Xdot,xxdot,ierr)) 188 PetscCall(VecGetArray(F,ff,ierr)) 189 190 PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx,xxdot,ff,user(user_a),user(user_k),user(user_s),ierr)) 191 192 PetscCall(VecRestoreArrayRead(X,xx,ierr)) 193 PetscCall(VecRestoreArrayRead(Xdot,xxdot,ierr)) 194 PetscCall(VecRestoreArray(F,ff,ierr)) 195 end subroutine 196 197 subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr) 198 implicit none 199 PetscInt mx,xs,xe,gxs,gxe 200 PetscReal t 201 PetscScalar x(2,gxs:gxe),f(2,xs:xe) 202 PetscReal a(2),k(2),s(2) 203 PetscErrorCode ierr 204 PetscInt i,j 205 PetscReal hx,u0t(2) 206 PetscReal one,two,three,four,six,twelve 207 PetscReal half,third,twothird,sixth 208 PetscReal twelfth 209 210 one = 1.0 211 two = 2.0 212 three = 3.0 213 four = 4.0 214 six = 6.0 215 twelve = 12.0 216 hx = one / mx 217! The Fortran standard only allows positive base for power functions; Nag compiler fails on this 218 u0t(1) = one - abs(sin(twelve*t))**four 219 u0t(2) = 0.0 220 half = one/two 221 third = one / three 222 twothird = two / three 223 sixth = one / six 224 twelfth = one / twelve 225 do 20 i = xs,xe 226 do 10 j = 1,2 227 if (i .eq. 1) then 228 f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) & 229 & + sixth*x(j,i+2)) 230 else if (i .eq. 2) then 231 f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) & 232 & - twothird*x(j,i+1) + twelfth*x(j,i+2)) 233 else if (i .eq. mx-1) then 234 f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) & 235 & - half*x(j,i) -third*x(j,i+1)) 236 else if (i .eq. mx) then 237 f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1)) 238 else 239 f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) & 240 & + twothird*x(j,i-1) & 241 & - twothird*x(j,i+1) + twelfth*x(j,i+2)) 242 end if 243 10 continue 244 20 continue 245 end subroutine 246 247 subroutine FormRHSFunction(ts,t,X,F,user,ierr) 248 use petscts 249 implicit none 250 251 TS ts 252 PetscReal t 253 Vec X,F 254 PetscReal user(6) 255 PetscErrorCode ierr 256 integer user_a,user_k,user_s 257 parameter (user_a = 1,user_k = 3,user_s = 5) 258 DM da 259 Vec Xloc 260 PetscInt mx,xs,xe,gxs,gxe 261 PetscScalar,pointer :: xx(:),ff(:) 262 263 PetscCall(TSGetDM(ts,da,ierr)) 264 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 265 266! Scatter ghost points to local vector,using the 2-step process 267! DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 268! By placing code between these two statements, computations can be 269! done while messages are in transition. 270 PetscCall(DMGetLocalVector(da,Xloc,ierr)) 271 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)) 272 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)) 273 274! Get access to vector data 275 PetscCall(VecGetArrayRead(Xloc,xx,ierr)) 276 PetscCall(VecGetArray(F,ff,ierr)) 277 278 PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx,ff,user(user_a),user(user_k),user(user_s),ierr)) 279 280 PetscCall(VecRestoreArrayRead(Xloc,xx,ierr)) 281 PetscCall(VecRestoreArray(F,ff,ierr)) 282 PetscCall(DMRestoreLocalVector(da,Xloc,ierr)) 283 end subroutine 284 285! --------------------------------------------------------------------- 286! 287! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 288! 289 subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 290 use petscts 291 implicit none 292 293 TS ts 294 PetscReal t,shift 295 Vec X,Xdot 296 Mat J,Jpre 297 PetscReal user(6) 298 PetscErrorCode ierr 299 integer user_a,user_k,user_s 300 parameter (user_a = 0,user_k = 2,user_s = 4) 301 302 DM da 303 PetscInt mx,xs,xe,gxs,gxe 304 PetscInt i,i1,row,col 305 PetscReal k1,k2; 306 PetscScalar val(4) 307 308 PetscCall(TSGetDM(ts,da,ierr)) 309 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 310 311 i1 = 1 312 k1 = user(user_k+1) 313 k2 = user(user_k+2) 314 do 10 i = xs,xe 315 row = i-gxs 316 col = i-gxs 317 val(1) = shift + k1 318 val(2) = -k2 319 val(3) = -k1 320 val(4) = shift + k2 321 PetscCall(MatSetValuesBlockedLocal(Jpre,i1,[row],i1,[col],val,INSERT_VALUES,ierr)) 322 10 continue 323 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 324 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 325 if (J /= Jpre) then 326 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)) 327 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)) 328 end if 329 end subroutine 330 331 subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr) 332 implicit none 333 PetscInt mx,xs,xe,gxs,gxe 334 PetscScalar x(2,xs:xe) 335 PetscReal a(2),k(2),s(2) 336 PetscErrorCode ierr 337 338 PetscInt i 339 PetscReal one,hx,r,ik 340 one = 1.0 341 hx = one / mx 342 do 10 i=xs,xe 343 r = i*hx 344 if (k(2) .ne. 0.0) then 345 ik = one/k(2) 346 else 347 ik = one 348 end if 349 x(1,i) = one + s(2)*r 350 x(2,i) = k(1)*ik*x(1,i) + s(2)*ik 351 10 continue 352 end subroutine 353 354 subroutine FormInitialSolution(ts,X,user,ierr) 355 use petscts 356 implicit none 357 358 TS ts 359 Vec X 360 PetscReal user(6) 361 PetscErrorCode ierr 362 integer user_a,user_k,user_s 363 parameter (user_a = 1,user_k = 3,user_s = 5) 364 365 DM da 366 PetscInt mx,xs,xe,gxs,gxe 367 PetscScalar,pointer :: xx(:) 368 369 PetscCall(TSGetDM(ts,da,ierr)) 370 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 371 372! Get access to vector data 373 PetscCall(VecGetArray(X,xx,ierr)) 374 375 PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx,user(user_a),user(user_k),user(user_s),ierr)) 376 377 PetscCall(VecRestoreArray(X,xx,ierr)) 378 end subroutine 379 380!/*TEST 381! 382! test: 383! args: -da_grid_x 200 -ts_arkimex_type 4 384! requires: !single 385! 386!TEST*/ 387