UC Berkeley Math 228B, Spring 2024, Problem Set 2

1. Consider the PDE \[\begin{aligned} u_t = \kappa u_{xx} - \gamma u, &&(1) \end{aligned}\] which models a diffusion with decay provided \(\kappa>0\) and \(\gamma>0\). Consider methods of the form \[\begin{aligned} U_j^{n+1} =U_j^n+{\frac{k\kappa}{2h^2}} [U_{j-1}^n-2U_j^n +U_{j+1}^n +U_{j-1}^{n+1} - 2U_j^{n+1} + U_{j+1}^{n+1}] - k\gamma[(1-\theta)U_j^n + \theta U_j^{n+1}] &&(2) \end{aligned}\] where \(\theta\) is a parameter. In particular, if \(\theta=1/2\) then the decay term is modeled with the same centered-in-time approach as the diffusion term and the method can be obtained by applying the Trapezoidal method to the MOL formulation of the PDE. If \(\theta=0\) then the decay term is handled explicitly. For more general reaction-diffusion equations it may be advantageous to handle the reaction terms explicitly since these terms are generally nonlinear, so making them implicit would require solving nonlinear systems in each time step (whereas handling the diffusion term implicitly only gives a linear system to solve in each time step).

   1. (In Problem Set 1)

   2. Using von Neumann analysis, show that this method is unconditionally stable if \(\theta \geq 1/2\).

   3. Show that if \(\theta = 0\) then the method is stable provided \(k\leq 2/\gamma\), independent of \(h\).

2. Suppose \(a>0\) and consider the following skewed leapfrog method for solving the advection equation \(u_t + au_x = 0\): \[\begin{aligned} U_j^{n+1} = U_{j-2}^{n-1} - \left(\frac{ak}{h} - 1\right) (U_j^n - U_{j-2}^n). &&(3) \end{aligned}\] Note that if \(ak/h \approx 1\) then this stencil roughly follows the characteristic of the advection equation and might be expected to be more accurate than standard leapfrog. (If \(ak/h = 1\) the method is exact.)

   

   1. What is the order of accuracy of this method?

   2. For what range of Courant number \(ak/h\) does this method satisfy the CFL condition?

   3. Show that the method is in fact stable for this range of Courant numbers by doing von Neumann analysis. Hint: Let \(\gamma(\xi) = e^{i\xi h}g(\xi)\) and show that \(\gamma\) satisfies a quadratic equation closely related to the equation (10.34) that arises from a von Neumann analysis of the leapfrog method.

3. Consider Euler’s equations of compressible gas dynamics in two space dimensions: \[\begin{aligned} u_t + \nabla \cdot F = 0,\quad\text{where } u=\begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{bmatrix} \text{ and } F=\begin{bmatrix} \rho u & \rho v \\ \rho u^2+p & \rho u v \\ \rho u v & \rho v^2 + p \\ u (\rho E+p) & v(\rho E+p) \end{bmatrix} && (4) \end{aligned}\] Here, \(\rho\) is the fluid density, \(u, v\) are the velocity components, and \(E\) is the total energy. For an ideal gas, the pressure \(p\) has the form \(p = (\gamma-1) \rho ( E - (u^2 + v^2)/2 )\), where \(\gamma\) is the adiabatic gas constant. We will solve these on a square domain with periodic boundary conditions, for \(0\le t \le T\). The spatial derivatives will be discretized with a fourth order compact Padé scheme, and the solution will be filtered using a sixth order compact filter. A standard RK4 scheme will be used for time integration.

   1. Write a function with the syntax

         Frx, Fry, Frux, Fruy, Frvx, Frvy, FrEx, FrEy =  euler_fluxes(r, ru, rv, rE)

      which returns the 8 flux functions in (4) for the 4 solution components. Note that all the inputs and outputs are arrays of the same size, and the fluxes are evaluated elementwise. Assume \(\gamma = 7/5\).

   2. Write a function with the syntax
       
      divF = compact_div(Fx, Fy, h)
       
      which calculates the divergence of a grid function field \(F=[F_x,F_y]\) using the 4th order compact Padé scheme with periodic boundary conditions and grid spacing \(h\): \[\begin{aligned} f'_{i-1} + 4f'_i + f'_{i+1} = 3 \frac{f_{i+1}-f_{i-1}}{h} \end{aligned}\]

   3. Write a function with the syntax
       
      u = compact_filter(u,\(\alpha\))
       
      which filters the grid solution \(u\) using the 6th order compact filter with parameter \(\alpha\): \[\begin{aligned} \alpha \hat{f}_{i-1} + \hat{f}_i + \alpha \hat{f}_{i+1} = af_i + \frac{c}{2} (f_{i+2}+f_{i-2}) + \frac{b}{2}(f_{i+1}+f_{i-1}) \end{aligned}\] where \(a=5/8 + 3\alpha/4\), \(b=\alpha + 1/2\), \(c=\alpha/4 - 1/8\)

   4. Write a function with the syntax

         fr, fru, frv, frE = euler_rhs(r, ru, rv, rE, h)

      which computes the right-hand side of the discretized \(-\nabla\cdot F\) (essentially just calling euler_fluxes and compact_div).

   5. Write a function with the syntax
       
      r, ru, rv, rE = euler_rk4step(r, ru, rv, rE, h, k,\(\alpha\))
       
      which takes one RK4 step using euler_rhs and filters each solution component using compact_filter.

   6. Validate your solver with the function euler_vortex (in a notebook on the course web page). Use a square domain \(0\le x,y \le 10\) with grid spacings \(h=10/N\) and \(N=32, 64, 128\). Use the time step \(k\le 0.3h\), adjusted so the final time \(T=5\sqrt{2}\) is a multiple of \(k\). Use the initial solution:

         pars = [0.5, 1, 0.5, pi/4, 2.5, 2.5]
         r, ru, rv, rE = euler_vortex(x, y, 0, pars)

      and compare with the exact final solution:

         r0, ru0, rv0, rE0 = euler_vortex(x, y, 5sqrt(2), pars)

      Calculate the errors in the infinity norm over all solution components. Plot the errors vs. \(h\) in a log-log plot, for the 3 grid spacings \(h\) and the two filter coefficients \(\alpha = 0.499\), \(\alpha=0.48\). Estimate the slopes of the two curves.

   7. Simulate a Kelvin-Helmholtz instability, using the unit square domain \(0\le x,y\le 1\) with grid spacing \(h=1/N\) and \(N=256\) grid points in each coordinate direction, \(\alpha=0.48\), time step \(k\le 0.3h\), final time \(T=1\), and the initial condition: \[\begin{aligned} \rho &= \begin{cases} 2 & \text{if } |y-0.5| < 0.15 + \sin(2\pi x) / 200, \\ 1 & \text{otherwise}. \end{cases} \\ u &= \rho - 1,\quad v = 0,\quad p = 3 \end{aligned}\] Plot the final solution using a contour or color plot of the density \(\rho\).

Code Submission: Your Julia file needs to define the functions euler_fluxes, compact_div, compact_filter, euler_rhs, and euler_rh4step, and any other supporting functions and variables.