3 | | |
| 3 | This program solves either the inviscid Burgers equation or the linear advection equation with the given boundary or initial conditions. First, input is taken from a file: |
| 4 | |
| 5 | {{{#!fortran |
| 6 | NAMELIST/ProblemData/ starttime,finaltime,frames,xL,xR,cells,alpha,beta,BC,Eq,IC,advectionConst,CFL |
| 7 | |
| 8 | ! Make sure problem.data exists. If it does, read in. If not, output error and quit. |
| 9 | inquire(file='GodunovProblem.data', exist=ex) |
| 10 | if (ex) then |
| 11 | OPEN(UNIT=30, FILE='GodunovProblem.data', STATUS="OLD") |
| 12 | READ(30,NML=ProblemData) ! read into variables in namelist |
| 13 | CLOSE(30) |
| 14 | else |
| 15 | print *, 'Missing GodunovProblem.data. Exiting.' |
| 16 | stop ! exit if missing |
| 17 | end if |
| 18 | }}} |
| 19 | |
| 20 | The domain is then initialized according to the selected conditions: |
| 21 | |
| 22 | {{{#!fortran |
| 23 | select case (IC) |
| 24 | case (Gaussian) |
| 25 | do i = 1, cells |
| 26 | ProblemRegion%u(i) = alpha*exp(-1.*beta*((i-0.5)*dx+xL)**2) |
| 27 | end do |
| 28 | case (SquareWave) |
| 29 | do i = 1, cells |
| 30 | x=(i-0.5)*dx + xL |
| 31 | if (x <= 0.3) then |
| 32 | ProblemRegion%u(i) = 0 |
| 33 | else if (x >= 0.7) then |
| 34 | ProblemRegion%u(i) = 0 |
| 35 | else |
| 36 | ProblemRegion%u(i) = 1 |
| 37 | end if |
| 38 | end do |
| 39 | case (BurgersShock) |
| 40 | do i = 1, cells |
| 41 | x=(i-1)*dx + xL |
| 42 | if (x <= 0.5) then |
| 43 | ProblemRegion%u(i) = -0.5 |
| 44 | else if (x >= 1.) then |
| 45 | ProblemRegion%u(i) = 0 |
| 46 | else |
| 47 | ProblemRegion%u(i) = 1 |
| 48 | end if |
| 49 | end do |
| 50 | end select |
| 51 | }}} |
| 52 | |
| 53 | The boundary conditions are then applied: |
| 54 | |
| 55 | {{{#!fortran |
| 56 | select case (BC) |
| 57 | case(periodic) ! Wraps around |
| 58 | ProblemRegion%u(0) = ProblemRegion%u(cells) |
| 59 | ProblemRegion%u(cells+1) = ProblemRegion%u(1) |
| 60 | case(transparent) ! Wave flows out |
| 61 | ProblemRegion%u(0) = ProblemRegion%u(1) |
| 62 | ProblemRegion%u(cells+1) = ProblemRegion%u(cells) |
| 63 | end select |
| 64 | }}} |
| 65 | |
| 66 | After this, the time step is determined using the maximum speed in the domain and the given CFL coefficient. |
| 67 | |
| 68 | {{{#!fortran |
| 69 | ! Determine maximum speed in domain |
| 70 | maxSpeed = 0. |
| 71 | do i = 0, cells+1 |
| 72 | if (abs(ProblemRegion%u(i)) > maxSpeed) maxSpeed = abs(ProblemRegion%u(i)) ! (See Toro, Ch. 5.6) |
| 73 | end do |
| 74 | |
| 75 | ! Calculate delta t |
| 76 | dt = CFL*dx/maxSpeed |
| 77 | |
| 78 | ! If time step would pass the next output frame, recalculate |
| 79 | if (time + dt > nexttime) dt = nexttime - time |
| 80 | }}} |
| 81 | |
| 82 | The intercell flux is then calculated and the velocity array is updated to the next time. |
| 83 | |
| 84 | {{{#!fortran |
| 85 | do i = 0, cells |
| 86 | uL = ProblemRegion%u(i) |
| 87 | uR = ProblemRegion%u(i+1) |
| 88 | select case (Eq) |
| 89 | case (advection) |
| 90 | ProblemRegion%flux(i) = LinearAdvectionFlux(RiemannSoln(uL,uR)) |
| 91 | |
| 92 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 93 | function RiemannSoln(uL,uR) |
| 94 | |
| 95 | if (uL > uR) then |
| 96 | S = 0.5*(uL + uR) |
| 97 | if (S <= 0) then |
| 98 | RiemannSoln = uR |
| 99 | else |
| 100 | RiemannSoln = uL |
| 101 | end if |
| 102 | else |
| 103 | if (uL >= 0) then |
| 104 | RiemannSoln = uL |
| 105 | else if (uR <= 0) then |
| 106 | RiemannSoln = uR |
| 107 | else |
| 108 | RiemannSoln = 0 |
| 109 | end if |
| 110 | end if |
| 111 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 112 | |
| 113 | case (inviscidBurgers) |
| 114 | ProblemRegion%flux(i) = BurgersFlux(RiemannSoln(uL,uR)) |
| 115 | end select |
| 116 | end do |
| 117 | |
| 118 | ! Loop through cells using previously calculated flux |
| 119 | do i = 1, cells |
| 120 | ProblemRegion%u(i) = ProblemRegion%u(i) + (dt/dx)*(ProblemRegion%flux(i-1) - ProblemRegion%flux(i)) |
| 121 | end do |
| 122 | }}} |
| 123 | |
| 124 | Finally, at each requested output frame, the current system state is written to a file: |
| 125 | |
| 126 | {{{#!fortran |
| 127 | write (40,*) 'x u' ! write column headers |
| 128 | do i = 1, cells |
| 129 | x = xL + (i-0.5)*dx ! Data is located at center of each cell |
| 130 | write (40,*) x, ProblemRegion%u(i) |
| 131 | end do |
| 132 | }}} |