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 module ex22f_mfmodule 18#include <petsc/finclude/petscts.h> 19 use petscts 20 PetscScalar::PETSC_SHIFT 21 TS::tscontext 22 Mat::Jmat 23 PetscReal::MFuser(6) 24 end module ex22f_mfmodule 25 26program main 27 use ex22f_mfmodule 28 use petscdmda 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 PetscReal user(6) 38 integer user_a,user_k,user_s 39 parameter (user_a = 0,user_k = 2,user_s = 4) 40 41 external FormRHSFunction,FormIFunction 42 external FormInitialSolution 43 external FormIJacobian 44 external MyMult,FormIJacobianMF 45 46 TS ts 47 Vec X 48 Mat J 49 PetscInt mx 50 PetscBool OptionSaveToDisk 51 PetscErrorCode ierr 52 DM da 53 PetscReal ftime,dt 54 PetscReal one,pone 55 PetscInt im11,i2 56 PetscBool flg 57 58 PetscInt xs,xe,gxs,gxe,dof,gdof 59 PetscScalar shell_shift 60 Mat A 61 62 im11 = 11 63 i2 = 2 64 one = 1.0 65 pone = one / 10 66 67 PetscCallA(PetscInitialize(ierr)) 68 69 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 70 ! Create distributed array (DMDA) to manage parallel grid and vectors 71 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 72 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr)) 73 PetscCallA(DMSetFromOptions(da,ierr)) 74 PetscCallA(DMSetUp(da,ierr)) 75 76 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 ! Extract global vectors from DMDA; 78 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 PetscCallA(DMCreateGlobalVector(da,X,ierr)) 80 81 ! Initialize user application context 82 ! Use zero-based indexing for command line parameters to match ex22.c 83 user(user_a+1) = 1.0 84 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a0',user(user_a+1),flg,ierr)) 85 user(user_a+2) = 0.0 86 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-a1',user(user_a+2),flg,ierr)) 87 user(user_k+1) = 1000000.0 88 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k0', user(user_k+1),flg,ierr)) 89 user(user_k+2) = 2*user(user_k+1) 90 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-k1', user(user_k+2),flg,ierr)) 91 user(user_s+1) = 0.0 92 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s0',user(user_s+1),flg,ierr)) 93 user(user_s+2) = 1.0 94 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-s1',user(user_s+2),flg,ierr)) 95 96 OptionSaveToDisk=.FALSE. 97 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sdisk',OptionSaveToDisk,flg,ierr)) 98 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 ! Create timestepping solver context 100 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr)) 102 tscontext=ts 103 PetscCallA(TSSetDM(ts,da,ierr)) 104 PetscCallA(TSSetType(ts,TSARKIMEX,ierr)) 105 PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormRHSFunction,user,ierr)) 106 107 ! - - - - - - - - -- - - - - 108 ! Matrix free setup 109 PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 110 dof=i2*(xe-xs+1) 111 gdof=i2*(gxe-gxs+1) 112 PetscCallA(MatCreateShell(PETSC_COMM_WORLD,dof,dof,gdof,gdof,shell_shift,A,ierr)) 113 114 PetscCallA(MatShellSetOperation(A,MATOP_MULT,MyMult,ierr)) 115 ! - - - - - - - - - - - - 116 117 PetscCallA(TSSetIFunction(ts,PETSC_NULL_VEC,FormIFunction,user,ierr)) 118 PetscCallA(DMSetMatType(da,MATAIJ,ierr)) 119 PetscCallA(DMCreateMatrix(da,J,ierr)) 120 121 Jmat=J 122 123 PetscCallA(TSSetIJacobian(ts,J,J,FormIJacobian,user,ierr)) 124 PetscCallA(TSSetIJacobian(ts,A,A,FormIJacobianMF,user,ierr)) 125 126 ftime = 1.0 127 PetscCallA(TSSetMaxTime(ts,ftime,ierr)) 128 PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)) 129 130 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 ! Set initial conditions 132 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 PetscCallA(FormInitialSolution(ts,X,user,ierr)) 134 PetscCallA(TSSetSolution(ts,X,ierr)) 135 PetscCallA(VecGetSize(X,mx,ierr)) 136 ! Advective CFL, I don't know why it needs so much safety factor. 137 dt = pone * max(user(user_a+1),user(user_a+2)) / mx; 138 PetscCallA(TSSetTimeStep(ts,dt,ierr)) 139 140 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 ! Set runtime options 142 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 PetscCallA(TSSetFromOptions(ts,ierr)) 144 145 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 ! Solve nonlinear system 147 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 148 PetscCallA(TSSolve(ts,X,ierr)) 149 150 if (OptionSaveToDisk) then 151 PetscCallA(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 152 dof=i2*(xe-xs+1) 153 gdof=i2*(gxe-gxs+1) 154 call SaveSolutionToDisk(da,X,gdof,xs,xe) 155 end if 156 157 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 ! Free work space. 159! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 PetscCallA(MatDestroy(A,ierr)) 161 PetscCallA(MatDestroy(J,ierr)) 162 PetscCallA(VecDestroy(X,ierr)) 163 PetscCallA(TSDestroy(ts,ierr)) 164 PetscCallA(DMDestroy(da,ierr)) 165 PetscCallA(PetscFinalize(ierr)) 166end program main 167 168! Small helper to extract the layout, result uses 1-based indexing. 169 subroutine GetLayout(da,mx,xs,xe,gxs,gxe,ierr) 170 use petscdmda 171 implicit none 172 173 DM da 174 PetscInt mx,xs,xe,gxs,gxe 175 PetscErrorCode ierr 176 PetscInt xm,gxm 177 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)) 178 PetscCall(DMDAGetCorners(da,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 179 PetscCall(DMDAGetGhostCorners(da,gxs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,gxm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 180 xs = xs + 1 181 gxs = gxs + 1 182 xe = xs + xm - 1 183 gxe = gxs + gxm - 1 184end subroutine GetLayout 185 186subroutine FormIFunctionLocal(mx,xs,xe,gxs,gxe,x,xdot,f,a,k,s,ierr) 187 implicit none 188 PetscInt mx,xs,xe,gxs,gxe 189 PetscScalar x(2,xs:xe) 190 PetscScalar xdot(2,xs:xe) 191 PetscScalar f(2,xs:xe) 192 PetscReal a(2),k(2),s(2) 193 PetscErrorCode ierr 194 PetscInt i 195 do i = xs,xe 196 f(1,i) = xdot(1,i) + k(1)*x(1,i) - k(2)*x(2,i) - s(1) 197 f(2,i) = xdot(2,i) - k(1)*x(1,i) + k(2)*x(2,i) - s(2) 198 end do 199end subroutine FormIFunctionLocal 200 201subroutine FormIFunction(ts,t,X,Xdot,F,user,ierr) 202 use petscdmda 203 use petscts 204 implicit none 205 206 TS ts 207 PetscReal t 208 Vec X,Xdot,F 209 PetscReal user(6) 210 PetscErrorCode ierr 211 integer user_a,user_k,user_s 212 parameter (user_a = 1,user_k = 3,user_s = 5) 213 214 DM da 215 PetscInt mx,xs,xe,gxs,gxe 216 PetscOffset ixx,ixxdot,iff 217 PetscScalar xx(0:1),xxdot(0:1),ff(0:1) 218 219 PetscCall(TSGetDM(ts,da,ierr)) 220 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 221 222 ! Get access to vector data 223 PetscCall(VecGetArrayRead(X,xx,ixx,ierr)) 224 PetscCall(VecGetArrayRead(Xdot,xxdot,ixxdot,ierr)) 225 PetscCall(VecGetArray(F,ff,iff,ierr)) 226 227 PetscCall(FormIFunctionLocal(mx,xs,xe,gxs,gxe,xx(ixx),xxdot(ixxdot),ff(iff),user(user_a),user(user_k),user(user_s),ierr)) 228 229 PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr)) 230 PetscCall(VecRestoreArrayRead(Xdot,xxdot,ixxdot,ierr)) 231 PetscCall(VecRestoreArray(F,ff,iff,ierr)) 232end subroutine FormIFunction 233 234subroutine FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,x,f,a,k,s,ierr) 235 implicit none 236 PetscInt mx,xs,xe,gxs,gxe 237 PetscReal t 238 PetscScalar x(2,gxs:gxe),f(2,xs:xe) 239 PetscReal a(2),k(2),s(2) 240 PetscErrorCode ierr 241 PetscInt i,j 242 PetscReal hx,u0t(2) 243 PetscReal one,two,three,four,six,twelve 244 PetscReal half,third,twothird,sixth 245 PetscReal twelfth 246 247 one = 1.0 248 two = 2.0 249 three = 3.0 250 four = 4.0 251 six = 6.0 252 twelve = 12.0 253 hx = one / mx 254 u0t(1) = one - sin(twelve*t)**four 255 u0t(2) = 0.0 256 half = one/two 257 third = one / three 258 twothird = two / three 259 sixth = one / six 260 twelfth = one / twelve 261 do i = xs,xe 262 do j = 1,2 263 if (i .eq. 1) then 264 f(j,i) = a(j)/hx*(third*u0t(j) + half*x(j,i) - x(j,i+1) + sixth*x(j,i+2)) 265 else if (i .eq. 2) then 266 f(j,i) = a(j)/hx*(-twelfth*u0t(j) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2)) 267 else if (i .eq. mx-1) then 268 f(j,i) = a(j)/hx*(-sixth*x(j,i-2) + x(j,i-1) - half*x(j,i) -third*x(j,i+1)) 269 else if (i .eq. mx) then 270 f(j,i) = a(j)/hx*(-x(j,i) + x(j,i-1)) 271 else 272 f(j,i) = a(j)/hx*(-twelfth*x(j,i-2) + twothird*x(j,i-1) - twothird*x(j,i+1) + twelfth*x(j,i+2)) 273 end if 274 end do 275 end do 276 277#ifdef EXPLICIT_INTEGRATOR22 278 do i = xs,xe 279 f(1,i) = f(1,i) -( k(1)*x(1,i) - k(2)*x(2,i) - s(1)) 280 f(2,i) = f(2,i) -(- k(1)*x(1,i) + k(2)*x(2,i) - s(2)) 281 end do 282#endif 283 284end subroutine FormRHSFunctionLocal 285 286subroutine FormRHSFunction(ts,t,X,F,user,ierr) 287 use petscts 288 use petscdmda 289 implicit none 290 291 TS ts 292 PetscReal t 293 Vec X,F 294 PetscReal user(6) 295 PetscErrorCode ierr 296 integer user_a,user_k,user_s 297 parameter (user_a = 1,user_k = 3,user_s = 5) 298 DM da 299 Vec Xloc 300 PetscInt mx,xs,xe,gxs,gxe 301 PetscOffset ixx,iff 302 PetscScalar xx(0:1),ff(0:1) 303 304 PetscCall(TSGetDM(ts,da,ierr)) 305 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 306 307 ! Scatter ghost points to local vector,using the 2-step process 308 ! DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 309 ! By placing code between these two statements, computations can be 310 ! done while messages are in transition. 311 PetscCall(DMGetLocalVector(da,Xloc,ierr)) 312 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc,ierr)) 313 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc,ierr)) 314 315 ! Get access to vector data 316 PetscCall(VecGetArrayRead(Xloc,xx,ixx,ierr)) 317 PetscCall(VecGetArray(F,ff,iff,ierr)) 318 319 PetscCall(FormRHSFunctionLocal(mx,xs,xe,gxs,gxe,t,xx(ixx),ff(iff),user(user_a),user(user_k),user(user_s),ierr)) 320 321 PetscCall(VecRestoreArrayRead(Xloc,xx,ixx,ierr)) 322 PetscCall(VecRestoreArray(F,ff,iff,ierr)) 323 PetscCall(DMRestoreLocalVector(da,Xloc,ierr)) 324end subroutine FormRHSFunction 325 326! --------------------------------------------------------------------- 327! 328! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 329! 330subroutine FormIJacobian(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 331 use petscts 332 use petscdmda 333 implicit none 334 335 TS ts 336 PetscReal t,shift 337 Vec X,Xdot 338 Mat J,Jpre 339 PetscReal user(6) 340 PetscErrorCode ierr 341 integer user_a,user_k,user_s 342 parameter (user_a = 0,user_k = 2,user_s = 4) 343 344 DM da 345 PetscInt mx,xs,xe,gxs,gxe 346 PetscInt i,i1,row,col 347 PetscReal k1,k2; 348 PetscScalar val(4) 349 350 PetscCall(TSGetDM(ts,da,ierr)) 351 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 352 353 i1 = 1 354 k1 = user(user_k+1) 355 k2 = user(user_k+2) 356 do i = xs,xe 357 row = i-gxs 358 col = i-gxs 359 val(1) = shift + k1 360 val(2) = -k2 361 val(3) = -k1 362 val(4) = shift + k2 363 PetscCall(MatSetValuesBlockedLocal(Jpre,i1,row,i1,col,val,INSERT_VALUES,ierr)) 364 end do 365 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 366 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 367 if (J /= Jpre) then 368 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr)) 369 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY,ierr)) 370 end if 371end subroutine FormIJacobian 372 373subroutine FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,x,a,k,s,ierr) 374 implicit none 375 PetscInt mx,xs,xe,gxs,gxe 376 PetscScalar x(2,xs:xe) 377 PetscReal a(2),k(2),s(2) 378 PetscErrorCode ierr 379 380 PetscInt i 381 PetscReal one,hx,r,ik 382 one = 1.0 383 hx = one / mx 384 do i=xs,xe 385 r = i*hx 386 if (k(2) .ne. 0.0) then 387 ik = one/k(2) 388 else 389 ik = one 390 end if 391 x(1,i) = one + s(2)*r 392 x(2,i) = k(1)*ik*x(1,i) + s(2)*ik 393 end do 394end subroutine FormInitialSolutionLocal 395 396subroutine FormInitialSolution(ts,X,user,ierr) 397 use petscts 398 use petscdmda 399 implicit none 400 401 TS ts 402 Vec X 403 PetscReal user(6) 404 PetscErrorCode ierr 405 integer user_a,user_k,user_s 406 parameter (user_a = 1,user_k = 3,user_s = 5) 407 408 DM da 409 PetscInt mx,xs,xe,gxs,gxe 410 PetscOffset ixx 411 PetscScalar xx(0:1) 412 413 PetscCall(TSGetDM(ts,da,ierr)) 414 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 415 416 ! Get access to vector data 417 PetscCall(VecGetArray(X,xx,ixx,ierr)) 418 419 PetscCall(FormInitialSolutionLocal(mx,xs,xe,gxs,gxe,xx(ixx),user(user_a),user(user_k),user(user_s),ierr)) 420 421 PetscCall(VecRestoreArray(X,xx,ixx,ierr)) 422end subroutine FormInitialSolution 423 424! --------------------------------------------------------------------- 425! 426! IJacobian - Compute IJacobian = dF/dU + shift*dF/dUdot 427! 428subroutine FormIJacobianMF(ts,t,X,Xdot,shift,J,Jpre,user,ierr) 429 use ex22f_mfmodule 430 implicit none 431 TS ts 432 PetscReal t,shift 433 Vec X,Xdot 434 Mat J,Jpre 435 PetscReal user(6) 436 PetscErrorCode ierr 437 438 PETSC_SHIFT=shift 439 MFuser=user 440 441end subroutine FormIJacobianMF 442 443! ------------------------------------------------------------------- 444! 445! MyMult - user provided matrix multiply 446! 447! Input Parameters: 448!. X - input vector 449! 450! Output Parameter: 451!. F - function vector 452! 453subroutine MyMult(A,X,F,ierr) 454 use ex22f_mfmodule 455 implicit none 456 457 Mat A 458 Vec X,F 459 460 PetscErrorCode ierr 461 PetscScalar shift 462 463! Mat J,Jpre 464 465 PetscReal user(6) 466 467 integer user_a,user_k,user_s 468 parameter (user_a = 0,user_k = 2,user_s = 4) 469 470 DM da 471 PetscInt mx,xs,xe,gxs,gxe 472 PetscInt i,i1,row,col 473 PetscReal k1,k2; 474 PetscScalar val(4) 475 476 shift=PETSC_SHIFT 477 user=MFuser 478 479 PetscCall(TSGetDM(tscontext,da,ierr)) 480 PetscCall(GetLayout(da,mx,xs,xe,gxs,gxe,ierr)) 481 482 i1 = 1 483 k1 = user(user_k+1) 484 k2 = user(user_k+2) 485 486 do i = xs,xe 487 row = i-gxs 488 col = i-gxs 489 val(1) = shift + k1 490 val(2) = -k2 491 val(3) = -k1 492 val(4) = shift + k2 493 PetscCall(MatSetValuesBlockedLocal(Jmat,i1,row,i1,col,val,INSERT_VALUES,ierr)) 494 end do 495 496! PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 497! PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY,ierr)) 498! if (J /= Jpre) then 499 PetscCall(MatAssemblyBegin(Jmat,MAT_FINAL_ASSEMBLY,ierr)) 500 PetscCall(MatAssemblyEnd(Jmat,MAT_FINAL_ASSEMBLY,ierr)) 501! end if 502 503 PetscCall(MatMult(Jmat,X,F,ierr)) 504 505 return 506end subroutine MyMult 507 508! 509subroutine SaveSolutionToDisk(da,X,gdof,xs,xe) 510 use petscdmda 511 implicit none 512 513 Vec X 514 DM da 515 PetscInt xs,xe,two 516 PetscInt gdof,i 517 PetscErrorCode ierr 518 PetscOffset ixx 519 PetscScalar data2(2,xs:xe),data(gdof) 520 PetscScalar xx(0:1) 521 522 PetscCall(VecGetArrayRead(X,xx,ixx,ierr)) 523 524 two = 2 525 data2=reshape(xx(ixx:ixx+gdof),(/two,xe-xs+1/)) 526 data=reshape(data2,(/gdof/)) 527 open(1020,file='solution_out_ex22f_mf.txt') 528 do i=1,gdof 529 write(1020,'(e24.16,1x)') data(i) 530 end do 531 close(1020) 532 533 PetscCall(VecRestoreArrayRead(X,xx,ixx,ierr)) 534end subroutine SaveSolutionToDisk 535 536!/*TEST 537! 538! test: 539! args: -da_grid_x 200 -ts_arkimex_type 4 540! requires: !single 541! output_file: output/ex22f_mf_1.out 542! 543!TEST*/ 544