dragonfly_sim.core.reconstruction ================================= .. py:module:: dragonfly_sim.core.reconstruction .. autoapi-nested-parse:: Green-Gauss gradient reconstruction for the p = 0 discretization. At p = 0, grad U vanishes inside every cell, which leaves no viscous volume term and no vorticity for the turbulence sources, and gives only first-order face states. As a cell-centred finite-volume code does, a DG0 gradient G is reconstructed by Green-Gauss over each cell's facets, with facet averages {U} inside and the registered boundary state U_b on the boundary, |K| G_K = sum_f |f| U_f (x) n_f. G is used in three places: * MUSCL: the interior-facet traces of the inviscid flux are reconstructed linearly, U_K + G_K . (x_f - x_K), optionally with a smooth van Albada limiter (see muscl_traces). * The viscous facet gradient, {G} with its component along the line between the two cell centres replaced by the direct difference (which keeps the two-point coupling that suppresses odd-even modes), or on a boundary facet with the normal component replaced by (U_b - U)/d_n. * The turbulence sources (vorticity, gradients of the transported scalars). G is eliminated from the unknowns exactly (utils/gradient_condensation.py), so Newton sees the true Jacobian. For adjoints and shape derivatives G and the cell centres are DerivedFields (utils/derived_fields.py), with exact derivatives with respect to U, alpha and the mesh coordinates. .. !! processed by numpydoc !! Attributes ---------- .. autoapisummary:: dragonfly_sim.core.reconstruction._JIT Classes ------- .. autoapisummary:: dragonfly_sim.core.reconstruction.CellCentres dragonfly_sim.core.reconstruction.GreenGaussGradient dragonfly_sim.core.reconstruction.GreenGaussReconstruction dragonfly_sim.core.reconstruction.GreenGaussStateOperator Module Contents --------------- .. py:class:: CellCentres(mesh_obj) Bases: :py:obj:`dragonfly_sim.utils.derived_fields.DerivedField` DG0 cell centres, each cell's reference midpoint mapped to the cell. For the degree-1 geometry used here that is the plain average of the cell's vertices, c = C x, so dc/dx is the constant averaging matrix C. .. !! processed by numpydoc !! .. py:method:: _mesh_jacobian_loop() .. py:method:: _mesh_jacobian_vectorized() The same matrix as _mesh_jacobian_loop, as one CSR build: row (cell, component i) holds 1/nv at the cell's vertex dofs' component i. .. !! processed by numpydoc !! .. py:method:: discard(obj) Free a matrix or vector returned by state_jacobian / alpha_vector / mesh_jacobian once the caller is done with it. A field that caches them keeps its own (override); a python-context operator also frees what its context holds. .. !! processed by numpydoc !! .. py:method:: mesh_jacobian() .. py:method:: release_cache() Free whatever the field caches between calls (nothing by default). .. !! processed by numpydoc !! .. py:method:: update_geometry() .. py:attribute:: _mesh_jac :value: None .. py:attribute:: function .. py:attribute:: mesh_obj .. py:attribute:: rec :value: None .. py:class:: GreenGaussGradient(reconstruction) Bases: :py:obj:`dragonfly_sim.utils.derived_fields.DerivedField` G = M^-1 b(U, alpha, x), as a DerivedField (see the module docstring). Options (attributes of the reconstruction): cache_jacobians keep dG/dU, dG/dalpha and dG/dx between calls while the state, alpha and the mesh are unchanged (several functionals and products of one totals evaluation share them); freed by release_cache structured_state_jacobian dG/dU as the operator W (x) I plus a small boundary part (GreenGaussStateOperator) instead of an assembled matrix .. !! processed by numpydoc !! .. py:method:: _build_alpha_vector() .. py:method:: _build_mesh_jacobian() .. py:method:: _build_state_jacobian() .. py:method:: _cached(name, build) build() once per (state, alpha, mesh), or every call if caching is off. .. !! processed by numpydoc !! .. py:method:: _form(key, factory) .. py:method:: _key() .. py:method:: _same_key(key) .. py:method:: alpha_vector() .. py:method:: discard(obj) Free a matrix or vector returned by state_jacobian / alpha_vector / mesh_jacobian once the caller is done with it. A field that caches them keeps its own (override); a python-context operator also frees what its context holds. .. !! processed by numpydoc !! .. py:method:: mesh_jacobian() .. py:method:: release_cache() Free whatever the field caches between calls (nothing by default). .. !! processed by numpydoc !! .. py:method:: state_jacobian() .. py:method:: update() .. py:method:: update_geometry() .. py:attribute:: _cache :value: None .. py:attribute:: _forms .. py:attribute:: depends_on_state :value: True .. py:attribute:: function .. py:attribute:: rec .. py:class:: GreenGaussReconstruction(model, muscl=True, limiter=None, limiter_eps=0.01, limiter_scales=None) p = 0 gradient reconstruction attached to a CompressibleEulerModel host. install() plugs it in: the MUSCL face traces (if muscl), the derived fields (cell centres and G), and the condensed nonlinear problem. Call it after the function spaces exist and before the boundary conditions are defined; the Green-Gauss forms are built lazily, once every boundary condition has registered its state. limiter: None (unlimited MUSCL) or "vanalbada". limiter_scales() returns one magnitude per state component for the limiter's smoothing parameter. .. !! processed by numpydoc !! .. py:method:: _weights() W, the scalar operator of the interior-facet part of dG/dU before M_G^-1: (3 n_cells) x n_cells, assembled once per mesh. .. !! processed by numpydoc !! .. py:method:: compact_forms(F) (a_UU, a_UG): the compiled forms of dR_U/dU and dR_U/dG, cached for F. .. !! processed by numpydoc !! .. py:method:: condensed_jacobian(F) The exact dR/dU of the residual form F, as a CondensedJacobian. .. !! processed by numpydoc !! .. py:method:: face_gradient_boundary(U_b, U=None) The cell's G with its normal component replaced by (U_b - U)/d_n, d_n the cell centre's normal distance from the facet. .. !! processed by numpydoc !! .. py:method:: face_gradient_interior() .. py:method:: green_gauss_rhs() b(U; tau) with |K| G_K = b: facet averages inside, the registered boundary states on the boundary. .. !! processed by numpydoc !! .. py:method:: install() .. py:method:: make_problem(F_newton, J_newton) .. py:method:: muscl_traces() Linear reconstruction of both interior-facet traces from G at the facet's quadrature points: U_K + G_K . (x_f - x_K) unlimited, or with limiter="vanalbada" the edge-based van Albada form (Blazek, Computational Fluid Dynamics, ch. 5) generalised to an off-centre facet. For the cell K and its neighbour N across the facet, with s = |x_f - x_K| / |x_N - x_K|, D+ = s (U_N - U_K) (central increment to the facet) D- = 2 G_K . (x_f - x_K) - D+ (upwind increment) U_f = U_K + Psi(D-, D+), Psi(a, b) = ((a^2 + e) b + (b^2 + e) a) / (a^2 + b^2 + 2e), componentwise. Psi(a, a) = a, so smooth regions keep the unlimited reconstruction; at a shock D- ~ 0 and the trace falls back to first order. Psi is smooth, so the exact Jacobian survives. e = (limiter_eps * S_k)^2 with S_k the component's scale. .. !! processed by numpydoc !! .. py:method:: release_cache() .. py:method:: set_up() .. py:method:: state_operator() dG/dU at the current state, as a python-context Mat. .. !! processed by numpydoc !! .. py:method:: update_gradient() G <- M_G^-1 b(U) at the current state and mesh. Collective. .. !! processed by numpydoc !! .. py:method:: update_mass() Refresh M_G^-1 in place after mesh motion. .. !! processed by numpydoc !! .. py:attribute:: G .. py:attribute:: _W :value: None .. py:attribute:: _W_key :value: None .. py:attribute:: _gg_form :value: None .. py:attribute:: cache_jacobians :value: False .. py:attribute:: centres .. py:attribute:: fast_cell_centre_jacobian :value: False .. py:attribute:: gg_rhs :value: None .. py:attribute:: gradient .. py:attribute:: inv_mass :value: None .. py:attribute:: limiter :value: None .. py:attribute:: limiter_eps :value: 0.01 .. py:attribute:: limiter_scales :value: None .. py:attribute:: mesh_obj .. py:attribute:: model .. py:attribute:: muscl :value: True .. py:attribute:: structured_state_jacobian :value: False .. py:attribute:: xc .. py:class:: GreenGaussStateOperator(rec, W, B_bdry) dG/dU = M_G^-1 db/dU as a python-context Mat (G layout x state layout) with mult and multTranspose, without assembling the 18-by-6 blocks. The interior-facet part of b is a facet average of U times the area vector, component by component: its derivative is W (x) I_n, with W the scalar operator (3 area-vector components per cell pair, assembled once per mesh) and I_n the identity over the n state components. Only the boundary-facet part, U_b(U) of the registered boundary states, couples components; it is assembled per state (boundary cells only). .. !! processed by numpydoc !! .. py:method:: destroy(mat=None) .. py:method:: mult(mat, x, y) .. py:method:: multTranspose(mat, v, y) .. py:attribute:: _t .. py:attribute:: n .. py:data:: _JIT