Source code for pyro.burgers_viscous.interface

from pyro.multigrid import MG


[docs] def get_lap(grid, a): r""" Parameters ---------- grid: grid object a: ndarray the variable that we want to find the laplacian of Returns ------- out : ndarray (laplacian of state a) """ lap = grid.scratch_array() lap.v(buf=2)[:, :] = (a.ip(1, buf=2) - 2.0*a.v(buf=2) + a.ip(-1, buf=2)) / grid.dx**2 \ + (a.jp(1, buf=2) - 2.0*a.v(buf=2) + a.jp(-1, buf=2)) / grid.dy**2 return lap
[docs] def diffuse(my_data, rp, dt, scalar_name, A): r""" A routine to solve the Helmhotlz equation with constant coefficient and update the state. (a + b \lap) phi = f using Crank-Nicolson discretization with multigrid V-cycle. Parameters ---------- my_data : CellCenterData2d object The data object containing the grid and advective scalar that we are advecting. rp : RuntimeParameters object The runtime parameters for the simulation dt : float The timestep we are advancing through. scalar_name : str The name of the variable contained in my_data that we are advecting A: ndarray The advective source term for diffusion Returns ------- out : ndarray (solution of the Helmholtz equation) """ myg = my_data.grid a = my_data.get_var(scalar_name) eps = rp.get_param("diffusion.eps") # Create the multigrid with a constant diffusion coefficient # the equation has form: # (alpha - beta L) phi = f # # with alpha = 1 # beta = (dt/2) k # f = phi + (dt/2) k L phi - A # # Same as Crank-Nicolson discretization except with an extra # advection source term) mg = MG.CellCenterMG2d(myg.nx, myg.ny, xmin=myg.xmin, xmax=myg.xmax, ymin=myg.ymin, ymax=myg.ymax, xl_BC_type=my_data.BCs[scalar_name].xlb, xr_BC_type=my_data.BCs[scalar_name].xrb, yl_BC_type=my_data.BCs[scalar_name].ylb, yr_BC_type=my_data.BCs[scalar_name].yrb, alpha=1.0, beta=0.5*dt*eps, verbose=0) # Compute the RHS: f f = mg.soln_grid.scratch_array() lap = get_lap(myg, a) f.v()[:, :] = a.v() + 0.5*dt*eps * lap.v() - dt*A.v() mg.init_RHS(f) # initial guess is zeros mg.init_zeros() mg.solve(rtol=1.e-12) # perform the diffusion update a.v()[:, :] = mg.get_solution().v()
[docs] def apply_diffusion_corrections(grid, dt, eps, u, v, u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr): r""" Apply diffusion correction term to the interface state .. math:: u_t + u u_x + v u_y = eps (u_xx + u_yy) v_t + u v_x + v v_y = eps (v_xx + v_yy) We use a second-order (piecewise linear) unsplit Godunov method (following Colella 1990). Our convection is that the fluxes are going to be defined on the left edge of the computational zones:: | | | | | | | | -+------+------+------+------+------+------+-- | i-1 | i | i+1 | a_l,i a_r,i a_l,i+1 a_r,i and a_l,i+1 are computed using the information in zone i,j. Parameters ---------- grid : Grid2d The grid object dt : float The timestep eps : float The viscosity u_xl, u_xr : ndarray ndarray left and right states of x-velocity in x-interface. u_yl, u_yr : ndarray ndarray left and right states of x-velocity in y-interface. v_xl, v_xr : ndarray ndarray left and right states of y-velocity in x-interface. v_yl, u_yr : ndarray ndarray left and right states of y-velocity in y-interface. Returns ------- out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray unsplit predictions of the left and right states of u and v on both the x- and y-interfaces along with diffusion correction terms. """ #apply diffusion correction to the interface lap_u = get_lap(grid, u) lap_v = get_lap(grid, v) u_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2) u_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2) u_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2) u_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2) v_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2) v_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2) v_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2) v_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2) return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr