Source code for pyrodeo.roe

# -*- coding: utf-8 -*-
"""Defines the Roe class defining the Roe solver.
"""

from __future__ import print_function

import numpy as np
from pyrodeo.claw_solver import ClawSolver

[docs]class Roe(ClawSolver): """Construct class for the Roe solver. Constructor sets two basic attributes (sb and min_dens), after which a time step can be taken through the :meth:`.step` method. Args: flux_limiter (float): Flux limiter parameter. Should be >= 1 (minmod, most diffusive limiter) and <= 2 (superbee, least diffusive limiter). min_dens (float): Minimum density when to switch to HLL to preserve positivity. The following attributes and methods are available: Attributes: sb (float): Flux limiter parameter. Should be >= 1 (minmod, most diffusive limiter) and <= 2 (superbee, least diffusive limiter). min_dens (float): Minimum density when to switch to HLL to preserve positivity. """ def __init__(self, flux_limiter, min_dens): self.sb = flux_limiter self.min_dens = min_dens
[docs] def limit_flux(self, dens, dens_left, f1dens, f2dens, dtdx): """Limit second order flux to preserve positivity Calculate the maximum contribution of the second order flux in order for the density to remain positive. Args: dens (ndarray): Current density. dens_left (ndarray): Density in the cell to the left f1dens (ndarray): First order mass flux. f2dens (ndarray): Second order mass flux. dtdx (float): Time step / space step Returns: ndarray: Array with values >= 0 and <= 1 specifying the maximum contribution of second order flux for the density to remain positive. """ thetap = 0.0*dens + 1.0 rhop = 0.5*(dens_left - 2.0*dtdx*f2dens) epsd = np.minimum(0.0*dens + self.min_dens, np.minimum(dens_left, dens)) sel = np.where(rhop < epsd) rho1 = 0.5*(dens_left[sel] - 2.0*dtdx*f1dens[sel]) if len(rho1) >= 1: if np.min(rho1) < 0.0: dens[sel] = -1.0*dens[sel] thetap[sel] = np.maximum(0.0*rho1, (rho1 - epsd[sel])/(rho1 - rhop[sel])) if np.max(thetap) > 1.0: dens = -1.0*dens thetam = 0.0*dens + 1.0 rhom = 0.5*(dens + 2.0*dtdx*f2dens) sel = np.where(rhom < epsd) rho1 = 0.5*(dens[sel] + 2.0*dtdx*f1dens[sel]) if len(rho1) >= 1: if np.min(rho1) < 0.0: dens[sel] = -1.0*dens[sel] thetam[sel] = np.maximum(0.0*rho1, (rho1 - epsd[sel])/(rho1 - rhom[sel])) if np.max(thetam) > 1.0: dens = -1.0*dens return np.minimum(thetap, thetam)
[docs] def step(self, dt, dx, state, source, bc): """Update state for a single time step. Args: dt (float): Time step. dx (float): Space step. state (:class:`.State`): Current :class:`.State`, will be updated. source (ndarray): Geometric source terms, must have same shape as state.dens. bc (str, str): Boundary conditions (in and out): 'periodic', 'symmetric', or 'closed' (other boundary conditions are dealt with elsewhere). """ # Set periodic boundaries if bc[0] == 'periodic': state.dens[:2,:,:] = state.dens[-4:-2,:,:] state.velx[:2,:,:] = state.velx[-4:-2,:,:] state.vely[:2,:,:] = state.vely[-4:-2,:,:] state.velz[:2,:,:] = state.velz[-4:-2,:,:] state.dens[-2:,:,:] = state.dens[2:4,:,:] state.velx[-2:,:,:] = state.velx[2:4,:,:] state.vely[-2:,:,:] = state.vely[2:4,:,:] state.velz[-2:,:,:] = state.velz[2:4,:,:] # Symmetry boundary if bc[0] == 'symmetric': state.dens[:2,:,:] = state.dens[3:1:-1,:,:] state.velx[:2,:,:] = -state.velx[3:1:-1,:,:] state.vely[:2,:,:] = state.vely[3:1:-1,:,:] state.velz[:2,:,:] = state.velz[3:1:-1,:,:] if bc[1] == 'symmetric': state.dens[-2:,:,:] = state.dens[-3:-5:-1,:,:] state.velx[-2:,:,:] = -state.velx[-3:-5:-1,:,:] state.vely[-2:,:,:] = state.vely[-3:-5:-1,:,:] state.velz[-2:,:,:] = state.velz[-3:-5:-1,:,:] dens_left = np.roll(state.dens, 1, axis=0) velx_left = np.roll(state.velx, 1, axis=0) vely_left = np.roll(state.vely, 1, axis=0) velz_left = np.roll(state.velz, 1, axis=0) cs_left = np.roll(state.soundspeed, 1, axis = 0) # Average sound speed 0.5*((i) + (i-1)) cs_left = np.roll(state.soundspeed, 1, axis = 0) cs = 0.5*(state.soundspeed + cs_left) dtdx = dt/dx # Parameter vector z0 = np.sqrt(state.dens) z1 = z0*state.velx z2 = z0*state.vely z3 = z0*state.velz # Average parameter vector: 0.5*((i) + (i-1)) z0 = 0.5*(z0 + np.roll(z0, 1, axis=0)) z1 = 0.5*(z1 + np.roll(z1, 1, axis=0)) z2 = 0.5*(z2 + np.roll(z2, 1, axis=0)) z3 = 0.5*(z3 + np.roll(z3, 1, axis=0)) # Roe averages u = z1/z0 v = z2/z0 w = z3/z0 # HLL signal speeds bm = np.minimum(0.0*u, np.minimum(u - cs, velx_left - cs_left)) bp = np.maximum(0.0*u, np.maximum(u + cs, state.velx + state.soundspeed)) # Cell-centred fluxes fdens = state.dens*state.velx fmomx = fdens*state.velx + \ state.soundspeed*state.soundspeed*state.dens fmomy = fdens*state.vely fmomz = fdens*state.velz # Left and right fluxes (note stationary extrapolation) fldens = np.roll(fdens, 1, axis=0) flmomx = np.roll(fmomx + 0.5*dx*source, 1, axis=0) flmomy = np.roll(fmomy, 1, axis=0) flmomz = np.roll(fmomz, 1, axis=0) frdens = fdens frmomx = fmomx - 0.5*dx*source frmomy = fmomy frmomz = fmomz # Flux difference: (i) - (i-1) dfdens = frdens - fldens dfmomx = frmomx - flmomx dfmomy = frmomy - flmomy dfmomz = frmomz - flmomz # State difference: (i) - (i-1) ddens = state.dens - dens_left dmomx = state.dens*state.velx - dens_left*velx_left dmomy = state.dens*state.vely - dens_left*vely_left dmomz = state.dens*state.velz - dens_left*velz_left # HLL flux f1dens = (bp*fldens - bm*frdens + bp*bm*ddens)/(bp - bm); f1momx = (bp*flmomx - bm*frmomx + bp*bm*dmomx)/(bp - bm); f1momy = (bp*flmomy - bm*frmomy + bp*bm*dmomy)/(bp - bm); f1momz = (bp*flmomz - bm*frmomz + bp*bm*dmomz)/(bp - bm); # Projection coefficients b1 =-0.5*(dfmomx - (u + cs)*dfdens)/cs b2 = 0.5*(dfmomx - (u - cs)*dfdens)/cs b3 = dfmomy - v*dfdens; b4 = dfmomz - w*dfdens; b1 = b1/(u - cs + 1.0e-10) b2 = b2/(u + cs + 1.0e-10) b3 = b3/(u + 1.0e-10) b4 = b4/(u + 1.0e-10) # Upwind projection coefficients b1u = -0.5*(np.sign(u - cs) - 1.0)*np.roll(b1, -1, axis=0) + \ 0.5*(np.sign(u - cs) + 1.0)*np.roll(b1, 1, axis=0) b2u = -0.5*(np.sign(u + cs) - 1.0)*np.roll(b2, -1, axis=0) + \ 0.5*(np.sign(u + cs) + 1.0)*np.roll(b2, 1, axis=0) b3u = -0.5*(np.sign(u) - 1.0)*np.roll(b3, -1, axis=0) + \ 0.5*(np.sign(u) + 1.0)*np.roll(b3, 1, axis=0) b4u = -0.5*(np.sign(u) - 1.0)*np.roll(b4, -1, axis=0) + \ 0.5*(np.sign(u) + 1.0)*np.roll(b4, 1, axis=0) # Flux limiter fl1 = self.limiter(b1, b1u, self.sb) fl2 = self.limiter(b2, b2u, self.sb) fl3 = self.limiter(b3, b3u, self.sb) fl4 = self.limiter(b4, b4u, self.sb) # Second order projection coefficients b1 = -np.abs(u - cs)*b1 + fl1*(np.abs(u - cs) - dtdx*(u - cs)*(u - cs)) b2 = -np.abs(u + cs)*b2 + fl2*(np.abs(u + cs) - dtdx*(u + cs)*(u + cs)) b3 = -np.abs(u)*b3 + fl3*(np.abs(u) - dtdx*u*u) b4 = -np.abs(u)*b4 + fl4*(np.abs(u) - dtdx*u*u) # Interface flux (i - 1/2) f2dens = 0.5*(fldens + frdens + b1 + b2) f2momx = 0.5*(flmomx + frmomx + (u - cs)*b1 + (u + cs)*b2) f2momy = 0.5*(flmomy + frmomy + v*(b1 + b2) + b3) f2momz = 0.5*(flmomz + frmomz + w*(b1 + b2) + b3) # Limit flux to maintain positivity theta = self.limit_flux(state.dens, dens_left, f1dens, f2dens, dtdx) # Limited fluxes fdens = theta*f2dens - (theta - 1.0)*f1dens fmomx = theta*f2momx - (theta - 1.0)*f1momx fmomy = theta*f2momy - (theta - 1.0)*f1momy fmomz = theta*f2momz - (theta - 1.0)*f1momz # Adjust boundary fluxes for closed boundaries if bc[0] == 'closed': c = state.soundspeed[2,:,:] fdens[2,:,:] = 0.0 fmomx[2,:,:] = c*c*state.dens[2,:,:] - 0.5*dx*source[2,:,:] fmomy[2,:,:] = 0.0 fmomz[2,:,:] = 0.0 if bc[1] == 'closed': c = state.soundspeed[-3,:,:] fdens[-2,:,:] = 0.0 fmomx[-2,:,:] = c*c*state.dens[-3,:,:] + 0.5*dx*source[-3,:,:] fmomy[-2,:,:] = 0.0 fmomz[-2,:,:] = 0.0 if bc[0] == 'symmetric': fdens[2,:,:] = 0.0 fmomy[2,:,:] = 0.0 fmomz[2,:,:] = 0.0 if bc[1] == 'symmetric': fdens[-2,:,:] = 0.0 fmomy[-2,:,:] = 0.0 fmomz[-2,:,:] = 0.0 # Interface flux difference: (i+1) - (i) dfdens = np.roll(fdens, -1, axis=0) - fdens dfmomx = np.roll(fmomx, -1, axis=0) - fmomx dfmomy = np.roll(fmomy, -1, axis=0) - fmomy dfmomz = np.roll(fmomz, -1, axis=0) - fmomz # Change in density and momenta (ignoring ghost cells) ddens = -dtdx*dfdens*state.no_ghost dmomx = (dt*source - dtdx*dfmomx)*state.no_ghost dmomy = -dtdx*dfmomy*state.no_ghost dmomz = -dtdx*dfmomz*state.no_ghost state.dens += ddens state.velx += (dmomx - ddens*state.velx)/state.dens state.vely += (dmomy - ddens*state.vely)/state.dens state.velz += (dmomz - ddens*state.velz)/state.dens