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