Changes between Version 4 and Version 5 of u/adebrech/code


Ignore:
Timestamp:
03/17/17 15:47:17 (8 years ago)
Author:
adebrech
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • u/adebrech/code

    v4 v5  
    1 = [attachment:ExactSolver.tar Exact Riemann Solver] (no vacuum) =
    2 
    3 The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional ideal gases separated at the origin x=0, with initial densities, velocities, and pressures given as input. It gives a solution to the set of equations
     1= [attachment:ExactSolver.tar Exact Riemann Solver] =
     2
     3The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional gases (ideal or with a covolume parameter b) separated at the origin x=0, with initial densities, velocities, and pressures given as input (vacuum is allowed, both as input and as a result of system evolution). It gives a solution to the set of equations
    44
    55{{{#!latex
    … …  
    2121}}}
    2222
    23 The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state.
     23The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state. If vacuum is present, one or more rarefaction fans connect the unperturbed states with the region of vacuum.
    2424
    2525The program reads input from a file called problem.data, and stores many of the relevant parameters in global variables.
    … …  
    4848function ConstructRegion(rho,u,p)
    4949       
     50! Fill state
    5051ConstructRegion%W(irho) = rho
    5152ConstructRegion%W(iu) = u
    5253ConstructRegion%W(ip) = p
    53 ConstructRegion%A = 2/((gamma + 1)*ConstructRegion%W(irho))
    54 ConstructRegion%B = (gamma - 1)*ConstructRegion%W(ip)/(gamma + 1)
    55 ConstructRegion%c = sqrt(gamma*ConstructRegion%W(ip)/ConstructRegion%W(irho))
     54! Calculate useful constants
     55if (.not. (rho == 0 .and. p == 0)) then              ! state is not vacuum
     56 ConstructRegion%A = 2.*(1. - cov*ConstructRegion%W(irho))/((gamma + 1.)*ConstructRegion%W(irho))
     57 ConstructRegion%B = (gamma - 1.)*ConstructRegion%W(ip)/(gamma + 1.)
     58 ConstructRegion%c = sqrt(gamma*ConstructRegion%W(ip)/(ConstructRegion%W(irho)*(1. - cov*ConstructRegion%W(irho))))
     59else
     60 ConstructRegion%c = 0          ! Set speed of sound to zero if state is vacuum (to avoid NaN)          !!! could also set A,B to zero, but since they don't represent physical quantities, not much point
     61end if
    5662!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     63}}}
     64
     65The initial states are then checked for vacuum, and the state of the system (initial vacuum, vacuum created, or no vacuum) is set.
     66
     67{{{#!fortran
     68call check_vacuum()
     69
     70!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     71check_vacuum()
     72
     73critdiff = 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
    5774       
    58 call validate(L,R)
    59 
    60 !∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
    61 validate(L,R)
    62 
    63 critdiff = 2.*L%c/(gamma - 1.) + 2.*R%c/(gamma - 1.)      ! Difference between right and left velocities must be strictly less than this to prevent vacuum formation
    64        
    65 if (critdiff <= R%W(iu) - L%W(iu)) then
    66  print *, 'Vacuum results from these initial states. Exiting.'
    67  stop
     75if ((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
     76 vacuum = 1
     77 print *, 'One of initial states is vacuum.'
     78else if (critdiff <= R%W(iu) - L%W(iu)) then             ! vacuum is created by initial states
     79 vacuum = 2
     80 print *, 'Vacuum is created by initial states.'
     81else            ! no vacuum
     82 vacuum = 3
    6883end if
    6984!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    7085}}}
    7186
    72 These state vectors are then passed to the Newton-Raphson solver, which calls a subroutine to guess a good value for p,,0,, (currently not very effective - just uses tol) and then iterates until the tolerance for change in pressure between steps is met.
    73 
    74 {{{#!fortran
    75 pstar = pguess(L,R)
    76        
     87If there is no vacuum, the state vectors L and R are then passed to the Newton-Raphson solver, along with a guess at a good value for p,,0,, (currently not very effective - just uses the average value of the initial states). The solver takes a guess at a root for a function, the function itself, and its derivative, and returns the root of the function (to within specified tolerance). If a negative value is encountered at any point during iteration, the guess is reseeded with the tolerance of the solver (since the two values solved for are density and pressure, negative solutions would be unphysical - would be pretty easy to rework to avoid NaNs if negative values were allowed).
     88
     89{{{#!fortran
     90! Use Newton-Raphson solver to find pressure in star region
     91Lstar%W(ip) = NRSolver(pguess(L,R),covolumeF,covolumeDf)
     92
     93!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     94! Newton-Raphson iterative solver. Takes guess of root of equation, equation, and derivative and calculates root to tol.
     95function NRSolver(guess,f,df)
     96
    7797! do while loop (exits once change is less than tolerance)
    7898do
    79  pstar_old = pstar
    80  pstar = pstar_old - ((f(pstar_old,L) + f(pstar_old,R) + R%W(iu)-L%W(iu))/(df(pstar_old,L) + df(pstar_old,R)))          ! Newton-Raphson formula
    81  delta = 2.*abs((pstar - pstar_old)/(pstar + pstar_old))                                                                ! Relative change in pressure from previous iteration
    82  if (delta < tol) exit                                                                                                  ! Exit upon sufficient accuracy
     99 guess_old = guess
     100 guess = guess_old - f(guess_old)/df(guess_old)           ! Newton-Raphson formula
     101 delta = 2.*abs((guess - guess_old)/(guess + guess_old))  ! Relative change in pressure from previous iteration
     102 if (guess < 0) guess = tol                               ! If we get a negative value (unphysical for our applications), set guess to tol      !!! could use nan inequality for more general case
     103 if (delta < tol) then                                    ! Exit upon sufficient accuracy
     104  NRSolver = guess
     105  exit
     106 end if
    83107end do
     108!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    84109}}}
    85110
    … …  
    92117}}}
    93118
    94 See below for derivation of function f.
    95 
    96 The velocity and densities for the star regions can then be calculated:
    97 
    98 {{{#!fortran
     119See below for derivation of function f (in the ideal gas case. The derivation is identical for the covolume case, except as noted below).
     120
     121Once the pressure is solved for, velocity and densities for the star regions can then be calculated:
     122
     123{{{#!fortran
     124! Calculate velocity and densities in star regions
     125Lstar%W(iu) = exactVelocity(Lstar%W(ip))
     126
     127!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
    99128! Calculates velocity in star region
    100129exactVelocity(L,R,pstar)
    101130       
    102131exactVelocity = 0.5*(L%W(iu) + R%W(iu)) + 0.5*(f(pstar,R) - f(pstar,L)) ! Velocity in star region is the same on either side of CD
    103 
    104 !----------------------------------------
    105 
     132!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     133
     134Rstar%W = Lstar%W
     135Lstar%W(irho) = exactDensity(Lstar,L)
     136Rstar%W(irho) = exactDensity(Rstar,R)
     137
     138!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
    106139! Calculates densities in star regions
    107140exactDensity(starReg,reg)
    … …  
    113146 exactDensity = reg%W(irho)*(starReg%W(ip)/reg%W(ip))**(1./gamma)
    114147end if
    115 }}}
    116 
    117 The final value of pstar is printed to console, then output is written to files, one for each time state requested. During output, the speeds of shocks and rarefaction boundaries are calculated:
     148!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     149}}}
     150
     151The final value of pstar is then printed to console and output is written to files. If vacuum is present, no calculation is required; the output can just be written to files. One file is written for each time state requested. During output for the case with no vacuum, the speeds of shocks and rarefaction boundaries are calculated:
    118152
    119153{{{#!fortran
    … …  
    126160}}}
    127161
    128 Then a large conditional is entered to determine which region we are sampling at the given time and position:
     162Then a large conditional is entered to determine which region we are sampling at the given time and position. If we are located in a rarefaction fan, the density, velocity, and pressure are calculated, which requires the solution of an equation for the density using the Newton-Raphson solver if the fluid has a covolume parameter:
    129163
    130164{{{#!fortran
    … …  
    161195
    162196!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
    163 rarefaction(reg,x,t,Wfan,lright)
    164 
    165 ! Rarefaction fan differs depending on whether it's a right or left rarefaction
    166 if (.not. lright) then
    167  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
    168  Wfan(iu) = (2/(gamma + 1))*(reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
    169  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
     197subroutine rarefaction(reg,x,t,Wfan,lright)
     198
     199! If fluid has covolume constant, iterate to solve for density
     200if (cov > 0) then
     201 if (.not. lright) then
     202  beta = ((((gamma - 1.)*(L%W(iu) + 2.*L%c*(1. - cov*L%W(irho))/(gamma - 1.) - x/t))**2.)/(gamma*L%W(ip)*((1. - cov*L%W(irho))/L%W(irho))**gamma)) ! Constant related to Riemann invariant
     203  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf)
     204  Wfan(ip) = reg%W(ip)*(((1. - cov*reg%W(irho))/reg%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
     205  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho))))                                         ! Speed of sound at current location
     206  Wfan(iu) = x/t + c
     207 else
     208  beta = (((x/t - R%W(iu) + (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
     209  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf)
     210  Wfan(ip) = reg%W(ip)*(((1. - cov*reg%W(irho))/reg%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
     211  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho))))                                         ! Speed of sound at current location
     212  Wfan(iu) = x/t - c
     213 end if
     214! Otherwise solution is closed-form
    170215else
    171  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
    172  Wfan(iu) = (2/(gamma + 1))*(-reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
    173  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
     216 ! Rarefaction fan differs depending on whether it's a right or left rarefaction
     217 if (.not. lright) then
     218  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
     219  Wfan(iu) = (2/(gamma + 1))*(reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
     220  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
     221 else
     222  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
     223  Wfan(iu) = (2/(gamma + 1))*(-reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
     224  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
     225 end if
    174226end if
    175227!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    … …  
    190242These output files can then be processed by VisIt, Mathematica, etc. to make nice pictures.
    191243
    192 View [http://www.pas.rochester.edu/~adebrech/code/ online], or download [attachment:test1.mp4 1], [attachment:test2.mp4 2], [attachment:test3.mp4 3], [attachment:test4.mp4 4], [attachment:test5.mp4 5]
     244View [http://www.pas.rochester.edu/~adebrech/code/ online] or download [attachment:test1.mp4 1], [attachment:test2.mp4 2], [attachment:test3.mp4 3], [attachment:test4.mp4 4], [attachment:test5.mp4 5] the tests from Toro. The vacuum routines have only been sanity tested, while the covolume routines have been tested once against the case in Table III(b) in Toro, E.F. "A Fast Riemann Solver with Constant Covolume Applied to the Random Choice Method," Int. J. Num. Meth. in Fluids, Vol. 9 (1989).
    193245
    194246== Deriving the function f and the Newton-Raphson procedure ==
    … …  
    243295The same analysis can be performed on the right side of the system to yield the same results.
    244296
     297==== Covolume ====
     298
     299In the covolume case, the internal energy becomes
     300
     301{{{#!latex
     302$e = \frac{p(1 - b \rho)}{\rho(\gamma - 1)}$
     303}}}
     304
     305and the sound speed is
     306
     307{{{#!latex
     308$c = (\frac{p \gamma}{\rho (1 - b \rho)})^\frac{1}{2}$
     309}}}
     310
     311which reduce to the ideal gas case when $b = 0$.
     312
    245313=== Newton-Raphson ===
    246314