Changes between Version 7 and Version 8 of u/adebrech/code


Ignore:
Timestamp:
04/14/17 11:35:23 (8 years ago)
Author:
adebrech
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • u/adebrech/code

    v7 v8  
    11= [attachment:GodunovSolver.tar Godunov Solver] =
    22
    3 
     3This 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
     6NAMELIST/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
     20The domain is then initialized according to the selected conditions:
     21
     22{{{#!fortran
     23select 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
     50end select
     51}}}
     52
     53The boundary conditions are then applied:
     54
     55{{{#!fortran
     56select 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)
     63end select
     64}}}
     65
     66After 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
     70maxSpeed = 0.
     71do i = 0, cells+1
     72 if (abs(ProblemRegion%u(i)) > maxSpeed) maxSpeed = abs(ProblemRegion%u(i))     ! (See Toro, Ch. 5.6)
     73end do
     74
     75! Calculate delta t
     76dt = CFL*dx/maxSpeed
     77
     78! If time step would pass the next output frame, recalculate
     79if (time + dt > nexttime) dt = nexttime - time
     80}}}
     81
     82The intercell flux is then calculated and the velocity array is updated to the next time.
     83
     84{{{#!fortran
     85do 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!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
     93function RiemannSoln(uL,uR)
     94
     95if (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
     102else
     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
     110end if
     111!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     112
     113  case (inviscidBurgers)
     114   ProblemRegion%flux(i) = BurgersFlux(RiemannSoln(uL,uR))
     115 end select
     116end do
     117
     118! Loop through cells using previously calculated flux
     119do i = 1, cells
     120 ProblemRegion%u(i) = ProblemRegion%u(i) + (dt/dx)*(ProblemRegion%flux(i-1) - ProblemRegion%flux(i))
     121end do
     122}}}
     123
     124Finally, at each requested output frame, the current system state is written to a file:
     125
     126{{{#!fortran
     127write (40,*) 'x u'       ! write column headers
     128do 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)
     131end do
     132}}}
    4133
    5134= [attachment:ExactSolver.tar Exact Riemann Solver] =