submodule (mod_matrix_manager) smod_flow_matrix implicit none contains module procedure add_flow_matrix_terms real(dp) :: eps, deps real(dp) :: rho, drho real(dp) :: T0 real(dp) :: v01, dv01, drv01 real(dp) :: v02, dv02, drv02 real(dp) :: v03, dv03 real(dp) :: Vop real(dp) :: gamma_1 gamma_1 = settings%physics%get_gamma_1() ! grid variables eps = grid%get_eps(x) deps = grid%get_deps() ! density variables rho = background%density%rho0(x) drho = background%density%drho0(x) ! temperature variables T0 = background%temperature%T0(x) ! flow variables v01 = background%velocity%v01(x) dv01 = background%velocity%dv01(x) drv01 = deps * v01 + eps * dv01 v02 = background%velocity%v02(x) dv02 = background%velocity%dv02(x) drv02 = deps * v02 + eps * dv02 v03 = background%velocity%v03(x) dv03 = background%velocity%dv03(x) Vop = k2 * v02 / eps + k3 * v03 ! ==================== Quadratic * Quadratic ==================== call elements%add(Vop - ic * dv01, sv_rho1, sv_rho1) call elements%add( & -drv02 * ic * v01 / eps, sv_v2, sv_rho1) call elements%add( & rho * (eps * Vop - ic * deps * v01), sv_v2, sv_v2) call elements%add(-ic * v01 * dv03, sv_v3, sv_rho1) call elements%add( & rho * (Vop + ic * dv01) + (deps * rho / eps + drho) * ic * v01, sv_v3, sv_v3 & ) call elements%add(eps * Vop, sv_a1, sv_a1) ! ==================== Quadratic * dQuadratic ==================== call elements%add(-ic * v01, sv_rho1, sv_rho1, s2do=1) call elements%add(-ic * eps * rho * v01, sv_v2, sv_v2, s2do=1) ! ==================== Cubic * Quadratic ==================== call elements%add(v01 * dv01 - deps * v02**2 / eps, sv_v1, sv_rho1) call elements%add(-2.0d0 * deps * rho * v02, sv_v1, sv_v2) call elements%add(ic * v01 * k2, sv_a2, sv_a1) call elements%add(eps * k3 * ic * v01, sv_a3, sv_a1) ! ==================== Cubic * Cubic ==================== call elements%add(rho * Vop + (deps * rho / eps + drho) * ic * v01, sv_v1, sv_v1) call elements%add(k3 * v03, sv_a2, sv_a2) call elements%add(-k2 * v03, sv_a2, sv_a3) call elements%add(-k3 * v02, sv_a3, sv_a2) call elements%add(k2 * v02, sv_a3, sv_a3) ! ==================== dCubic * Cubic ==================== call elements%add(ic * rho * v01, sv_v1, sv_v1, s1do=1) ! ==================== Quadratic * Cubic ==================== call elements%add(-drv02 * rho / eps, sv_v2, sv_v1) call elements%add(-rho * dv03, sv_v3, sv_v1) ! ==================== dQuadratic * Quadratic ==================== call elements%add(ic * rho * v01, sv_v3, sv_v3, s1do=1) ! ==================== Quadratic * dCubic ==================== call elements%add(-v02, sv_a1, sv_a2, s2do=1) call elements%add(-eps * v03, sv_a1, sv_a3, s2do=1) ! ==================== Cubic * dCubic ==================== call elements%add(-ic * v01, sv_a2, sv_a2, s2do=1) call elements%add(-ic * v01, sv_a3, sv_a3, s2do=1) if (.not. settings%physics%is_incompressible) then ! ==================== Quadratic * Quadratic ==================== call elements%add(-ic * gamma_1 * drv01 * T0 / eps, sv_T1, sv_rho1) call elements%add( & rho * (Vop + ic * dv01 - ic * gamma_1 * drv01 / eps) & + ic * v01 * (deps * rho / eps + drho), & sv_T1, & sv_T1 & ) ! ==================== dQuadratic * Quadratic ==================== call elements%add(ic * rho * v01, sv_T1, sv_T1, s1do=1) end if end procedure add_flow_matrix_terms end submodule smod_flow_matrix