submodule (mod_matrix_manager) smod_conduction_matrix implicit none contains module procedure add_conduction_matrix_terms real(dp) :: eps, deps real(dp) :: dT0 real(dp) :: WVop real(dp) :: kappa_perp, dkappa_perp_drho, dkappa_perp_dT real(dp) :: gamma_1 type(matrix_elements_t) :: elements if (settings%physics%is_incompressible) return gamma_1 = settings%physics%get_gamma_1() eps = grid%get_eps(x_gauss) deps = grid%get_deps() dT0 = background%temperature%dT0(x_gauss) kappa_perp = physics%conduction%tcperp(x_gauss) dkappa_perp_drho = physics%conduction%dtcperpdrho(x_gauss) dkappa_perp_dT = physics%conduction%dtcperpdT(x_gauss) ! operators WVop = k2**2 / eps + eps * k3**2 elements = new_matrix_elements(state_vector=settings%get_state_vector()) ! ==================== Quadratic * Quadratic ==================== call elements%add( & -gamma_1 * ic * WVop * kappa_perp / eps, & "T", & "T", & spline1=h_quad, & spline2=h_quad & ) ! ==================== dQuadratic * Quadratic ==================== call elements%add( & -ic * gamma_1 * dT0 * dkappa_perp_drho, & "T", & "rho", & spline1=dh_quad, & spline2=h_quad & ) call elements%add( & gamma_1 * (ic * deps * kappa_perp / eps - ic * dT0 * dkappa_perp_dT), & "T", & "T", & spline1=dh_quad, & spline2=h_quad & ) ! ==================== dQuadratic * dQuadratic ==================== call elements%add( & -ic * gamma_1 * kappa_perp, "T", "T", spline1=dh_quad, spline2=dh_quad & ) if (settings%has_bfield()) then call add_conduction_matrix_terms_bfield( & x_gauss, settings, grid, background, physics, elements & ) end if call add_to_quadblock(quadblock, elements, weight, settings%dims) call elements%delete() end procedure add_conduction_matrix_terms subroutine add_conduction_matrix_terms_bfield( & x_gauss, settings, grid, background, physics, elements & ) real(dp), intent(in) :: x_gauss type(settings_t), intent(in) :: settings type(grid_t), intent(in) :: grid type(background_t), intent(in) :: background type(physics_t), intent(in) :: physics type(matrix_elements_t), intent(inout) :: elements real(dp) :: eps, deps real(dp) :: dT0, ddT0 real(dp) :: B0, B01, B02, B03, dB02, dB03 real(dp) :: diffKp, Kp, Kp_plus, Kp_plusplus real(dp) :: Fop_plus, dFop_plus, Gop_min real(dp) :: kappa_para, dkappa_para_dT real(dp) :: kappa_perp, dkappa_perp_drho, dkappa_perp_dT, dkappa_perp_dB2 real(dp) :: gamma_1 complex(dp) :: Fop_B01 gamma_1 = settings%physics%get_gamma_1() ! grid variables eps = grid%get_eps(x_gauss) deps = grid%get_deps() ! temperature variables dT0 = background%temperature%dT0(x_gauss) ddT0 = background%temperature%ddT0(x_gauss) ! magnetic field variables B0 = background%magnetic%get_B0(x_gauss) B01 = background%magnetic%B01(x_gauss) B02 = background%magnetic%B02(x_gauss) dB02 = background%magnetic%dB02(x_gauss) B03 = background%magnetic%B03(x_gauss) dB03 = background%magnetic%dB03(x_gauss) ! parallel thermal conduction variables kappa_para = physics%conduction%tcpara(x_gauss) dkappa_para_dT = physics%conduction%dtcparadT(x_gauss) ! perpendicular thermal conduction variables kappa_perp = physics%conduction%tcperp(x_gauss) dkappa_perp_drho = physics%conduction%dtcperpdrho(x_gauss) dkappa_perp_dT = physics%conduction%dtcperpdT(x_gauss) dkappa_perp_dB2 = physics%conduction%dtcperpdB2(x_gauss) ! prefactors Kp = physics%conduction%get_tcprefactor(x_gauss) diffKp = physics%conduction%get_dtcprefactordr(x_gauss) Kp_plus = Kp + dkappa_perp_dB2 Kp_plusplus = dkappa_perp_dB2 - (B01**2 * Kp_plus / B0**2) ! operators Fop_plus = k2 * B02 / eps + k3 * B03 dFop_plus = (k2 / eps) * (dB02 - deps * B02 / eps) + k3 * dB03 Gop_min = k3 * B02 - k2 * B03 / eps ! B01 modified F+ operator Fop_B01 = deps * ic * B01 / eps + Fop_plus ! ==================== Quadratic * Quadratic ==================== call elements%add( & gamma_1 * dT0 * (B01 / B0**2) * dkappa_perp_drho * Fop_B01, & "T", & "rho", & spline1=h_quad, & spline2=h_quad & ) call elements%add( & gamma_1 * ( & B01 * Fop_plus * diffKp & - B01 * dT0 * (dkappa_para_dT - dkappa_perp_dT) * Fop_B01 / B0**2 & + Kp * ( & B01 * deps * (2.0d0 * deps * ic * B01 / eps + 3.0d0 * Fop_plus) / eps & + B01 * dFop_plus & - ic * Fop_plus**2 & ) & ), & "T", & "T", & spline1=h_quad, & spline2=h_quad & ) call elements%add( & 2.0d0 * gamma_1 * eps * dT0 * Gop_min * Kp_plus * B01 * Fop_B01 / B0**2, & "T", & "a1", & spline1=h_quad, & spline2=h_quad & ) ! ==================== dQuadratic * Quadratic ==================== call elements%add( & ic * gamma_1 * dT0 * dkappa_perp_drho * B01**2 / B0**2, & "T", & "rho", & spline1=dh_quad, & spline2=h_quad & ) call elements%add( & gamma_1 * ( & B01 * Kp * (2.0d0 * deps * ic * B01 / eps + 3.0d0 * Fop_plus) & - ic * dT0 * ( & B01**2 * dkappa_para_dT / B0**2 - dkappa_perp_dT * B01**2 / B0**2 & ) & ), & "T", & "T", & spline1=dh_quad, & spline2=h_quad & ) call elements%add( & -2.0d0 * ic * gamma_1 * eps * dT0 * Gop_min * Kp_plusplus, & "T", & "a1", & spline1=dh_quad, & spline2=h_quad & ) ! ==================== Quadratic * dQuadratic ==================== call elements%add( & -2.0d0 * ic * gamma_1 * deps * B01**2 * Kp / eps, & "T", & "T", & spline1=h_quad, & spline2=dh_quad & ) ! ==================== dQuadratic * dQuadratic ==================== call elements%add( & -ic * gamma_1 * 2.0d0 * B01**2 * Kp, & "T", & "T", & spline1=dh_quad, & spline2=dh_quad & ) ! ==================== Quadratic * Cubic ==================== call elements%add( & gamma_1 * k3 * ( & Kp * (B01 * ddT0 + ic * dT0 * Fop_plus) & + B01 * dT0 * (diffKp - 2.0d0 * ic * B01 * Kp_plus * Fop_B01 / B0**2) & ), & "T", & "a2", & spline1=h_quad, & spline2=h_cubic & ) call elements%add( & -gamma_1 * k2 * ( & Kp * (B01 * ddT0 + ic * dT0 * Fop_plus) & + B01 * dT0 * (diffKp - 2.0d0 * ic * B01 * Kp_plus * Fop_B01 / B0**2) & ), & "T", & "a3", & spline1=h_quad, & spline2=h_cubic & ) ! ==================== Quadratic * dCubic ==================== call elements%add( & gamma_1 * dT0 * B01 * (2.0d0 * B03 * Kp_plus * Fop_B01 / B0**2 - k3 * Kp), & "T", & "a2", & spline1=h_quad, & spline2=dh_cubic & ) call elements%add( & gamma_1 * dT0 * B01 * ( & -2.0d0 * eps * B02 * Kp_plus * Fop_B01 / B0**2 + k2 * Kp & ), & "T", & "a3", & spline1=h_quad, & spline2=dh_cubic & ) ! ==================== dQuadratic * Cubic ==================== call elements%add( & -2.0d0 * gamma_1 * dT0 * B01 * k3 * Kp_plusplus, & "T", & "a2", & spline1=dh_quad, & spline2=h_cubic & ) call elements%add( & 2.0d0 * gamma_1 * dT0 * B01 * k2 * Kp_plusplus, & "T", & "a3", & spline1=dh_quad, & spline2=h_cubic & ) ! ==================== dQuadratic * dCubic ==================== call elements%add( & -2.0d0 * ic * gamma_1 * dT0 * B03 * Kp_plusplus, & "T", & "a2", & spline1=dh_quad, & spline2=dh_cubic & ) call elements%add( & 2.0d0 * ic * gamma_1 * dT0 * eps * B02 * Kp_plusplus, & "T", & "a3", & spline1=dh_quad, & spline2=dh_cubic & ) end subroutine add_conduction_matrix_terms_bfield end submodule smod_conduction_matrix