Changes between Version 10 and Version 11 of u/adebrech/code


Ignore:
Timestamp:
05/09/17 10:31:17 (8 years ago)
Author:
adebrech
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • u/adebrech/code

    v10 v11  
    11= [attachment:AllCode.tar Godunov Solver for Euler Equations] =
    22
    3 
     3This 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
     6NAMELIST/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.
     9inquire(file='problem.data', exist=ex)
     10if (ex) then
     11 OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
     12 READ(30,NML=ProblemData)                               ! read into variables in namelist
     13 CLOSE(30)
     14else
     15 print *, 'Missing problem.data. Exiting.'
     16 stop           ! exit if missing
     17end if
     18}}}
     19
     20The 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!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     96PrimToCons(irho) = W(irho)
     97PrimToCons(imom(ix)) = W(iu(ix))*W(irho)
     98PrimToCons(iE) = W(ip)*(1. - cov*W(irho))/(gamma - 1.) + 0.5*W(irho)*sum(W(iu)**2)
     99! Tangential directins
     100if (nDim > 1) PrimToCons(imom(iy)) = W(iu(iy))*W(irho)
     101if (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
     179After 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!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     192ConsToPrim(irho) = U(irho)
     193ConsToPrim(iu(ix)) = U(imom(ix))/U(irho)
     194ConsToPrim(ip) = (U(iE) - (sum(U(imom)**2))/(2*U(irho)))*(gamma - 1.)/(1. - cov*U(irho))
     195! Tangential directions
     196if (nDim > 1) ConsToPrim(iu(iy)) = U(imom(iy))/U(irho)
     197if (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
     227The 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
     274The 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
     277select 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)
     293end select
     294}}}
     295
     296and 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!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     317p0 = max(0d0, PrimitiveVariablePressure(L,R))
     318! Approximate shock or rarefaction speeds
     319if (p0 > L%W(ip)) then
     320 SL = L%W(iu(ix)) - L%c*sqrt(1 + (gamma + 1.)*(p0/L%W(ip) - 1.)/(2.*gamma))
     321else
     322 SL = L%W(iu(ix)) - L%c
     323end if
     324if (p0 > R%W(ip)) then
     325 SR = R%W(iu(ix)) + R%c*sqrt(1 + (gamma + 1.)*(p0/R%W(ip) - 1.)/(2.*gamma))
     326else
     327 SR = R%W(iu(ix)) + R%c
     328end if
     329if (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
     336if (0 <= SL) then
     337 HLLCflux = EulerFlux(L%W)
     338! Approximate with integral average of left star region
     339else 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
     348else 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
     357else
     358 HLLCflux = EulerFlux(R%W)
     359end 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
     368if (0 <= SL) then
     369 HLLflux = EulerFlux(L%W)
     370! Approximate with integral average of star region
     371else 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
     376else
     377 HLLflux = EulerFlux(R%W)
     378end 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!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     386function EulerRiemannSoln
     387
     388! Check for vacuum
     389call check_vacuum(L,R)
     390
     391!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     392critdiff = 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
     394if ((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.'
     397else 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.'
     400else    ! no vacuum
     401 vacuum = 3
     402end if
     403!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     404
     405! Vacuum present, don't need star region - exact solution is cheap
     406if (vacuum /= 3) then
     407 EulerRiemannSoln = sample(L,R,0d0)
     408
     409!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     410! Calculate shock and rarefaction speeds
     411SvacL = L%W(iu(ix)) + 2.*L%c/(gamma - 1.)   ! Speed of vacuum front (vacuum on right)
     412SvacR = R%W(iu(ix)) - 2.*R%c/(gamma - 1.)   ! Speed of vacuum front (vacuum on left)
     413
     414! Initial state is vacuum
     415if (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
     428if (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
     447else
     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
     464end 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
     482else 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
     502end if
     503!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     504
     505! No vacuum, calculate star region and sample along t axis
     506else
     507 ! Star region may be approximated
     508 call StarRegion(L,R,Lstar,Rstar)
     509
     510!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     511! Passively advected - no math, so no approximations
     512if (nDim > 1) then
     513 Lstar%W(iu(iy)) = L%W(iu(iy))
     514 Rstar%W(iu(iy)) = R%W(iu(iy))
     515end if
     516if (nDim > 2) then
     517 Lstar%W(iu(iz)) = L%W(iu(iz))
     518 Rstar%W(iu(iz)) = R%W(iu(iz))
     519end if
     520
     521select 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
     571end select
     572!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     573
     574 EulerRiemannSoln = sample(L,R,Lstar,Rstar,0d0)
     575
     576!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     577! Calculate shock and rarefaction speeds
     578SL = L%W(iu(ix)) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma))      ! Left shock speed
     579SHL = L%W(iu(ix)) - L%c     ! Head speed of left rarefaction
     580STL = Lstar%W(iu(ix)) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))  ! Tail speed of left rarefaction
     581SR = R%W(iu(ix)) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma))      ! right
     582SHR = R%W(iu(ix)) + R%c     ! right
     583STR = Rstar%W(iu(ix)) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))  ! right
     584
     585if (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
     603else !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
     621end if
     622!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     623
     624 deallocate(Lstar%W,Rstar%W)
     625end if
     626!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     627!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     628EulerFlux(irho) = W(irho)*W(iu(ix))
     629EulerFlux(imom(ix)) = W(irho)*W(iu(ix))**2 + W(ip)
     630EulerFlux(iE) = W(iu(ix))*((W(ip)*(1. - cov*W(ip)))/(gamma - 1.) + 0.5*W(irho)*W(iu(ix))**2 + W(ip))
     631if (nDim > 1) EulerFlux(imom(iy)) = W(irho)*W(iu(ix))*W(iu(iy))
     632if (nDim > 2) EulerFlux(imom(iz)) = W(irho)*W(iu(ix))*W(iu(iz))
     633!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     634
     635   end select
     636  end do
     637}}}
     638
     639The grid is then updated using Godunov's method to find the new state of the system:
     640
     641{{{#!fortran
     642do 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
     648end do
     649}}}
     650
     651The 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
     658The 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
     666Once 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
     669do 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
     676The 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
     687Once 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
     747and time is updated and output is done if requested.
    4748
    5749= [attachment:GodunovSolver.tar Godunov Solver] =