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