Coverage for dragonfly_sim/core/MUSCL_Euler_model.py: 86%

1""" 

2Euler equations at p = 0 with MUSCL reconstruction: second-order inviscid 

3fluxes. 

4 

5The plain p = 0 CompressibleEulerModel evaluates the HLL flux on the cell 

6averages, a first-order scheme whose numerical dissipation dominates the drag 

7of a wing (it scales with the wetted area). MUSCLEulerModel installs the same 

8Green-Gauss reconstruction the RANS model uses (core/reconstruction.py) on the 

9Euler model: the interior-facet traces become U_K + G_K . (x_f - x_K), the 

10gradient G is condensed out of Newton exactly, and the adjoint and shape 

11derivatives include the chain terms through G and the cell centres. 

12 

13The solver defaults are the RANS model's: pseudo-transient continuation from 

14the free stream (CFL0 5), one Newton step per pseudo step, FGMRES on the exact 

15condensed Jacobian, preconditioned with ASM(1)/ILU(2) of its compact, 

16gradient-frozen part (jacobian_mode "compact_pc"). 

17 

18Use it through DG_windtunnel_model(model_class=MUSCLEulerModel, 

19model_kwargs={...}). 

20""" 

21from dragonfly_sim.core.Euler_model import CompressibleEulerModel 

22from dragonfly_sim.core.reconstruction import GreenGaussReconstruction 

23from dragonfly_sim.utils.Euler_utils import hll_flux 

24 

25 

26class MUSCLEulerModel(CompressibleEulerModel): 

27 """CompressibleEulerModel with a MUSCL (Green-Gauss) reconstruction. 

28 

29 limiter: None (unlimited MUSCL, the RANS default) or "vanalbada" (smooth, 

30 so the exact Jacobian survives; see GreenGaussReconstruction). 

31 limiter_eps: the van Albada smoothing parameter, relative to the 

32 free-stream scale of each state component. 

33 jacobian_mode: what Newton's linear solve uses (utils/gradient_condensation.py). 

34 use_ptc, cfl0: pseudo-transient continuation from the free stream. 

35 """ 

36 

37 def __init__(self, mesh, poly_order=0, gamma=1.4, limiter=None, limiter_eps=1e-2, 

38 jacobian_mode="compact_pc", use_ptc=True, cfl0=5.0): 

39 super().__init__(mesh, poly_order=poly_order, gamma=gamma, flux_function=hll_flux) 

40 self.limiter = limiter 

41 self.limiter_eps = limiter_eps 

42 self.jacobian_mode = jacobian_mode 

43 self.reconstruction = None 

44 # solver defaults (see the module docstring) 

45 self.use_ptc = use_ptc 

46 self.ptc_settings = {"cfl0": cfl0} 

47 if use_ptc: 

48 self.fom_newton_inner_max_it = 1 

49 self.krylov_type = "fgmres" 

50 self.asm_overlap = 1 

51 self.ilu_levels = 2 

52 

53 def define_elements_functionspaces(self): 

54 super().define_elements_functionspaces() 

55 if self.reconstruction is None: 

56 GreenGaussReconstruction(self, muscl=True, limiter=self.limiter, 

57 limiter_eps=self.limiter_eps, 

58 limiter_scales=self._limiter_scales).install() 

59 

60 def _limiter_scales(self): 

61 """Free-stream magnitude of each state component (van Albada limiter).""" 

62 ff = self.farfield_conditions 

63 rho, p = ff['rho'], ff['p'] 

64 speed = ff['M'] * ff['c'] 

65 return ([rho] + [rho * speed] * self.dimensions 

66 + [p / (self.gamma - 1.0) + 0.5 * rho * speed**2]) 

67 

68 def set_up_solver(self): 

69 was_set_up = self.solver_is_set_up 

70 super().set_up_solver() 

71 if was_set_up: 

72 return 

73 mode = "exact" if self.direct_linear_solver else self.jacobian_mode 

74 if mode != "exact": 

75 self.problem.configure_snes(self.solver.snes, mode) 

76 self.solver.krylov_solver.setFromOptions()