Coverage for dragonfly_sim/utils/farfield_utils.py: 100%

1""" 

2Characteristic (Riemann-invariant) far-field state. 

3 

4Used by CompressibleEulerModel.define_farfield_bc. It returns the inflow 

5switch along with the state, so the caller can upwind any further 

6transported scalars (turbulence variables) with it. Unlike a fixed 

7free-stream state, it lets outgoing acoustic waves leave the domain: 

8exact at normal incidence, partly reflecting at oblique incidence (see 

9core/farfield.py for the transverse correction). 

10""" 

11import ufl 

12 

13from dragonfly_sim.utils.Euler_utils import pressure, primitives 

14 

15 

16def riemann_farfield_state(U, n, rho_inf, u_inf, p_inf, gamma, R_in=None): 

17 """ 

18 Subsonic far-field state from the Riemann invariants, and the inflow flag. 

19 

20 With n the outward normal, R+ = u.n + 2c/(g-1) comes from the interior 

21 and R- = u.n - 2c/(g-1) from the free stream (or R_in when given): 

22 

23 u.n_b = (R+ + R-)/2, c_b = (g-1)(R+ - R-)/4 

24 

25 Entropy p/rho^g and the tangential velocity are taken from the free 

26 stream where u.n_b < 0 (inflow) and from the interior otherwise. A 

27 uniform free-stream interior reproduces U_inf exactly. Exact for waves at 

28 normal incidence; oblique waves are partly reflected. 

29 

30 U is the mean flow (rho, rho u, rho E). Returns (U_b, inflow). 

31 """ 

32 rho_i, u_i, _ = primitives(U) 

33 p_i = pressure(U, gamma) 

34 c_i = ufl.sqrt(gamma * p_i / rho_i) 

35 c_inf = (gamma * p_inf / rho_inf) ** 0.5 

36 un_i = ufl.dot(u_i, n) 

37 un_inf = ufl.dot(u_inf, n) 

38 

39 R_out = un_i + 2.0 * c_i / (gamma - 1.0) 

40 if R_in is None: 

41 R_in = un_inf - 2.0 * c_inf / (gamma - 1.0) 

42 un_b = 0.5 * (R_out + R_in) 

43 c_b = 0.25 * (gamma - 1.0) * (R_out - R_in) 

44 

45 inflow = ufl.lt(un_b, 0.0) 

46 s_b = ufl.conditional(inflow, p_inf / rho_inf ** gamma, p_i / rho_i ** gamma) 

47 ut_b = ufl.conditional(inflow, u_inf - un_inf * n, u_i - un_i * n) 

48 

49 rho_b = (c_b ** 2 / (gamma * s_b)) ** (1.0 / (gamma - 1.0)) 

50 p_b = rho_b * c_b ** 2 / gamma 

51 u_b = ut_b + un_b * n 

52 d = u_b.ufl_shape[0] 

53 U_b = ufl.as_vector([rho_b, *[rho_b * u_b[i] for i in range(d)], 

54 p_b / (gamma - 1.0) + 0.5 * rho_b * ufl.dot(u_b, u_b)]) 

55 return U_b, inflow