3 | | |
| 3 | This program solves the Euler equations in 1, 2, or 3 dimensions using Godunov's first-order upwind method. (The attachment linked in this header contains source code and makefile for all three of the programs on this page.) Input is taken from a file: |
| 4 | |
| 5 | {{{#!fortran |
| 6 | NAMELIST/ProblemData/ starttime,finaltime,frames,xL,xR,xcells,yL,yR,ycells,zL,zR,zcells,BC,IC,CFL,cov,gamma,tol,maxiter,RiemannSolver,Q,alpha,beta,advectionConst,Eq,nDim,num_threads,lthreading |
| 7 | |
| 8 | ! Make sure problem.data exists. If it does, read in. If not, output error and quit. |
| 9 | inquire(file='problem.data', exist=ex) |
| 10 | if (ex) then |
| 11 | OPEN(UNIT=30, FILE='problem.data', STATUS="OLD") |
| 12 | READ(30,NML=ProblemData) ! read into variables in namelist |
| 13 | CLOSE(30) |
| 14 | else |
| 15 | print *, 'Missing problem.data. Exiting.' |
| 16 | stop ! exit if missing |
| 17 | end if |
| 18 | }}} |
| 19 | |
| 20 | The domain is initialized based on selected test (options depend on number of dimensions): |
| 21 | |
| 22 | {{{#!fortran |
| 23 | select case (nDim) |
| 24 | |
| 25 | ! 1D case |
| 26 | case (1) |
| 27 | ! Allocate domain arrays for conserved variables and intercell flux |
| 28 | allocate(ProblemRegion%U(0:xcells+1,NrVar),ProblemRegion%flux(0:xcells,NrVar)) |
| 29 | |
| 30 | select case (IC) |
| 31 | case (1) |
| 32 | WL(irho) = 1.0 |
| 33 | WR(irho) = 0.125 |
| 34 | WL(iu) = 0.75 |
| 35 | WR(iu) = 0.0 |
| 36 | WL(ip) = 1.0 |
| 37 | WR(ip) = 0.1 |
| 38 | x0(ix) = 0.3 |
| 39 | case (2) |
| 40 | WL(irho) = 1.0 |
| 41 | WR(irho) = 1.0 |
| 42 | WL(iu) = -2.0 |
| 43 | WR(iu) = 2.0 |
| 44 | WL(ip) = 0.4 |
| 45 | WR(ip) = 0.4 |
| 46 | x0(ix) = 0.5 |
| 47 | case (3) |
| 48 | WL(irho) = 1.0 |
| 49 | WR(irho) = 1.0 |
| 50 | WL(iu) = 0.0 |
| 51 | WR(iu) = 0.0 |
| 52 | WL(ip) = 1000.0 |
| 53 | WR(ip) = 0.01 |
| 54 | x0(ix) = 0.5 |
| 55 | case (4) |
| 56 | WL(irho) = 5.99924 |
| 57 | WR(irho) = 5.99242 |
| 58 | WL(iu) = 19.975 |
| 59 | WR(iu) = -6.19633 |
| 60 | WL(ip) = 460.894 |
| 61 | WR(ip) = 46.0950 |
| 62 | x0(ix) = 0.4 |
| 63 | case (5) |
| 64 | WL(irho) = 1.0 |
| 65 | WR(irho) = 1.0 |
| 66 | WL(iu) = -19.59745 |
| 67 | WR(iu) = -19.59745 |
| 68 | WL(ip) = 1000.0 |
| 69 | WR(ip) = 0.01 |
| 70 | x0(ix) = 0.8 |
| 71 | case (6) |
| 72 | WL(irho) = 1.4 |
| 73 | WR(irho) = 1.0 |
| 74 | WL(iu) = 0.0 |
| 75 | WR(iu) = 0.0 |
| 76 | WL(ip) = 1.0 |
| 77 | WR(ip) = 1.0 |
| 78 | x0(ix) = 0.5 |
| 79 | case (7) |
| 80 | WL(irho) = 1.4 |
| 81 | WR(irho) = 1.0 |
| 82 | WL(iu) = 0.1 |
| 83 | WR(iu) = 0.1 |
| 84 | WL(ip) = 1.0 |
| 85 | WR(ip) = 1.0 |
| 86 | x0(ix) = 0.5 |
| 87 | end select |
| 88 | |
| 89 | ! Initialize domain using set initial conditions (conserved variables) |
| 90 | do i = 1, xcells |
| 91 | x(ix) = (i-0.5)*dx(ix) + xL |
| 92 | if (x(ix) <= x0(ix)) then |
| 93 | ProblemRegion%U(i,:) = PrimToCons(WL) |
| 94 | |
| 95 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 96 | PrimToCons(irho) = W(irho) |
| 97 | PrimToCons(imom(ix)) = W(iu(ix))*W(irho) |
| 98 | PrimToCons(iE) = W(ip)*(1. - cov*W(irho))/(gamma - 1.) + 0.5*W(irho)*sum(W(iu)**2) |
| 99 | ! Tangential directins |
| 100 | if (nDim > 1) PrimToCons(imom(iy)) = W(iu(iy))*W(irho) |
| 101 | if (nDim > 2) PrimToCons(imom(iz)) = W(iu(iz))*W(irho) |
| 102 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 103 | |
| 104 | else |
| 105 | ProblemRegion%U(i,:) = PrimToCons(WR) |
| 106 | end if |
| 107 | end do |
| 108 | |
| 109 | ! 2D case |
| 110 | case (2) |
| 111 | allocate(ProblemRegion2D%U(0:xcells+1,0:ycells+1,NrVar),ProblemRegion%U(0:max(xcells,ycells)+1,NrVar),ProblemRegion%flux(0:max(xcells,ycells),NrVar)) |
| 112 | |
| 113 | select case (IC) |
| 114 | case (explosion) |
| 115 | ! "Left" is inside circular boundary, "Right" is outside circular boundary |
| 116 | WL(irho) = 1.0 |
| 117 | WL(iu(ix)) = 0.0 |
| 118 | WL(iu(iy)) = 0.0 |
| 119 | WL(ip) = 1.0 |
| 120 | WR(irho) = 0.125 |
| 121 | WR(iu(ix)) = 0.0 |
| 122 | WR(iu(iy)) = 0.0 |
| 123 | WR(ip) = 0.1 |
| 124 | x0 = (/ 0, 0 /) |
| 125 | radius = 0.4 |
| 126 | do i = 1, xcells |
| 127 | x(ix) = (i-0.5)*dx(ix) + xL |
| 128 | do j = 1, ycells |
| 129 | x(iy) = (j-0.5)*dx(iy) + yL |
| 130 | if (sum((x-x0)**2) < radius**2) then |
| 131 | ProblemRegion2D%U(i,j,:) = PrimToCons(WL) |
| 132 | print *, ProblemRegion2D%U(i,j,:) |
| 133 | else |
| 134 | ProblemRegion2D%U(i,j,:) = PrimToCons(WR) |
| 135 | print *, ProblemRegion2D%U(i,j,:) |
| 136 | end if |
| 137 | end do |
| 138 | end do |
| 139 | end select |
| 140 | |
| 141 | ! 3D case |
| 142 | case (3) |
| 143 | allocate(ProblemRegion3D%U(0:xcells+1,0:ycells+1,0:zcells+1,NrVar),ProblemRegion%U(0:max(xcells,ycells,zcells)+1,NrVar),ProblemRegion%flux(0:max(xcells,ycells,zcells),NrVar)) |
| 144 | |
| 145 | select case (IC) |
| 146 | case (explosion) |
| 147 | ! "Left" is inside spherical boundary, "Right" is outside spherical boundary |
| 148 | WL(irho) = 1.0 |
| 149 | WL(iu(ix)) = 0.0 |
| 150 | WL(iu(iy)) = 0.0 |
| 151 | WL(iu(iz)) = 0.0 |
| 152 | WL(ip) = 1.0 |
| 153 | WR(irho) = 0.125 |
| 154 | WR(iu(ix)) = 0.0 |
| 155 | WR(iu(iy)) = 0.0 |
| 156 | WR(iu(iz)) = 0.0 |
| 157 | WR(ip) = 0.1 |
| 158 | x0 = (/ 0, 0, 0 /) |
| 159 | radius = 0.4 |
| 160 | do i = 1, xcells |
| 161 | x(ix) = (i-0.5)*dx(ix) + xL |
| 162 | do j = 1, ycells |
| 163 | x(iy) = (j-0.5)*dx(iy) + yL |
| 164 | do k = 1, zcells |
| 165 | x(iz) = (k-0.5)*dx(iz) + zL |
| 166 | if (sum((x-x0)**2) < radius**2) then |
| 167 | ProblemRegion3D%U(i,j,k,:) = PrimToCons(WL) |
| 168 | else |
| 169 | ProblemRegion3D%U(i,j,k,:) = PrimToCons(WR) |
| 170 | end if |
| 171 | end do |
| 172 | end do |
| 173 | end do |
| 174 | end select |
| 175 | |
| 176 | end select |
| 177 | }}} |
| 178 | |
| 179 | After the domain is initialized, the first frame is written to file by |
| 180 | |
| 181 | {{{#!fortran |
| 182 | ! Euler solver with Godunov method |
| 183 | case (GodunovEuler) |
| 184 | if (nDim == 1) then |
| 185 | write (40,*) 'x rho u p' ! write column headers |
| 186 | do i = 1, xcells |
| 187 | x(ix) = xL + (i-0.5)*dx(ix) ! Data is located at center of each cell |
| 188 | U = ProblemRegion%U(i,:) |
| 189 | write (40,*) x, ConsToPrim(U) |
| 190 | |
| 191 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 192 | ConsToPrim(irho) = U(irho) |
| 193 | ConsToPrim(iu(ix)) = U(imom(ix))/U(irho) |
| 194 | ConsToPrim(ip) = (U(iE) - (sum(U(imom)**2))/(2*U(irho)))*(gamma - 1.)/(1. - cov*U(irho)) |
| 195 | ! Tangential directions |
| 196 | if (nDim > 1) ConsToPrim(iu(iy)) = U(imom(iy))/U(irho) |
| 197 | if (nDim > 2) ConsToPrim(iu(iz)) = U(imom(iz))/U(irho) |
| 198 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 199 | |
| 200 | end do |
| 201 | else if (nDim == 2) then |
| 202 | write (40,*) 'x y rho u v p' |
| 203 | do i = 1, xcells |
| 204 | x(ix) = xL + (i-0.5)*dx(ix) ! Data is located at center of each cell |
| 205 | do j = 1, ycells |
| 206 | x(iy) = yL + (j-0.5)*dx(iy) ! Data is located at center of each cell |
| 207 | U = ProblemRegion2D%U(i,j,:) |
| 208 | write (40,*) x, ConsToPrim(U) |
| 209 | end do |
| 210 | end do |
| 211 | else if (nDim == 3) then |
| 212 | write (40,*) 'x y z rho u v w p' |
| 213 | do i = 1, xcells |
| 214 | x(1) = xL + (i-0.5)*dx(1) ! Data is located at center of each cell |
| 215 | do j = 1, ycells |
| 216 | x(2) = yL + (j-0.5)*dx(2) ! Data is located at center of each cell |
| 217 | do k = 1, zcells |
| 218 | x(3) = zL + (k-0.5)*dx(3) ! Data is located at center of each cell |
| 219 | U = ProblemRegion3D%U(i,j,k,:) |
| 220 | write (40,*) x, ConsToPrim(U) |
| 221 | end do |
| 222 | end do |
| 223 | end do |
| 224 | end if |
| 225 | }}} |
| 226 | |
| 227 | The time-stepping procedure is then entered. The $\Delta t$ for the current step is calculated by determining the maximum speed in any direction on the grid: |
| 228 | |
| 229 | {{{#!fortran |
| 230 | if (nDim == 1) then |
| 231 | do i = 1, xcells |
| 232 | U = ProblemRegion%U(i,:) |
| 233 | cell = region(ConsToPrim(U)) |
| 234 | if (cell%c .ne. cell%c) then |
| 235 | print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.' |
| 236 | stop |
| 237 | end if |
| 238 | if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c |
| 239 | end do |
| 240 | ! 2D |
| 241 | else if (nDim == 2) then |
| 242 | do i = 1, xcells |
| 243 | do j = 1, ycells |
| 244 | U = ProblemRegion2D%U(i,j,:) |
| 245 | cell = region(ConsToPrim(U)) |
| 246 | if (cell%c .ne. cell%c) then |
| 247 | print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.' |
| 248 | stop |
| 249 | end if |
| 250 | if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c |
| 251 | if (abs(cell%W(iu(iy))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iy))) + cell%c |
| 252 | end do |
| 253 | end do |
| 254 | ! 3D |
| 255 | else if (nDim == 3) then |
| 256 | do i = 1, xcells |
| 257 | do j = 1, ycells |
| 258 | do k = 1, zcells |
| 259 | U = ProblemRegion3D%U(i,j,k,:) |
| 260 | cell = region(ConsToPrim(U)) |
| 261 | if (cell%c .ne. cell%c) then |
| 262 | print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.' |
| 263 | stop |
| 264 | end if |
| 265 | if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c |
| 266 | if (abs(cell%W(iu(iy))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iy))) + cell%c |
| 267 | if (abs(cell%W(iu(iz))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iz))) + cell%c |
| 268 | end do |
| 269 | end do |
| 270 | end do |
| 271 | end if |
| 272 | }}} |
| 273 | |
| 274 | The time-stepping procedure then proceeds in a slightly different manner, depending upon the number of dimensions in the problem. In the 1D case, boundary conditions are applied: |
| 275 | |
| 276 | {{{#!fortran |
| 277 | select case (BC) |
| 278 | case(periodic) ! Wraps around |
| 279 | ProblemRegion%U(0,1:NrVar) = ProblemRegion%U(cells,1:NrVar) |
| 280 | ProblemRegion%U(cells+1,1:NrVar) = ProblemRegion%U(1,1:NrVar) |
| 281 | case(transparent) ! Wave flows out |
| 282 | ProblemRegion%U(0,1:NrVar) = ProblemRegion%U(1,1:NrVar) |
| 283 | ProblemRegion%U(cells+1,1:NrVar) = ProblemRegion%U(cells,1:NrVar) |
| 284 | case(reflective) ! Wave bounces off hard boundary |
| 285 | ProblemRegion%U(cells+1,irho) = ProblemRegion%U(cells,irho) |
| 286 | ProblemRegion%U(cells+1,imom) = ProblemRegion%U(cells,imom) |
| 287 | ProblemRegion%U(cells+1,imom(ix)) = -ProblemRegion%U(cells,imom(ix)) ! Momentum normal to boundary is reversed |
| 288 | ProblemRegion%U(cells+1,iE) = ProblemRegion%U(cells,iE) |
| 289 | ProblemRegion%U(0,irho) = ProblemRegion%U(1,irho) |
| 290 | ProblemRegion%U(0,imom) = ProblemRegion%U(1,imom) |
| 291 | ProblemRegion%U(0,imom(ix)) = -ProblemRegion%U(1,imom(ix)) |
| 292 | ProblemRegion%U(0,iE) = ProblemRegion%U(1,iE) |
| 293 | end select |
| 294 | }}} |
| 295 | |
| 296 | and the flux across cell boundaries is then calculated. Depending on the method selected, the flux may be estimated directly, or the fluid state may be approximated, which can then be used to calculate the flux. Summaries of the various methods here can be found near the end of this page. |
| 297 | |
| 298 | {{{#!fortran |
| 299 | ! Loop through domain using equation set in problem.data |
| 300 | do i = 0, cells |
| 301 | |
| 302 | U = ProblemRegion%U(i,:) |
| 303 | L = region(ConsToPrim(U)) |
| 304 | if (L%c .ne. L%c) then |
| 305 | print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in CalcFlux.' |
| 306 | stop |
| 307 | end if |
| 308 | U = ProblemRegion%U(i+1,:) |
| 309 | R = region(ConsToPrim(U)) |
| 310 | |
| 311 | ! Some solvers approximate flux rather than fluid state |
| 312 | select case (RiemannSolver) |
| 313 | case (HLLC) |
| 314 | call WaveSpeed(L,R,SL,SR,S0) |
| 315 | |
| 316 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 317 | p0 = max(0d0, PrimitiveVariablePressure(L,R)) |
| 318 | ! Approximate shock or rarefaction speeds |
| 319 | if (p0 > L%W(ip)) then |
| 320 | SL = L%W(iu(ix)) - L%c*sqrt(1 + (gamma + 1.)*(p0/L%W(ip) - 1.)/(2.*gamma)) |
| 321 | else |
| 322 | SL = L%W(iu(ix)) - L%c |
| 323 | end if |
| 324 | if (p0 > R%W(ip)) then |
| 325 | SR = R%W(iu(ix)) + R%c*sqrt(1 + (gamma + 1.)*(p0/R%W(ip) - 1.)/(2.*gamma)) |
| 326 | else |
| 327 | SR = R%W(iu(ix)) + R%c |
| 328 | end if |
| 329 | if (RiemannSolver == HLLC) S0 = (R%W(ip) - L%W(ip) + L%W(irho)*L%W(iu(ix))*(SL - L%W(iu(ix))) - R%W(irho)*R%W(iu(ix))*(SR - R%W(iu(ix))))/(L%W(irho)*(SL - L%W(iu(ix))) - R%W(irho)*(SR - R%W(iu(ix)))) |
| 330 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 331 | |
| 332 | ProblemRegion%flux(i,:) = HLLCFlux(L,R,S0,SL,SR) |
| 333 | |
| 334 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 335 | ! Approximate with flux from left Riemann state |
| 336 | if (0 <= SL) then |
| 337 | HLLCflux = EulerFlux(L%W) |
| 338 | ! Approximate with integral average of left star region |
| 339 | else if (SL <= 0 .and. 0 <= S0) then |
| 340 | UL = PrimToCons(L%W) |
| 341 | U0L(irho) = L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0) |
| 342 | U0L(imom(ix)) = S0*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0) |
| 343 | if (nDim > 1) U0L(imom(iy)) = L%W(iu(iy))*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0) |
| 344 | if (nDim > 2) U0L(imom(iz)) = L%W(iu(iz))*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0) |
| 345 | U0L(iE) = L%W(irho)*((SL - L%W(iu(ix)))/(SL - S0))*(UL(iE)/L%W(irho) + (S0 - L%W(iu(ix)))*(S0 + L%W(ip)/(L%W(irho)*(SL - L%W(iu(ix)))))) |
| 346 | HLLCflux = EulerFlux(L%W) + SL*(U0L - UL) |
| 347 | ! Approximate with integral average of right star region |
| 348 | else if (S0 <= 0 .and. 0 <= SR) then |
| 349 | UR = PrimToCons(R%W) |
| 350 | U0R(irho) = R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0) |
| 351 | U0R(imom(ix)) = S0*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0) |
| 352 | if (nDim > 1) U0R(imom(iy)) = R%W(iu(iy))*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0) |
| 353 | if (nDim > 2) U0R(imom(iz)) = R%W(iu(iz))*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0) |
| 354 | U0R(iE) = R%W(irho)*((SR - R%W(iu(ix)))/(SR - S0))*(UR(iE)/R%W(irho) + (S0 - R%W(iu(ix)))*(S0 + R%W(ip)/(R%W(irho)*(SR - R%W(iu(ix)))))) |
| 355 | HLLCflux = EulerFlux(R%W) + SR*(U0R - UR) |
| 356 | ! Approximate with right Riemann state |
| 357 | else |
| 358 | HLLCflux = EulerFlux(R%W) |
| 359 | end if |
| 360 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 361 | |
| 362 | case (HLL) |
| 363 | call WaveSpeed(L,R,SL,SR) |
| 364 | ProblemRegion%flux(i,:) = HLLFlux(L,R,SL,SR) |
| 365 | |
| 366 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 367 | ! Approximate with left Riemann state |
| 368 | if (0 <= SL) then |
| 369 | HLLflux = EulerFlux(L%W) |
| 370 | ! Approximate with integral average of star region |
| 371 | else if (SL <= 0 .and. 0 <= SR) then |
| 372 | UL = PrimToCons(L%W) |
| 373 | UR = PrimToCons(R%W) |
| 374 | HLLflux = (SR*EulerFlux(L%W) - SL*EulerFlux(R%W) + SL*SR*(UR - UL))/(SR - SL) |
| 375 | ! Approximate with right Riemann state |
| 376 | else |
| 377 | HLLflux = EulerFlux(R%W) |
| 378 | end if |
| 379 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 380 | |
| 381 | ! If solver approximates fluid state, use it to calculate flux |
| 382 | case default |
| 383 | ProblemRegion%flux(i,:) = EulerFlux(RiemannSoln(L,R)) |
| 384 | |
| 385 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 386 | function EulerRiemannSoln |
| 387 | |
| 388 | ! Check for vacuum |
| 389 | call check_vacuum(L,R) |
| 390 | |
| 391 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 392 | critdiff = 2.*L%c/(gamma - 1.) + 2.*R%c/(gamma - 1.) ! Difference between right and left velocities must be strictly less than this for pressure to remain positive |
| 393 | |
| 394 | if ((L%W(ip) == 0 .and. L%W(irho) == 0) .or. (R%W(ip) == 0 .and. R%w(irho) == 0)) then ! one of the initial states is vacuum |
| 395 | vacuum = 1 |
| 396 | print *, 'One of initial states is vacuum.' |
| 397 | else if (critdiff <= R%W(iu(ix)) - L%W(iu(ix))) then ! vacuum is created by initial states |
| 398 | vacuum = 2 |
| 399 | print *, 'Vacuum is created by initial states.' |
| 400 | else ! no vacuum |
| 401 | vacuum = 3 |
| 402 | end if |
| 403 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 404 | |
| 405 | ! Vacuum present, don't need star region - exact solution is cheap |
| 406 | if (vacuum /= 3) then |
| 407 | EulerRiemannSoln = sample(L,R,0d0) |
| 408 | |
| 409 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 410 | ! Calculate shock and rarefaction speeds |
| 411 | SvacL = L%W(iu(ix)) + 2.*L%c/(gamma - 1.) ! Speed of vacuum front (vacuum on right) |
| 412 | SvacR = R%W(iu(ix)) - 2.*R%c/(gamma - 1.) ! Speed of vacuum front (vacuum on left) |
| 413 | |
| 414 | ! Initial state is vacuum |
| 415 | if (vacuum == 1) then |
| 416 | if (R%W(ip) == 0 .and. R%W(irho) == 0) then !vacuum on the right |
| 417 | lright = .false. |
| 418 | R%W(iu(ix)) = 0 |
| 419 | if (location <= L%W(iu(ix)) - L%c) then ! Left of rarefaction fan |
| 420 | sampleWithVacuum = L%W |
| 421 | else if (location >= SvacL) then! Right of rarefaction fan (in vacuum) |
| 422 | sampleWithVacuum = R%W |
| 423 | else ! In rarefaction fan |
| 424 | call rarefaction(L,R,Wfan,lright,location) |
| 425 | |
| 426 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 427 | ! If fluid has covolume constant, iterate to solve for density |
| 428 | if (cov > 0) then |
| 429 | if (.not. lright) then |
| 430 | beta = ((((gamma - 1.)*(L%W(iu(ix)) + 2.*L%c*(1. - cov*L%W(irho))/(gamma - 1.) - location))**2.)/(gamma*L%W(ip)*((1. - cov*L%W(irho))/L%W(irho))**gamma)) ! Constant related to Riemann invariant |
| 431 | Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf,L,R) |
| 432 | Wfan(ip) = L%W(ip)*(((1. - cov*L%W(irho))/L%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma) |
| 433 | c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho)))) ! Speed of sound at current location |
| 434 | Wfan(iu(ix)) = location + c |
| 435 | if (nDim > 1) Wfan(iu(iy)) = L%W(iu(iy)) |
| 436 | if (nDim > 2) Wfan(iu(iz)) = L%W(iu(iz)) |
| 437 | else |
| 438 | beta = (((location - R%W(iu(ix)) + (2.*R%c*(1. - cov*R%W(irho))/(gamma - 1.)))**2.)/(gamma*R%W(ip)*((1. - cov*R%W(irho))/R%W(irho))**gamma)) ! Constant related to Riemann invariant |
| 439 | Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf,L,R) |
| 440 | Wfan(ip) = R%W(ip)*(((1. - cov*R%W(irho))/R%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma) |
| 441 | c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho)))) ! Speed of sound at current location |
| 442 | Wfan(iu(ix)) = location - c |
| 443 | if (nDim > 1) Wfan(iu(iy)) = R%W(iu(iy)) |
| 444 | if (nDim > 2) Wfan(iu(iz)) = R%W(iu(iz)) |
| 445 | end if |
| 446 | ! Otherwise solution is closed-form |
| 447 | else |
| 448 | ! Rarefaction fan differs depending on whether it's a right or left rarefaction |
| 449 | if (.not. lright) then |
| 450 | Wfan(irho) = L%W(irho)*(2/(gamma + 1) + (gamma - 1)*(L%W(iu(ix)) - location)/((gamma + 1)*L%c))**(2/(gamma - 1)) |
| 451 | Wfan(iu(ix)) = (2/(gamma + 1))*(L%c + (gamma - 1)*L%W(iu(ix))/2 + location) |
| 452 | Wfan(ip) = L%W(ip)*(2/(gamma + 1) + (gamma - 1)*(L%W(iu(ix)) - location)/((gamma + 1)*L%c))**(2*gamma/(gamma - 1)) |
| 453 | ! Tangential directions unchanged |
| 454 | if (nDim > 1) Wfan(iu(iy)) = L%W(iu(iy)) |
| 455 | if (nDim > 2) Wfan(iu(iz)) = L%W(iu(iz)) |
| 456 | else |
| 457 | Wfan(irho) = R%W(irho)*(2/(gamma + 1) - (gamma - 1)*(R%W(iu(ix)) - location)/((gamma + 1)*R%c))**(2/(gamma - 1)) |
| 458 | Wfan(iu(ix)) = (2/(gamma + 1))*(-R%c + (gamma - 1)*R%W(iu(ix))/2 + location) |
| 459 | Wfan(ip) = R%W(ip)*(2/(gamma + 1) - (gamma - 1)*(R%W(iu(ix)) - location)/((gamma + 1)*R%c))**(2*gamma/(gamma - 1)) |
| 460 | ! Tangential directions unchanged |
| 461 | if (nDim > 1) Wfan(iu(iy)) = R%W(iu(iy)) |
| 462 | if (nDim > 2) Wfan(iu(iz)) = R%W(iu(iz)) |
| 463 | end if |
| 464 | end if |
| 465 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 466 | |
| 467 | sampleWithVacuum = Wfan |
| 468 | end if |
| 469 | else!vacuum on the left |
| 470 | lright = .true. |
| 471 | L%W(iu(ix)) = 0 |
| 472 | if (location >= R%W(iu(ix)) - R%c) then ! Right of rarefaction fan |
| 473 | sampleWithVacuum = R%W |
| 474 | else if (location <= SvacR) then! Left of rarefaction fan (in vacuum) |
| 475 | sampleWithVacuum = L%W |
| 476 | else |
| 477 | call rarefaction(L,R,Wfan,lright,location) ! In rarefaction fan |
| 478 | sampleWithVacuum = Wfan |
| 479 | end if |
| 480 | end if |
| 481 | ! Initial regions result in vacuum between |
| 482 | else if (vacuum == 2) then |
| 483 | if (location <= L%W(iu(ix)) - L%c) then ! Left of left rarefaction fan |
| 484 | sampleWithVacuum = L%W |
| 485 | else if (location <= SvacL) then ! In left rarefaction fan |
| 486 | lright = .false. |
| 487 | call rarefaction(L,R,Wfan,lright,location) |
| 488 | sampleWithVacuum = Wfan |
| 489 | else if (location >= R%W(iu(ix)) + R%c) then ! Right of right rarefaction fan |
| 490 | sampleWithVacuum = R%W |
| 491 | else if (location >= SvacR) then ! In right rarefaction fan |
| 492 | lright = .true. |
| 493 | call rarefaction(L,R,Wfan,lright,location) |
| 494 | sampleWithVacuum = Wfan |
| 495 | else! In vacuum (between rarefaction fans) |
| 496 | sampleWithVacuum(irho) = 0 |
| 497 | sampleWithVacuum(iu(ix)) = 0 |
| 498 | sampleWithVacuum(ip) = 0 |
| 499 | if (nDim > 1) sampleWithVacuum(iu(iy)) = 0 |
| 500 | if (nDim > 2) sampleWithVacuum(iu(iz)) = 0 |
| 501 | endif |
| 502 | end if |
| 503 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 504 | |
| 505 | ! No vacuum, calculate star region and sample along t axis |
| 506 | else |
| 507 | ! Star region may be approximated |
| 508 | call StarRegion(L,R,Lstar,Rstar) |
| 509 | |
| 510 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 511 | ! Passively advected - no math, so no approximations |
| 512 | if (nDim > 1) then |
| 513 | Lstar%W(iu(iy)) = L%W(iu(iy)) |
| 514 | Rstar%W(iu(iy)) = R%W(iu(iy)) |
| 515 | end if |
| 516 | if (nDim > 2) then |
| 517 | Lstar%W(iu(iz)) = L%W(iu(iz)) |
| 518 | Rstar%W(iu(iz)) = R%W(iu(iz)) |
| 519 | end if |
| 520 | |
| 521 | select case (RiemannSolver) |
| 522 | ! Exact Riemann solver |
| 523 | case (exact) |
| 524 | ! Use Newton-Raphson solver to find pressure in star region |
| 525 | Lstar%W(ip) = NRSolver(pguess(L,R),covolumeF,covolumeDf,L,R) |
| 526 | ! Calculate velocity and densities in star regions |
| 527 | Lstar%W(iu(ix)) = exactVelocity(Lstar%W(ip),L,R) |
| 528 | Rstar%W(ip) = Lstar%W(ip) |
| 529 | Rstar%W(iu(ix)) = Rstar%W(iu(ix)) |
| 530 | Lstar%W(irho) = exactDensity(Lstar,L) |
| 531 | Rstar%W(irho) = exactDensity(Rstar,R) |
| 532 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 533 | ! See below for exact Riemann solver summary ! |
| 534 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 535 | ! Two-shock approximation |
| 536 | case (TSRS) |
| 537 | p0 = max(PrimitiveVariablePressure(L,R),0d0) |
| 538 | Lstar%W(ip) = TwoShockPressure(L,R,p0) |
| 539 | Lstar%W(iu(ix)) = TwoShockVelocity(L,R,Lstar,p0) |
| 540 | Rstar%W(ip) = Lstar%W(ip) |
| 541 | Rstar%W(iu(ix)) = Rstar%W(iu(ix)) |
| 542 | Lstar%W(irho) = TwoShockDensity(L,Lstar) |
| 543 | Rstar%W(irho) = TwoShockDensity(R,Rstar) |
| 544 | ! Adaptive non-iterative method |
| 545 | case (ANRS) |
| 546 | p0 = PrimitiveVariablePressure(L,R) |
| 547 | pmax = max(L%W(ip),R%W(ip)) |
| 548 | pmin = min(L%W(ip),R%W(ip)) |
| 549 | if (pmax/pmin < Q .and. pmin < p0 .and. pmax > p0) then |
| 550 | Lstar%W(ip) = p0 |
| 551 | Lstar%W(iu(ix)) = PrimitiveVariableVelocity(L,R) |
| 552 | Rstar%W(ip) = Lstar%W(ip) |
| 553 | Rstar%W(iu(ix)) = Rstar%W(iu(ix)) |
| 554 | Lstar%W(irho) = PrimitiveVariableDensity(L,Lstar,R) |
| 555 | Rstar%W(irho) = PrimitiveVariableDensity(R,Rstar,L) |
| 556 | else if (p0 < pmin) then |
| 557 | Lstar%W(ip) = TwoRarefactionPressure(L,R) |
| 558 | Lstar%W(iu(ix)) = TwoRarefactionVelocity(L,R) |
| 559 | Rstar%W(ip) = Lstar%W(ip) |
| 560 | Rstar%W(iu(ix)) = Rstar%W(iu(ix)) |
| 561 | Lstar%W(irho) = TwoRarefactionDensity(L,Lstar) |
| 562 | Rstar%W(irho) = TwoRarefactionDensity(R,Rstar) |
| 563 | else |
| 564 | Lstar%W(ip) = TwoShockPressure(L,R,p0) |
| 565 | Lstar%W(iu(ix)) = TwoShockVelocity(L,R,Lstar,p0) |
| 566 | Rstar%W(ip) = Lstar%W(ip) |
| 567 | Rstar%W(iu(ix)) = Rstar%W(iu(ix)) |
| 568 | Lstar%W(irho) = TwoShockDensity(L,Lstar) |
| 569 | Rstar%W(irho) = TwoShockDensity(R,Rstar) |
| 570 | end if |
| 571 | end select |
| 572 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 573 | |
| 574 | EulerRiemannSoln = sample(L,R,Lstar,Rstar,0d0) |
| 575 | |
| 576 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 577 | ! Calculate shock and rarefaction speeds |
| 578 | SL = L%W(iu(ix)) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma)) ! Left shock speed |
| 579 | SHL = L%W(iu(ix)) - L%c ! Head speed of left rarefaction |
| 580 | STL = Lstar%W(iu(ix)) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma)) ! Tail speed of left rarefaction |
| 581 | SR = R%W(iu(ix)) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma)) ! right |
| 582 | SHR = R%W(iu(ix)) + R%c ! right |
| 583 | STR = Rstar%W(iu(ix)) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma)) ! right |
| 584 | |
| 585 | if (location <= Lstar%W(iu(ix))) then !left side of contact |
| 586 | lright = .false. |
| 587 | if (Lstar%W(ip) > L%W(ip)) then !shock on left side |
| 588 | if (location <= SL) then !left side of shock |
| 589 | sampleWithoutVacuum = L%W |
| 590 | else |
| 591 | sampleWithoutVacuum = Lstar%W |
| 592 | end if |
| 593 | else !rarefaction on left side |
| 594 | if (location <= SHL) then !left of the rarefaction fan |
| 595 | sampleWithoutVacuum = L%W |
| 596 | else if (location >= STL) then !right of the rarefaction fan |
| 597 | sampleWithoutVacuum = Lstar%W |
| 598 | else !in rarefaction fan |
| 599 | call rarefaction(L,R,Wfan,lright,location) |
| 600 | sampleWithoutVacuum = Wfan |
| 601 | end if |
| 602 | end if |
| 603 | else !right side of contact |
| 604 | lright = .true. |
| 605 | if (Rstar%W(ip) > R%W(ip)) then !shock on right side |
| 606 | if (location >= SR) then !right side of shock |
| 607 | sampleWithoutVacuum = R%W |
| 608 | else |
| 609 | sampleWithoutVacuum = Rstar%W |
| 610 | end if |
| 611 | else !rarefaction on right side |
| 612 | if (location >= SHR) then !right of the rarefaction fan |
| 613 | sampleWithoutVacuum = R%W |
| 614 | else if (location <= STR) then !left of the rarefaction fan |
| 615 | sampleWithoutVacuum = Rstar%W |
| 616 | else !in rarefaction fan |
| 617 | call rarefaction(L,R,Wfan,lright,location) |
| 618 | sampleWithoutVacuum = Wfan |
| 619 | end if |
| 620 | end if |
| 621 | end if |
| 622 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 623 | |
| 624 | deallocate(Lstar%W,Rstar%W) |
| 625 | end if |
| 626 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 627 | !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ |
| 628 | EulerFlux(irho) = W(irho)*W(iu(ix)) |
| 629 | EulerFlux(imom(ix)) = W(irho)*W(iu(ix))**2 + W(ip) |
| 630 | EulerFlux(iE) = W(iu(ix))*((W(ip)*(1. - cov*W(ip)))/(gamma - 1.) + 0.5*W(irho)*W(iu(ix))**2 + W(ip)) |
| 631 | if (nDim > 1) EulerFlux(imom(iy)) = W(irho)*W(iu(ix))*W(iu(iy)) |
| 632 | if (nDim > 2) EulerFlux(imom(iz)) = W(irho)*W(iu(ix))*W(iu(iz)) |
| 633 | !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 634 | |
| 635 | end select |
| 636 | end do |
| 637 | }}} |
| 638 | |
| 639 | The grid is then updated using Godunov's method to find the new state of the system: |
| 640 | |
| 641 | {{{#!fortran |
| 642 | do i = 1, cells |
| 643 | ProblemRegion%U(i,:) = ProblemRegion%U(i,:) + (dt/dl)*(ProblemRegion%flux(i-1,:) - ProblemRegion%flux(i,:)) |
| 644 | if (any(ProblemRegion%U(i,:) .ne. ProblemRegion%U(i,:))) then |
| 645 | print *,' NaN at time ', time, 'in cell ',i, ' in UpdateCells.' |
| 646 | stop |
| 647 | end if |
| 648 | end do |
| 649 | }}} |
| 650 | |
| 651 | The time is updated, and if the program is at one of the times requested for output, information is written to file (as above). In 2D, the problem is first reduced to a 1D problem along each strip of cells in the x direction: |
| 652 | |
| 653 | {{{#!fortran |
| 654 | do j = 1, ycells |
| 655 | ProblemRegion%U(:,:) = ProblemRegion2D%U(:,j,:) |
| 656 | }}} |
| 657 | |
| 658 | The solution then proceeds as above, and the current strip is updated with the 1D solution: |
| 659 | |
| 660 | {{{#!fortran |
| 661 | ! Update 2D grid from 1D solution |
| 662 | ProblemRegion2D%U(:,j,:) = ProblemRegion%U(:,:) |
| 663 | end do |
| 664 | }}} |
| 665 | |
| 666 | Once the solution in the x direction has been obtained, the problem is reduced to 1D along each strip in the y direction: |
| 667 | |
| 668 | {{{#!fortran |
| 669 | do i = 1, xcells |
| 670 | ProblemRegion%U(:,:) = ProblemRegion2D%U(i,:,:) |
| 671 | ! x is considered the direction normal to the cells in 1D routines - so switch v to u for y sweep |
| 672 | ProblemRegion%U(:,imom(ix)) = ProblemRegion2D%U(i,:,imom(iy)) |
| 673 | ProblemRegion%U(:,imom(iy)) = ProblemRegion2D%U(i,:,imom(ix)) |
| 674 | }}} |
| 675 | |
| 676 | The solution proceeds once again as the 1D problem, and the grid is updated with the new solution in the y direction. |
| 677 | |
| 678 | {{{#!fortran |
| 679 | ! Update 2D grid |
| 680 | ProblemRegion2D%U(i,:,:) = ProblemRegion%U(:,:) |
| 681 | ! Switch speeds back |
| 682 | ProblemRegion2D%U(i,:,imom(iy)) = ProblemRegion%U(:,imom(ix)) |
| 683 | ProblemRegion2D%U(i,:,imom(ix)) = ProblemRegion%U(:,imom(iy)) |
| 684 | end do |
| 685 | }}} |
| 686 | |
| 687 | Once complete, the time is advanced and output may once again occur. In 3D, the procedure is the same as in 2D, incorporating the z direction as appropriate: |
| 688 | |
| 689 | {{{#!fortran |
| 690 | ! For each strip in x direction, copy to 1D region, then solve |
| 691 | do j = 1, ycells |
| 692 | do k = 1, zcells |
| 693 | ProblemRegion%U(:,:) = ProblemRegion3D%U(:,j,k,:) |
| 694 | |
| 695 | !!!!!!!!!!!!!!!!!!!!!!! |
| 696 | ! As before ! |
| 697 | !!!!!!!!!!!!!!!!!!!!!!! |
| 698 | |
| 699 | ! Update 3D grid |
| 700 | ProblemRegion3D%U(:,j,k,:) = ProblemRegion%U(:,:) |
| 701 | end do |
| 702 | end do |
| 703 | |
| 704 | ! For each strip in y direction, copy to 1D region, then solve |
| 705 | do i = 1, xcells |
| 706 | do k = 1, zcells |
| 707 | ProblemRegion%U(:,:) = ProblemRegion3D%U(i,:,k,:) |
| 708 | ! x is considered the direction normal to the cells in 1D routines - so switch v to u for y sweep (w doesn't change - still passively advected) |
| 709 | ProblemRegion%U(:,imom(ix)) = ProblemRegion3D%U(i,:,k,imom(iy)) |
| 710 | ProblemRegion%U(:,imom(iy)) = ProblemRegion3D%U(i,:,k,imom(ix)) |
| 711 | |
| 712 | !!!!!!!!!!!!!!!!!!!!!!! |
| 713 | ! As before ! |
| 714 | !!!!!!!!!!!!!!!!!!!!!!! |
| 715 | |
| 716 | ! Update 3D grid |
| 717 | ProblemRegion3D%U(i,:,k,:) = ProblemRegion%U(:,:) |
| 718 | ! switch speeds back |
| 719 | ProblemRegion3D%U(i,:,k,imom(ix)) = ProblemRegion%U(:,imom(iy)) |
| 720 | ProblemRegion3D%U(i,:,k,imom(iy)) = ProblemRegion%U(:,imom(ix)) |
| 721 | end do |
| 722 | end do |
| 723 | |
| 724 | ! For each strip in z direction, copy to 1D region, then solve |
| 725 | do i = 1, xcells |
| 726 | do j = 1, ycells |
| 727 | |
| 728 | ProblemRegion%U(:,:) = ProblemRegion3D%U(i,j,:,:) |
| 729 | ! x is considered the direction normal to the cells for 1D routines - so switch w to u for y sweep (v doesn't change - still passively advected) |
| 730 | ProblemRegion%U(:,imom(ix)) = ProblemRegion3D%U(i,j,:,imom(iz)) |
| 731 | ProblemRegion%U(:,imom(iz)) = ProblemRegion3D%U(i,j,:,imom(ix)) |
| 732 | |
| 733 | !!!!!!!!!!!!!!!!!!!!!!! |
| 734 | ! As before ! |
| 735 | !!!!!!!!!!!!!!!!!!!!!!! |
| 736 | |
| 737 | ! Update 3D grid |
| 738 | ProblemRegion3D%U(i,j,:,:) = ProblemRegion%U(:,:) |
| 739 | ! switch speeds back |
| 740 | ProblemRegion3D%U(i,j,:,imom(ix)) = ProblemRegion%U(:,imom(iz)) |
| 741 | ProblemRegion3D%U(i,j,:,imom(iz)) = ProblemRegion%U(:,imom(ix)) |
| 742 | |
| 743 | end do |
| 744 | end do |
| 745 | }}} |
| 746 | |
| 747 | and time is updated and output is done if requested. |