Changes between Version 2 and Version 3 of u/adebrech/code

Timestamp:

03/15/17 18:21:33 (8 years ago)

Author:

adebrech

Comment:

—

u/adebrech/code

@@ -40 +40 @@
 }}}
 
-The state vectors W,,L,, and W,,R,, are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program.
-
+The state vectors W,,L,, and W,,R,, are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program. (Note code between comment lines is excerpted from its actual location.)
+
+{{{#!fortran
+L = region(rhoL,uL,pL)
+R = region(rhoR,uR,pR)
+
+!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
+function ConstructRegion(rho,u,p)
+        
+ConstructRegion%W(irho) = rho
+ConstructRegion%W(iu) = u
+ConstructRegion%W(ip) = p
+ConstructRegion%A = 2/((gamma + 1)*ConstructRegion%W(irho))
+ConstructRegion%B = (gamma - 1)*ConstructRegion%W(ip)/(gamma + 1)
+ConstructRegion%c = sqrt(gamma*ConstructRegion%W(ip)/ConstructRegion%W(irho))
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        
+call validate(L,R)
+
+!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
+validate(L,R)
+
+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
+        
+if (critdiff <= R%W(iu) - L%W(iu)) then
+ print *, 'Vacuum results from these initial states. Exiting.'
+ stop
+end if
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+}}}
 
 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.
 
-
-{{{#!latex
-$f_reg(p_*) = (p_* - p_reg) [\frac{2}{(\gamma + 1) \rho_reg (p_* + \frac{\gamma - 1}{\gamma + 1} p_reg)}]$
-}}}
-
+{{{#!fortran
+pstar = pguess(L,R)
+        
+! do while loop (exits once change is less than tolerance)
+do
+ pstar_old = pstar
+ 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
+ delta = 2.*abs((pstar - pstar_old)/(pstar + pstar_old))                                                                ! Relative change in pressure from previous iteration
+ if (delta < tol) exit                                                                                                  ! Exit upon sufficient accuracy
+end do
+}}}
+
+The functions f and df ($f'$) for either the left or right regions are given by
+
+{{{#!latex
+$f_{reg}(p_*) = \begin{cases} (p_* - p_{reg}) [\frac{A_{reg}}{p_* + B_{reg})}], & p_* > p_{reg} \\ f_{reg}(p_*) = \frac{2 c_{reg}}{\gamma - 1}[(\frac{p_*}{p_{reg}})^{\frac{\gamma - 1}{2 \gamma}} - 1], & p_* \leq p_{reg} \end{cases}$
+
+$f'_{reg}(p_*) = \begin{cases} (\frac{A_{reg}}{B_{reg} + p_*})^{\frac{1}{2}}[1 - \frac{p_* - p_{reg}}{2(B_{reg} + p_*)}], & p_* > p_{reg} \\ \frac{1}{\rho_{reg} c_{reg}}(\frac{p_*}{p_{reg}})^{-\frac{\gamma + 1}{2 \gamma}}, & p_* \leq p_{reg} \end{cases}$
+}}}
 
 See below for derivation of function f.
 
+The velocity and densities for the star regions can then be calculated:
+
+{{{#!fortran
+! Calculates velocity in star region
+exactVelocity(L,R,pstar)
+        
+exactVelocity = 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
+
+!----------------------------------------
+
+! Calculates densities in star regions
+exactDensity(starReg,reg)
+        
+! Density is different for shocks and rarefactions - may differ on either side of CD
+if (starReg%W(ip) > reg%W(ip)) then
+ exactDensity = reg%W(irho)*((gamma - 1.)/(gamma + 1.) + (starReg%W(ip)/reg%W(ip)))/(((gamma - 1.)/(gamma + 1.))*(starReg%W(ip)/reg%W(ip)) + 1.)
+else
+ exactDensity = reg%W(irho)*(starReg%W(ip)/reg%W(ip))**(1./gamma)
+end if
+}}}
+
+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:
+
+{{{#!fortran
+SL = L%W(iu) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma))      ! Left shock speed
+SHL = L%W(iu) - L%c                                                                             ! Head speed of left rarefaction
+STL = Lstar%W(iu) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! Tail speed of left rarefaction
+SR = R%W(iu) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma))      ! right
+SHR = R%W(iu) + R%c                                                                             ! right
+STR = Rstar%W(iu) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! right
+}}}
+
+Then a large conditional is entered to determine which region we are sampling at the given time and position:
+
+{{{#!fortran
+! End once t > finaltime (last output will be t = finaltime)
+do while (t <= finaltime)
+ !print fileformat, "out/result",frame,".dat"
+ write(filename,fileformat) "out/result",frame,".dat"   ! filename to output - goes to out/resultxxxxx.dat
+ open(unit=40,file=filename,iostat=stat,status='replace')
+ write (40,*) 'x rho u p'                               ! write column headers
+ x = xL
+ ! Calculate values for each position - similar to time above
+ !!! could parallelize if I changed from while to for loop
+ do while (x <= xR)
+  if (t == 0) then      !initial state (so we don't divide by zero)
+   if (x < 0) write (40,*) x, L%W
+   if (x > 0) write (40,*) x, R%W
+   if (x == 0) write (40,*) x, 0, 0, 0
+  else
+   if (x/t <= Lstar%W(iu)) then !left side of contact
+   lright = .false.
+    if (Lstar%W(ip) > L%W(ip)) then !shock on left side
+     if (x/t <= SL) then !left side of shock
+      write (40,*) x, L%W
+     else
+      write (40,*) x, Lstar%W
+     end if
+    else !rarefaction on left side
+     if (x/t <= SHL) then !left of the rarefaction fan
+      write (40,*) x, L%W
+     else if (x/t >= STL) then !right of the rarefaction fan
+      write (40,*) x, Lstar%W
+     else !in rarefaction fan
+      call rarefaction(L,x,t,Wfan,lright)
+
+!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
+rarefaction(reg,x,t,Wfan,lright)
+
+! Rarefaction fan differs depending on whether it's a right or left rarefaction
+if (.not. lright) then
+ Wfan(irho) = reg%W(irho)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
+ Wfan(iu) = (2/(gamma + 1))*(reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
+ Wfan(ip) = reg%W(ip)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
+else
+ Wfan(irho) = reg%W(irho)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
+ Wfan(iu) = (2/(gamma + 1))*(-reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
+ Wfan(ip) = reg%W(ip)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
+end if
+!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+      write (40,*) x, Wfan
+     end if
+    end if
+   ! Repeat entire sequence above for right side (lright = .true.)
+   end if
+  end if
+  x = x + dx            ! increment x to next position
+ end do
+ t = t + dt             ! increment t to next time
+ frame = frame + 1      ! update frame counter
+end do
+}}}
+
+These output files can then be processed by VisIt, Mathematica, etc. to make nice pictures.
+
+== Deriving the function f and the Newton-Raphson procedure ==
+{{{#!comment
+Really these should be two separate sections, but I already had a nice transition between them.
+}}}
+
+=== f ===
+
 From the shock jump conditions 
 
@@ -97 +242 @@
 The same analysis can be performed on the right side of the system to yield the same results.
 
+=== Newton-Raphson ===
+
 The Riemann problem can then be solved by finding the value $p = p_*$ such that