From f924530319d02bbd11a3e0939faf3ebbb39d972e Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Thu, 16 Sep 2021 11:04:32 +0200 Subject: [PATCH 1/9] remove old modules (reworked/moved) --- src/CMakeLists.txt | 2 - src/mod_boundary_conditions.f08 | 282 -------------- src/mod_hallmhd.f08 | 632 -------------------------------- 3 files changed, 916 deletions(-) delete mode 100644 src/mod_boundary_conditions.f08 delete mode 100644 src/mod_hallmhd.f08 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ed59885e..248fe69a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -73,8 +73,6 @@ set (sources matrices/smod_conduction_matrix.f08 matrices/smod_viscosity_matrix.f08 matrices/smod_hall_matrix.f08 - # physics/mod_viscosity.f08 - # mod_boundary_conditions.f08 solvers/mod_solvers.f08 solvers/mod_arpack_type.f08 solvers/smod_qr_invert.f08 diff --git a/src/mod_boundary_conditions.f08 b/src/mod_boundary_conditions.f08 deleted file mode 100644 index c2d49498..00000000 --- a/src/mod_boundary_conditions.f08 +++ /dev/null @@ -1,282 +0,0 @@ -! ============================================================================= -!> Imposes boundary conditions on the finite element matrices after assembly. -!! Both natural and essential conditions are handled. -module mod_boundary_conditions - use mod_global_variables, only: dp, matrix_gridpts, dim_quadblock, dim_subblock - implicit none - - private - - !> boolean to check presence of perpendicular thermal conduction - logical, save :: kappa_perp_is_zero - - public :: apply_boundary_conditions - -contains - - - !> Main routine to handle the boundary conditions. This subroutine first performs some checks - !! and then handles the natural and essential boundary conditions. - !! @note Perpendicular thermal conduction is hard-checked, that is, not just the - !! global logical. There must be a value in the corresponding array that is non-zero. - subroutine apply_boundary_conditions(matrix_A, matrix_B) - use mod_equilibrium, only: kappa_field - ! use mod_viscosity, only: viscosity_boundaries - use mod_global_variables, only: dp_LIMIT - use mod_equilibrium, only: kappa_field - - !> the A-matrix with boundary conditions imposed on exit - complex(dp), intent(inout) :: matrix_A(matrix_gridpts, matrix_gridpts) - !> the B-matrix with boundary conditions imposed on exit - real(dp), intent(inout) :: matrix_B(matrix_gridpts, matrix_gridpts) - complex(dp) :: quadblock(dim_quadblock, dim_quadblock) - ! complex(dp) :: quadblock_visc(dim_quadblock, dim_quadblock) - ! complex(dp) :: quadblock_Hall(dim_quadblock, dim_quadblock) - integer :: idx_end_left, idx_start_right - ! real(dp) :: hf - - ! check if perpendicular thermal conduction is present - kappa_perp_is_zero = .true. - if (any(abs(kappa_field % kappa_perp) > dp_LIMIT)) then - kappa_perp_is_zero = .false. - end if - - ! end of index first-gridpoint quadblock - idx_end_left = dim_quadblock - ! start of index last-gridpoint quadblock - idx_start_right = matrix_gridpts - dim_quadblock + 1 - - ! matrix B left-edge quadblock - quadblock = matrix_B(1:idx_end_left, 1:idx_end_left) - call essential_boundaries(quadblock, edge='l_edge', matrix='B') - matrix_B(1:idx_end_left, 1:idx_end_left) = real(quadblock) - ! matrix B right-edge quadblock - quadblock = matrix_B(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) - call essential_boundaries(quadblock, edge='r_edge', matrix='B') - matrix_B(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) = real(quadblock) - - ! matrix A left-edge quadblock - quadblock = matrix_A(1:idx_end_left, 1:idx_end_left) - call essential_boundaries(quadblock, edge='l_edge', matrix='A') - call natural_boundaries(quadblock, edge='l_edge') - ! if (viscosity) then - ! call viscosity_boundaries(quadblock_visc, edge='l_edge') - ! quadblock = quadblock + quadblock_visc - ! end if - ! if (hall_mhd) then - ! call hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge='l_edge') - ! hf = hall_field % hallfactor(1) - ! quadblock = quadblock - hf * quadblock_Hall - ! end if - matrix_A(1:idx_end_left, 1:idx_end_left) = quadblock - ! matrix A right-edge quadblock - quadblock = matrix_A(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) - call essential_boundaries(quadblock, edge='r_edge', matrix='A') - call natural_boundaries(quadblock, edge='r_edge') - ! if (hall_mhd) then - ! call hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge='r_edge') - ! hf = hall_field % hallfactor(gauss_gridpts) - ! quadblock = quadblock - hf * quadblock_Hall - ! end if - matrix_A(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) = quadblock - end subroutine apply_boundary_conditions - - - !> Imposes essential boundary conditions on a given matrix, that is, - !! the ones that have to be implemented explicitly. - !! For a wall, that means handling \(v1, a2\) and \(a3\) (and \(T1\) if there is - !! perpendicular thermal conduction). - !! @note The contribution from the 0 quadratic basis function automatically - !! zeroes out the odd rows/columns for the quadratic variables on the - !! left edge. These are explicitly handled. @endnote - !! @note In the case of the Arnoldi solver, setting boundary conditions to - !! 1e20 throws off the solver such that it returns completely wrong results. - !! Setting them to a "normal" value, 0 here, seems to fix it. We now have - !! an exact match with the QR/QZ solvers in this case. - !! Eigenvalues will be introduced at 0 instead of at 1, but this does not - !! affect the final spectrum in any way. @endnote - !! @warning Throws an error if boundary_type is not known, - !! or if edge is not known. - subroutine essential_boundaries(quadblock, edge, matrix) - use mod_global_variables, only: boundary_type, solver, viscosity, geometry, & - coaxial - use mod_logging, only: log_message - - !> the quadblock corresponding to the left/right edge - complex(dp), intent(inout) :: quadblock(dim_quadblock, dim_quadblock) - !> the edge, either "l_edge" or "r_edge" - character(len=6), intent(in) :: edge - !> the matrix, either "A" or "B" - character, intent(in) :: matrix - - complex(dp) :: diagonal_factor - integer :: i, j, qua_zeroes(5), wall_idx_left(4), wall_idx_right(4) - integer :: noslip_idx_left(2), noslip_idx_right(2) - - if (matrix == 'B') then - diagonal_factor = (1.0d0, 0.0d0) - else if (matrix == 'A') then - if (solver == "arnoldi") then - diagonal_factor = (0.0d0, 0.0d0) - else - diagonal_factor = (1.0d20, 0.0d0) - end if - else - call log_message("essential boundaries: invalid matrix argument", level='error') - end if - - ! explicitly handle zeroed out rows/cols due to 0 term quadratic basis function - qua_zeroes = [1, 5, 7, 9, 11] - if (edge == 'l_edge') then - do i = 1, size(qua_zeroes) - j = qua_zeroes(i) - quadblock(j, j) = diagonal_factor - end do - end if - - ! wall/regularity conditions: handling of v1, a2 and a3 (and T if conduction). - ! v1, a2 and a3 are cubic elements, so omit non-zero basis functions (odd rows/columns) - ! T is a quadratic element, so omit even row/columns - wall_idx_left = [3, 13, 15, 10] - wall_idx_right = [19, 29, 31, 26] - ! For a no-slip condition for viscous fluids v2 and v3 should equal the wall's - ! tangential velocities, here zero. (Quadratic elements, so omit even row/columns) - noslip_idx_left = [6, 8] - noslip_idx_right = [22, 24] - - select case(boundary_type) - case('wall') - if (edge == 'l_edge') then - ! left regularity/wall conditions - do i = 1, size(wall_idx_left) - j = wall_idx_left(i) - if (j == 10 .and. kappa_perp_is_zero) then - cycle - end if - quadblock(j, :) = (0.0d0, 0.0d0) - quadblock(:, j) = (0.0d0, 0.0d0) - quadblock(j, j) = diagonal_factor - end do - ! No-slip condition does not apply at left edge for a cylindrical - ! geometry unless two coaxial walls are used - if (viscosity .and. (geometry == 'Cartesian' .or. coaxial)) then - do i = 1, size(noslip_idx_left) - j = noslip_idx_left(i) - quadblock(j, :) = (0.0d0, 0.0d0) - quadblock(:, j) = (0.0d0, 0.0d0) - quadblock(j, j) = diagonal_factor - end do - end if - else if (edge == 'r_edge') then - do i = 1, size(wall_idx_right) - j = wall_idx_right(i) - if ((j == 26) .and. kappa_perp_is_zero) then - cycle - end if - quadblock(j, :) = (0.0d0, 0.0d0) - quadblock(:, j) = (0.0d0, 0.0d0) - quadblock(j, j) = diagonal_factor - end do - if (viscosity) then - do i = 1, size(noslip_idx_right) - j = noslip_idx_right(i) - quadblock(j, :) = (0.0d0, 0.0d0) - quadblock(:, j) = (0.0d0, 0.0d0) - quadblock(j, j) = diagonal_factor - end do - end if - else - call log_message("essential boundaries: invalid edge argument", level='error') - end if - - case default - call log_message( "essential boundaries: invalid boundary_type", level='error') - end select - end subroutine essential_boundaries - - - !> Imposes the natural boundary conditions, that is, the ones originating - !! from the weak formulation. This is only applicable for the A-matrix. - !! @note For now only the terms of a solid wall boundary are implemented, hence - !! we only have additional resistive and conductive terms in the energy equation. - !! If perpendicular thermal conduction is included we have \(T1 = 0\), such that - !! then there are no surface terms to be added and we simply return. - !! If free boundary conditions are considered (e.g. vacuum-wall) additional terms - !! have to be added. @endnote - !! @warning Throws an error if edge is unknown. - subroutine natural_boundaries(quadblock, edge) - use mod_global_variables, only: boundary_type, gauss_gridpts, gamma_1, ic, dim_subblock - use mod_logging, only: log_message - use mod_equilibrium, only: B_field, eta_field - use mod_grid, only: eps_grid, d_eps_grid_dr - use mod_equilibrium_params, only: k2, k3 - - !> the quadblock corresponding to the left/right edge - complex(dp), intent(inout) :: quadblock(dim_quadblock, dim_quadblock) - !> the edge, either "l_edge" or "r_edge" - character(len=6), intent(in) :: edge - - complex(dp), allocatable :: surface_terms(:) - real(dp) :: eps, d_eps_dr, eta, B02, dB02, B03, dB03, drB02 - integer, allocatable :: positions(:, :) - integer :: idx, i - - if (.not. boundary_type == 'wall') then - call log_message('natural boundaries: only wall is implemented!', level='error') - end if - - ! For a wall with perpendicular thermal conduction we also have the condition T1 = 0. - ! Hence, in that case there are no surface terms to be added, and we simply return. - ! So for now we only have resistive terms in the energy equation. - if (.not. kappa_perp_is_zero) then - return - end if - - if (edge == 'l_edge') then - idx = 1 - else if (edge == 'r_edge') then - idx = gauss_gridpts - else - call log_message('natural boundaries: wrong edge supplied' // edge, level='error') - end if - - ! retrieve variables at current edge - eps = eps_grid(idx) - d_eps_dr = d_eps_grid_dr(idx) - B02 = B_field % B02(idx) - dB02 = B_field % d_B02_dr(idx) - B03 = B_field % B03(idx) - dB03 = B_field % d_B03_dr(idx) - drB02 = d_eps_dr * B02 + eps * dB02 - eta = eta_field % eta(idx) - - allocate(positions(3, 2)) - allocate(surface_terms(3)) - - ! surface term for element (5, 6) - surface_terms(1) = 2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * (k3 * drB02 - k2 * dB03) - positions(1, :) = [5, 6] - ! surface term for element (5, 7) - surface_terms(2) = 2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * dB03 - positions(2, :) = [5, 7] - ! surface term for element (5, 8) - surface_terms(3) = -2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * drB02 - positions(3, :) = [5, 8] - - ! l_edge: add to bottom-right of 2x2 block, for top-left subblock only - ! r_edge: add to bottom-right of 2x2 block, for bottom-right subblock only - if (edge == 'l_edge') then - positions = 2 * positions - else if (edge == 'r_edge') then - positions = 2 * positions + dim_subblock - end if - - do i = 1, size(surface_terms) - quadblock(positions(i, 1), positions(i, 2)) = quadblock(positions(i, 1), positions(i, 2)) + surface_terms(i) - end do - - deallocate(positions) - deallocate(surface_terms) - end subroutine natural_boundaries - -end module mod_boundary_conditions diff --git a/src/mod_hallmhd.f08 b/src/mod_hallmhd.f08 deleted file mode 100644 index ea1ee552..00000000 --- a/src/mod_hallmhd.f08 +++ /dev/null @@ -1,632 +0,0 @@ -! ============================================================================= -!> This module handles everything related to the Hall term, computing -!! the normalised Hall factor, the contributions to the A matrix, and natural -!! boundary conditions, when necessary. -module mod_hallmhd - use mod_global_variables, only: dp, dim_quadblock - -implicit none - -private - -public :: set_hallfactor -public :: get_hallterm -public :: get_hallterm_Bmat -public :: hall_boundaries - -contains - - !> Retrieves the normalised Hall factor as described by Porth et al. (2014), - !! with a dropoff at the boundary - subroutine set_hallfactor(hall_field) - use mod_grid, only: grid_gauss - use mod_physical_constants, only: dpi - use mod_global_variables, only: cgs_units, gauss_gridpts, dropoff_edge_dist, & - dropoff_width, hall_dropoff, inertia_dropoff - use mod_physical_constants, only: mp_cgs, mp_si, ec_cgs, ec_si, me_cgs, me_si - use mod_units, only: unit_velocity, unit_length, unit_magneticfield - use mod_types, only: hall_type - - type (hall_type), intent(inout) :: hall_field - - real(dp) :: sleft, sright, width, hallval, inertiaval - real(dp) :: x, shift, stretch, shift2, stretch2 - integer :: i - - width = dropoff_width - if (cgs_units) then - hallval = (mp_cgs * unit_velocity) / (ec_cgs * unit_length * unit_magneticfield) - inertiaval = (mp_cgs * me_cgs * unit_velocity**2) / (ec_cgs * unit_length * unit_magneticfield)**2 - else - hallval = (mp_si * unit_velocity) / (ec_si * unit_length * unit_magneticfield) - inertiaval = (mp_si * me_si * unit_velocity**2) / (ec_si * unit_length * unit_magneticfield)**2 - end if - - sleft = grid_gauss(1) + 0.5d0 * width + dropoff_edge_dist - sright = grid_gauss(gauss_gridpts) - dropoff_edge_dist - 0.5d0 * width - - if (hall_dropoff) then - shift = hallval * tanh(-dpi) / (tanh(-dpi) - tanh(dpi)) - stretch = hallval / (tanh(dpi) - tanh(-dpi)) - - do i = 1, gauss_gridpts - x = grid_gauss(i) - if (sleft - 0.5d0*width <= x .and. x <= sleft + 0.5d0*width) then - hall_field % hallfactor(i) = shift + stretch * tanh(2.0d0 * dpi * (x - sleft) / width) - else if (sleft + 0.5d0*width < x .and. x < sright - 0.5d0*width) then - hall_field % hallfactor(i) = hallval - else if (sright - 0.5d0*width <= x .and. x <= sright + 0.5d0*width) then - hall_field % hallfactor(i) = shift + stretch * tanh(2.0d0 * dpi * (sright - x) / width) - else - hall_field % hallfactor(i) = 0.0d0 - end if - end do - else - hall_field % hallfactor = hallval - end if - - if (inertia_dropoff) then - shift2 = inertiaval * tanh(-dpi) / (tanh(-dpi) - tanh(dpi)) - stretch2 = inertiaval / (tanh(dpi) - tanh(-dpi)) - - do i = 1, gauss_gridpts - x = grid_gauss(i) - if (sleft - 0.5d0*width <= x .and. x <= sleft + 0.5d0*width) then - hall_field % inertiafactor(i) = shift2 + stretch2 * tanh(2.0d0 * dpi * (x - sleft) / width) - else if (sleft + 0.5d0*width < x .and. x < sright - 0.5d0*width) then - hall_field % inertiafactor(i) = inertiaval - else if (sright - 0.5d0*width <= x .and. x <= sright + 0.5d0*width) then - hall_field % inertiafactor(i) = shift2 + stretch2 * tanh(2.0d0 * dpi * (sright - x) / width) - else - hall_field % inertiafactor(i) = 0.0d0 - end if - end do - else - hall_field % inertiafactor = inertiaval - end if - end subroutine set_hallfactor - - !> Retrieves the A-matrix Hall contribution for a given quadblock corresponding - !! to a particular point in the grid and a particular Gaussian weight. - !! If hall_mhd = True, this routine is called n_gauss times - !! for every grid interval. - subroutine get_hallterm(gauss_idx, eps, d_eps_dr, curr_weight, quadblock_Hall, & - h_quadratic, dh_quadratic_dr, h_cubic, dh_cubic_dr) - use mod_global_variables, only: dp_LIMIT, elec_pressure - use mod_equilibrium_params, only: k2, k3 - use mod_equilibrium, only: rho_field, B_field, T_field - use mod_make_subblock, only: subblock, reset_factors, reset_positions - - !> current index in the Gaussian grid - integer, intent(in) :: gauss_idx - !> current value for the scale factor epsilon - real(dp), intent(in) :: eps - !> current value for the derivative of epsilon - real(dp), intent(in) :: d_eps_dr - !> current weight in the Gaussian quadrature - real(dp), intent(in) :: curr_weight - !> the A-quadblock for a particular grid interval - complex(dp), intent(out) :: quadblock_Hall(dim_quadblock, dim_quadblock) - !> array containing the 4 quadratic basis functions - real(dp), intent(in) :: h_quadratic(4) - !> array containing the derivatives of the 4 quadratic basis functions - real(dp), intent(in) :: dh_quadratic_dr(4) - !> array containing the 4 cubic basis functions - real(dp), intent(in) :: h_cubic(4) - !> array containing the derivatives of the 4 cubic basis functions - real(dp), intent(in) :: dh_cubic_dr(4) - - complex(dp), allocatable :: factors(:) - integer, allocatable :: positions(:, :) - - real(dp) :: eps_inv, rho0_inv - real(dp) :: rho0, B01, B02, B03, T0 - real(dp) :: drho0, drB02, dB02_r, dB03, dB02, dT0 - - quadblock_Hall = (0.0d0, 0.0d0) - eps_inv = 1.0d0 / eps - - ! Equilibrium quantities - rho0 = rho_field % rho0(gauss_idx) - drho0 = rho_field % d_rho0_dr(gauss_idx) - rho0_inv = 1.0d0 / rho0 - - T0 = T_field % T0(gauss_idx) - dT0 = T_field % d_T0_dr(gauss_idx) - - ! To be uncommented when B01 is implemented - B01 = 0.0d0 !B_field % B01(gauss_idx) - B02 = B_field % B02(gauss_idx) - dB02 = B_field % d_B02_dr(gauss_idx) - B03 = B_field % B03(gauss_idx) - dB03 = B_field % d_B03_dr(gauss_idx) - - ! Calculate derivatives eps*B02, B02/eps - drB02 = d_eps_dr * B02 + eps * dB02 - dB02_r = eps_inv * dB02 - d_eps_dr * eps_inv**2 * B02 - - ! Setup of matrix elements - - ! Quadratic * Quadratic - call reset_factors(factors, 3) - call reset_positions(positions, 3) - ! A(6, 1) - factors(1) = - eps_inv * rho0_inv**2 * (B02 * drB02 * eps_inv + B03 * dB03) - positions(1, :) = [6, 1] - ! A(6, 5) - factors(2) = 0.0d0 - positions(2, :) = [6, 5] - ! A(6, 6) - factors(3) = - drho0 * rho0_inv**2 * (k2 * B03 * eps_inv - k3 * B02) & - + k3 * rho0_inv * (drB02 * eps_inv - eps * dB02_r) - positions(3, :) = [6, 6] - if (elec_pressure) then - factors(1) = factors(1) - dT0 * eps_inv * rho0_inv - factors(2) = factors(2) + drho0 * eps_inv * rho0_inv - end if - call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, h_quadratic) - - ! d(Quadratic)/dr * Quadratic - call reset_factors(factors, 3) - call reset_positions(positions, 3) - ! A(6, 1) - factors(1) = 0.0d0 - positions(1, :) = [6, 1] - ! A(6, 5) - factors(2) = 0.0d0 - positions(2, :) = [6, 5] - ! A(6, 6) - factors(3) = rho0_inv * (k2 * B03 * eps_inv - k3 * B02) - positions(3, :) = [6, 6] - if (elec_pressure) then - factors(1) = factors(1) - T0 * eps_inv * rho0_inv - factors(2) = factors(2) - eps_inv - end if - call subblock(quadblock_Hall, factors, positions, curr_weight, dh_quadratic_dr, h_quadratic) - - ! Quadratic * Cubic - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(6, 7) - factors(1) = - k3 * eps_inv * rho0_inv * (k2 * B02 * eps_inv + k3 * B03) - positions(1, :) = [6, 7] - ! A(6, 8) - factors(2) = k2 * eps_inv * rho0_inv * (k2 * B02 * eps_inv + k3 * B03) - positions(2, :) = [6, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, h_cubic) - - ! Quadratic * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(6, 7) - factors(1) = drho0 * rho0_inv**2 * B03 * eps_inv - positions(1, :) = [6, 7] - ! A(6, 8) - factors(2) = rho0_inv * (eps * dB02_r - drB02 * eps_inv - drho0 * rho0_inv * B02) - positions(2, :) = [6, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, dh_cubic_dr) - - ! d(Quadratic)/dr * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(6, 7) - factors(1) = - B03 * eps_inv * rho0_inv - positions(1, :) = [6, 7] - ! A(6, 8) - factors(2) = B02 * rho0_inv - positions(2, :) = [6, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, dh_quadratic_dr, dh_cubic_dr) - - ! Cubic * Quadratic - call reset_factors(factors, 6) - call reset_positions(positions, 6) - ! A(7, 1) - factors(1) = 0.0d0 - positions(1, :) = [7, 1] - ! A(7, 5) - factors(2) = 0.0d0 - positions(2, :) = [7, 5] - ! A(7, 6) - factors(3) = - B03 * rho0_inv * ((k2 * eps_inv)**2 + k3**2) - positions(3, :) = [7, 6] - ! A(8, 1) - factors(4) = 0.0d0 - positions(4, :) = [8, 1] - ! A(8, 5) - factors(5) = 0.0d0 - positions(5, :) = [8, 5] - ! A(8, 6) - factors(6) = B02 * rho0_inv * ((k2 * eps_inv)**2 + k3**2) - positions(6, :) = [8, 6] - if (elec_pressure) then - factors(1) = factors(1) + k2 * T0 * rho0_inv * eps_inv**2 - factors(2) = factors(2) + k2 * eps_inv**2 - factors(4) = factors(4) + k3 * T0 * rho0_inv * eps_inv - factors(5) = factors(5) + k3 * eps_inv - end if - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_quadratic) - - ! Cubic * Cubic - call reset_factors(factors, 4) - call reset_positions(positions, 4) - ! A(7, 7) - factors(1) = k3 * drB02 * eps_inv**2 * rho0_inv - positions(1, :) = [7, 7] - ! A(7, 8) - factors(2) = - k2 * drB02 * eps_inv**2 * rho0_inv - positions(2, :) = [7, 8] - ! A(8, 7) - factors(3) = k3 * dB03 * eps_inv * rho0_inv - positions(3, :) = [8, 7] - ! A(8, 8) - factors(4) = - k2 * dB03 * eps_inv * rho0_inv - positions(4, :) = [8, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_cubic) - - ! Cubic * d(Cubic)/dr - call reset_factors(factors, 4) - call reset_positions(positions, 4) - ! A(7, 7) - factors(1) = k2 * B03 * eps_inv**2 * rho0_inv - positions(1, :) = [7, 7] - ! A(7, 8) - factors(2) = k3 * B03 * rho0_inv - positions(2, :) = [7, 8] - ! A(8, 7) - factors(3) = - k2 * B02 * eps_inv**2 * rho0_inv - positions(3, :) = [8, 7] - ! A(8, 8) - factors(4) = - k3 * B02 * rho0_inv - positions(4, :) = [8, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, dh_cubic_dr) - - deallocate(factors) - deallocate(positions) - - if (B01 > dp_LIMIT) then - call get_hallterm_B01(gauss_idx, eps, d_eps_dr, curr_weight, & - quadblock_Hall, h_quadratic, h_cubic, dh_cubic_dr) - end if - end subroutine get_hallterm - - !> Retrieves the A-matrix Hall contribution for a given quadblock corresponding - !! to a particular point in the grid and a particular Gaussian weight. - !! If hall_mhd = True, this routine is called n_gauss times - !! for every grid interval. - subroutine get_hallterm_Bmat(gauss_idx, eps, d_eps_dr, curr_weight, quadblock_HallB, & - h_quadratic, h_cubic, dh_cubic_dr) - use mod_equilibrium_params, only: k2, k3 - use mod_equilibrium, only: rho_field - use mod_make_subblock, only: subblock, reset_factors, reset_positions - - !> current index in the Gaussian grid - integer, intent(in) :: gauss_idx - !> current value for the scale factor epsilon - real(dp), intent(in) :: eps - !> current value for the derivative of epsilon - real(dp), intent(in) :: d_eps_dr - !> current weight in the Gaussian quadrature - real(dp), intent(in) :: curr_weight - !> the A-quadblock for a particular grid interval - complex(dp), intent(out) :: quadblock_HallB(dim_quadblock, dim_quadblock) - !> array containing the 4 quadratic basis functions - real(dp), intent(in) :: h_quadratic(4) - !> array containing the 4 cubic basis functions - real(dp), intent(in) :: h_cubic(4) - !> array containing the derivatives of the 4 cubic basis functions - real(dp), intent(in) :: dh_cubic_dr(4) - - complex(dp), allocatable :: factors(:) - integer, allocatable :: positions(:, :) - - real(dp) :: eps_inv, rho0_inv - real(dp) :: rho0, drho0 - - quadblock_HallB = (0.0d0, 0.0d0) - eps_inv = 1.0d0 / eps - - ! Equilibrium quantities - rho0 = rho_field % rho0(gauss_idx) - drho0 = rho_field % d_rho0_dr(gauss_idx) - rho0_inv = 1.0d0 / rho0 - - ! Setup of matrix elements - - ! Quadratic * Quadratic - call reset_factors(factors, 1) - call reset_positions(positions, 1) - ! A(6, 6) - factors(1) = rho0_inv * ((k2 * eps_inv)**2 + k3**2) - positions(1, :) = [6, 6] - call subblock(quadblock_HallB, factors, positions, curr_weight, h_quadratic, h_quadratic) - - ! Quadratic * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(6, 7) - factors(1) = - k2 * rho0_inv * eps_inv**2 - positions(1, :) = [6, 7] - ! A(6, 8) - factors(2) = - k3 * rho0_inv - positions(2, :) = [6, 8] - call subblock(quadblock_HallB, factors, positions, curr_weight, h_quadratic, dh_cubic_dr) - - ! Cubic * Quadratic - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 6) - factors(1) = drho0 * k2 * eps_inv * rho0_inv**2 - positions(1, :) = [7, 6] - ! A(8, 6) - factors(2) = drho0 * k3 * rho0_inv**2 + d_eps_dr * k3 * eps_inv * rho0_inv - positions(2, :) = [8, 6] - call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, h_quadratic) - - ! d(Cubic)/dr * Quadratic - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 6) - factors(1) = - k2 * eps_inv * rho0_inv - positions(1, :) = [7, 6] - ! A(8, 6) - factors(2) = - k3 * rho0_inv - positions(2, :) = [8, 6] - call subblock(quadblock_HallB, factors, positions, curr_weight, dh_cubic_dr, h_quadratic) - - ! Cubic * Cubic - call reset_factors(factors, 4) - call reset_positions(positions, 4) - ! A(7, 7) - factors(1) = k3**2 * eps_inv * rho0_inv - positions(1, :) = [7, 7] - ! A(7, 8) - factors(2) = - k2 * k3 * eps_inv * rho0_inv - positions(2, :) = [7, 8] - ! A(8, 7) - factors(3) = - k2 * k3 * eps_inv**2 * rho0_inv - positions(3, :) = [8, 7] - ! A(8, 8) - factors(4) = (k2 * eps_inv)**2 * rho0_inv - positions(4, :) = [8, 8] - call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, h_cubic) - - ! Cubic * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 7) - factors(1) = - drho0 * eps_inv * rho0_inv**2 - positions(1, :) = [7, 7] - ! A(8, 8) - factors(2) = - rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv) - positions(2, :) = [8, 8] - call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, dh_cubic_dr) - - ! d(Cubic)/dr * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 7) - factors(1) = eps_inv * rho0_inv - positions(1, :) = [7, 7] - ! A(8, 8) - factors(2) = rho0_inv - positions(2, :) = [8, 8] - call subblock(quadblock_HallB, factors, positions, curr_weight, dh_cubic_dr, dh_cubic_dr) - - deallocate(factors) - deallocate(positions) - end subroutine get_hallterm_Bmat - - !> Creates a quadblock for the A matrix containing the natural boundary - !! conditions coming from the Hall term, depending on the supplied edge. - !! @note Currently, no boundary conditions are applied due to the use of a - !! dropoff profile for the Hall parameter, which is zero at the - !! edges. @endnote - subroutine hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge) - use mod_global_variables, only: gauss_gridpts, dim_subblock, elec_pressure - use mod_logging, only: log_message - use mod_equilibrium, only: rho_field, B_field, T_field - use mod_grid, only: eps_grid - use mod_equilibrium_params, only: k2, k3 - - !> the quadblock corresponding to the left/right edge - complex(dp), intent(out) :: quadblock_Hall(dim_quadblock, dim_quadblock) - !> the edge, either "l_edge" or "r_edge" - character(len=6), intent(in) :: edge - !> boolean indicating if perpendicular thermal conduction is included or not - logical, intent(in) :: kappa_perp_is_zero - - complex(dp), allocatable :: surface_terms(:) - real(dp) :: eps, eps_inv, rho0_inv - real(dp) :: rho0, B02, B03, T0 - integer, allocatable :: positions(:, :) - integer :: idx, i - - quadblock_Hall = (0.0d0, 0.0d0) - - if (edge == 'l_edge') then - idx = 1 - else if (edge == 'r_edge') then - idx = gauss_gridpts - else - call log_message('Hall boundaries: wrong edge supplied' // edge, level='error') - end if - - ! retrieve variables at current edge - eps = eps_grid(idx) - eps_inv = 1.0d0 / eps - rho0 = rho_field % rho0(idx) - rho0_inv = 1.0d0 / rho0 - T0 = T_field % T0(idx) - B02 = B_field % B02(idx) - B03 = B_field % B03(idx) - - allocate(positions(4, 2)) - allocate(surface_terms(4)) - - ! surface term for element (6, 5) - surface_terms(1) = 0.0d0 - positions(1, :) = [6, 5] - ! surface term for element (6, 6) - surface_terms(2) = - rho0_inv * (k2 * B03 * eps_inv - k3 * B02) - positions(2, :) = [6, 6] - ! surface term for element (6, 7) - surface_terms(3) = B03 * eps_inv * rho0_inv - positions(3, :) = [6, 7] - ! surface term for element (6, 8) - surface_terms(4) = - B02 * rho0_inv - positions(4, :) = [6, 8] - - ! T1 is zero at the wall if perpendicular thermal conduction is included - if (elec_pressure .and. kappa_perp_is_zero) then - surface_terms(2) = surface_terms(2) + eps_inv - end if - - ! l_edge: add to bottom-right of 2x2 block, for top-left subblock only - ! r_edge: add to bottom-right of 2x2 block, for bottom-right subblock only - if (edge == 'l_edge') then - positions = 2 * positions - else if (edge == 'r_edge') then - positions = 2 * positions + dim_subblock - end if - - do i = 1, 2 - quadblock_Hall(positions(i, 1), positions(i, 2)) = quadblock_Hall(positions(i, 1), positions(i, 2)) + surface_terms(i) - end do - do i = 3, size(surface_terms) - quadblock_Hall(positions(i, 1), positions(i, 2)) = quadblock_Hall(positions(i, 1), positions(i, 2)) + surface_terms(i) - end do - - deallocate(positions) - deallocate(surface_terms) - end subroutine hall_boundaries - - !> Retrieves the B01 (constant) contributions to the A-matrix Hall - !! contribution for a given quadblock corresponding - !! to a particular point in the grid and a particular Gaussian weight. - !! If hall_mhd = True, this routine is called n_gauss times - !! for every grid interval. - !! @note Usage of a constant B01 component is currently not implemented. This - !! subroutine supports a future extension of the code. @endnote - subroutine get_hallterm_B01(gauss_idx, eps, d_eps_dr, curr_weight, & - quadblock_Hall, h_quadratic, h_cubic, dh_cubic_dr) - use mod_global_variables, only: ic - use mod_equilibrium_params, only: k2, k3 - use mod_equilibrium, only: rho_field, B_field - use mod_make_subblock, only: subblock, reset_factors, reset_positions - - !> current index in the Gaussian grid - integer, intent(in) :: gauss_idx - !> current value for the scale factor epsilon - real(dp), intent(in) :: eps - !> current value for the derivative of epsilon - real(dp), intent(in) :: d_eps_dr - !> current weight in the Gaussian quadrature - real(dp), intent(in) :: curr_weight - !> the A-quadblock for a particular grid interval - complex(dp), intent(inout):: quadblock_Hall(dim_quadblock, dim_quadblock) - !> array containing the 4 quadratic basis functions - real(dp), intent(in) :: h_quadratic(4) - !> array containing the 4 cubic basis functions - real(dp), intent(in) :: h_cubic(4) - !> array containing the derivatives of the 4 cubic basis functions - real(dp), intent(in) :: dh_cubic_dr(4) - - complex(dp), allocatable :: factors(:) - integer, allocatable :: positions(:, :) - - real(dp) :: eps_inv, rho0_inv - real(dp) :: rho0, B01, B02, B03 - real(dp) :: drho0, drB02, dB02_r, dB03, dB02 - - eps_inv = 1.0d0 / eps - - ! Equilibrium quantities - rho0 = rho_field % rho0(gauss_idx) - drho0 = rho_field % d_rho0_dr(gauss_idx) - rho0_inv = 1.0d0 / rho0 - - ! To be uncommented when B01 is implemented - B01 = 0.0d0 !B_field % B01(gauss_idx) - B02 = B_field % B02(gauss_idx) - dB02 = B_field % d_B02_dr(gauss_idx) - B03 = B_field % B03(gauss_idx) - dB03 = B_field % d_B03_dr(gauss_idx) - - ! Calculate derivatives eps*B02, B02/eps - drB02 = d_eps_dr * B02 + eps * dB02 - dB02_r = eps_inv * dB02 - d_eps_dr * eps_inv**2 * B02 - - ! Setup of matrix elements - - ! Cubic * Quadratic - call reset_factors(factors, 4) - call reset_positions(positions, 4) - ! A(7, 1) - factors(1) = - ic * eps_inv**2 * rho0_inv**2 * B01 * drB02 - positions(1, :) = [7, 1] - ! A(7, 6) - factors(2) = ic * k3 * B01 * rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv) - positions(2, :) = [7, 6] - ! A(8, 1) - factors(3) = - ic * eps_inv * rho0_inv**2 * B01 * dB03 - positions(3, :) = [8, 1] - ! A(8, 6) - factors(4) = - ic * drho0 * rho0_inv**2 * k2 * B01 * eps_inv - positions(4, :) = [8, 6] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_quadratic) - - ! d(Cubic)/dr * Quadratic - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 6) - factors(1) = - ic * k3 * B01 * rho0_inv - positions(1, :) = [7, 6] - ! A(8, 6) - factors(2) = ic * k2 * B01 * eps_inv * rho0_inv - positions(2, :) = [8, 6] - call subblock(quadblock_Hall, factors, positions, curr_weight, dh_cubic_dr, h_quadratic) - - ! Cubic * Cubic - call reset_factors(factors, 4) - call reset_positions(positions, 4) - ! A(7, 7) - factors(1) = - ic * k2 * k3 * B01 * eps_inv**2 * rho0_inv - positions(1, :) = [7, 7] - ! A(7, 8) - factors(2) = ic * k2**2 * B01 * eps_inv**2 * rho0_inv - positions(2, :) = [7, 8] - ! A(8, 7) - factors(3) = - ic * k3**2 * B01 * eps_inv * rho0_inv - positions(3, :) = [8, 7] - ! A(8, 8) - factors(4) = ic * k2 * k3 * B01 * eps_inv * rho0_inv - positions(4, :) = [8, 8] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_cubic) - - ! Cubic * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 8) - factors(1) = ic * B01 * rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv) - positions(1, :) = [7, 8] - ! A(8, 7) - factors(2) = ic * B01 * eps_inv * drho0 * rho0_inv**2 - positions(2, :) = [8, 7] - call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, dh_cubic_dr) - - ! d(Cubic)/dr * d(Cubic)/dr - call reset_factors(factors, 2) - call reset_positions(positions, 2) - ! A(7, 8) - factors(1) = - ic * B01 * rho0_inv - positions(1, :) = [7, 8] - ! A(8, 7) - factors(2) = - ic * B01 * eps_inv * rho0_inv - positions(2, :) = [8, 7] - call subblock(quadblock_Hall, factors, positions, curr_weight, dh_cubic_dr, dh_cubic_dr) - - deallocate(factors) - deallocate(positions) - end subroutine get_hallterm_B01 - -end module mod_hallmhd From b1bf3cef1e62984abfc84d2e378c621deb7bc6ea Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Thu, 16 Sep 2021 11:40:41 +0200 Subject: [PATCH 2/9] pylbo: cleanup continua --- .../pylbo/visualisation/continua.py | 2 + .../pylbo/visualisation/spectra.py | 63 +++++-------------- 2 files changed, 18 insertions(+), 47 deletions(-) diff --git a/post_processing/pylbo/visualisation/continua.py b/post_processing/pylbo/visualisation/continua.py index 02c01121..4c1048e3 100644 --- a/post_processing/pylbo/visualisation/continua.py +++ b/post_processing/pylbo/visualisation/continua.py @@ -162,6 +162,8 @@ def __init__(self, interactive): super().__init__(interactive) self.continua_names = CONTINUA_NAMES self._continua_colors = CONTINUA_COLORS + self.marker = "." + self.markersize = 6 @property def continua_colors(self): diff --git a/post_processing/pylbo/visualisation/spectra.py b/post_processing/pylbo/visualisation/spectra.py index 180ba47d..1bde09b4 100644 --- a/post_processing/pylbo/visualisation/spectra.py +++ b/post_processing/pylbo/visualisation/spectra.py @@ -187,45 +187,20 @@ def add_continua(self, interactive=True): continuum = self.dataset.continua[key] if self._c_handler.check_if_all_zero(continuum=continuum): continue - min_value = np.min(continuum) - max_value = np.max(continuum) - # check if continua are complex - if np.any(np.iscomplex(continuum)): - item = self.ax.scatter( - continuum.real * self.x_scaling, - continuum.imag * self.y_scaling, - marker=".", - color=color, - linewidth=self.markersize, - alpha=self.alpha / 1.5, - label=key, - ) - elif self._c_handler.check_if_collapsed(continuum=continuum): - item = self.ax.scatter( - min_value * self.x_scaling, - 0, - marker=self._c_handler.marker, - s=self._c_handler.markersize, - c=color, - alpha=self._c_handler.alpha_point, - label=key, - ) - else: - props = dict( - facecolor=colors.to_rgba(color, self._c_handler.alpha_region), - edgecolor=colors.to_rgba(color, self._c_handler.alpha_point), - label=key, - ) - if key == "thermal": - item = self.ax.axhspan( - min_value.imag * self.y_scaling, - max_value * self.y_scaling, - **props, - ) - else: - item = self.ax.axvspan( - min_value * self.x_scaling, max_value * self.x_scaling, **props - ) + realvals = continuum.real * self.x_scaling + imagvals = continuum.imag * self.y_scaling + if self._c_handler.check_if_collapsed(continuum=continuum): + realvals = np.min(realvals) + imagvals = 0 + item = self.ax.scatter( + realvals, + imagvals, + marker=self._c_handler.marker, + linewidth=self._c_handler.markersize, + c=color, + alpha=self._c_handler.alpha_point, + label=key, + ) self._c_handler.add(item) self._c_handler.legend = self.ax.legend(**self._c_handler.legend_properties) super().add_continua(self._c_handler.interactive) @@ -251,11 +226,6 @@ def add_derived_eigenfunctions(self): ) super()._enable_interface(handle=self._def_handler) - def _ensure_smooth_continuum(self, continuum): - # TODO: this method should split the continuum into multiple parts, such that - # regions that lie far apart are not connected by a line. - pass - class MultiSpectrumPlot(SpectrumFigure): """ @@ -432,9 +402,8 @@ def add_continua(self, interactive=True): self._c_handler.continua_names, self._c_handler.continua_colors ): # we skip duplicates if eigenvalues are squared - if self.use_squared_omega: - if key in ("slow-", "alfven-"): - continue + if self.use_squared_omega and key in ("slow-", "alfven-"): + continue # retrieve continuum, calculate region boundaries continuum = self.dataseries.continua[key] ** self._w_pow if self.use_real_parts: From 1bdc60aa04638e8aebf807c2f17b68d45c328175 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Mon, 27 Sep 2021 12:09:09 +0200 Subject: [PATCH 3/9] updated parfile template --- .gitignore | 1 + legolas_config.par | 68 +++++++++++++++++++++++----------------------- 2 files changed, 35 insertions(+), 34 deletions(-) diff --git a/.gitignore b/.gitignore index b6662d91..c4be27f1 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ output/ build/ legolas test_legolas +legolas_config.par *.DS_Store *.aux diff --git a/legolas_config.par b/legolas_config.par index bd099444..8a74e3e2 100644 --- a/legolas_config.par +++ b/legolas_config.par @@ -1,58 +1,58 @@ &gridlist - ! geometry = 'Cartesian' - ! geometry = 'cylindrical' - ! x_start = 0.0d0 - ! x_end = 1.0d0 - gridpoints = 51 - mesh_accumulation = .false. + geometry = "cylindrical" + x_start = 0.0d0 + x_end = 1.0d0 + gridpoints = 50 / &solvelist solver = "QR-invert" ! solver = "arnoldi" ! arpack_mode = "shift-invert" - ! number_of_eigenvalues = 10 + ! number_of_eigenvalues = 50 ! which_eigenvalues = "LI" - ! sigma = (0.0d0, 0.01d0) + ! sigma = (-0.13d0, 0.0d0) / &physicslist - ! hall_mhd = .true. + ! flow = .true. / &equilibriumlist - ! equilibrium_type = "adiabatic_homo" ! Cartesian - ! equilibrium_type = "resistive_homo" ! Cartesian - ! equilibrium_type = "gravito_acoustic" ! Cartesian - ! equilibrium_type = "gravito_mhd" ! Cartesian - equilibrium_type = "resistive_tearing" ! Cartesian - ! equilibrium_type = "resistive_tearing_flow" ! Cartesian - ! equilibrium_type = "suydam_cluster" ! cylindrical - ! equilibrium_type = "rotating_plasma_cylinder" ! cylindrical - ! equilibrium_type = "kelvin_helmholtz_cd" ! cylindrical - ! equilibrium_type = "internal_kink" ! cylindrical - ! equilibrium_type = "RTI_theta_pinch" ! cylindrical - ! equilibrium_type = "discrete_alfven" ! cylindrical - ! equilibrium_type = "MRI_accretion" ! cylindrical - ! equilibrium_type = "interchange_modes" ! Cartesian - ! equilibrium_type = "constant_current_tokamak" ! cylindrical - ! equilibrium_type = "magnetothermal_instabilities" - ! equilibrium_type = "photospheric_flux_tube" - ! equilibrium_type = "coronal_flux_tube" - ! equilibrium_type = "resonant_absorption" + ! equilibrium_type = "adiabatic_homo" + ! equilibrium_type = "constant_current_tokamak" + ! equilibrium_type = "couette_flow" + ! equilibrium_type = "discrete_alfven" ! equilibrium_type = "gold_hoyle" - ! equilibrium_type = "rayleigh_taylor" + ! equilibrium_type = "gravito_acoustic" + ! equilibrium_type = "gravito_mhd" + ! equilibrium_type = "coronal_flux_tube" + ! equilibrium_type = "photospheric_flux_tube" + ! equilibrium_type = "interchange_modes" + ! equilibrium_type = "internal_kink" ! equilibrium_type = "kelvin_helmholtz" + ! equilibrium_type = "kelvin_helmholtz_cd" + ! equilibrium_type = "MRI_accretion" + ! equilibrium_type = "magnetothermal_instabilities" + ! equilibrium_type = "rayleigh_taylor" ! equilibrium_type = "RTI_KHI" + ! equilibrium_type = "RTI_theta_pinch" + ! equilibrium_type = "resistive_homo" + ! equilibrium_type = "resistive_tearing" + ! equilibrium_type = "resistive_tearing_flow" + ! equilibrium_type = "resonant_absorption" + ! equilibrium_type = "rotating_plasma_cylinder" ! equilibrium_type = "isothermal_atmosphere" - ! equilibrium_type = "harris_sheet" - ! equilibrium_type = "user_defined" + equilibrium_type = "suydam_cluster" - ! boundary_type = "wall" + ! equilibrium_type = "user_defined" / &savelist - write_matrices = .false. - write_eigenfunctions = .true. show_results = .true. + write_eigenfunctions = .true. + write_derived_eigenfunctions = .true. + ! write_eigenfunction_subset = .true. + ! eigenfunction_subset_center = (-0.13d0, 0.0d0) + ! eigenfunction_subset_radius = 0.02d0 / From 1fa18f0da6d5e60c8f84d828ba1ccc66c06ac520 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Wed, 29 Sep 2021 13:40:15 +0200 Subject: [PATCH 4/9] updated PR and issue templates --- .github/issue_template.md | 27 +++++++++++++++++++ .github/pull_request_template.md | 45 ++++++++------------------------ 2 files changed, 38 insertions(+), 34 deletions(-) create mode 100644 .github/issue_template.md diff --git a/.github/issue_template.md b/.github/issue_template.md new file mode 100644 index 00000000..13916e48 --- /dev/null +++ b/.github/issue_template.md @@ -0,0 +1,27 @@ +## Issue description + + +## Feature request + + + +## Bug report +**Minimal example for reproduction** + +``` + +``` + +**Actual result** + + +**Expected result** + + +**Version info** +- Legolas version: +- Pylbo version: +- Python version: +- Version of relevant libraries (matplotlib, numpy, etc.): \ No newline at end of file diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 8205178b..91922b56 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -1,42 +1,19 @@ -Hi, thanks for the pull request! To make sure everything stays nice and consistent, here is a default template which you can use. -It also includes a checklist of small things that are usually forgotten or overlooked. You may remove parts that are not applicable to the current PR -(including these sentences). - -Cheers, - -Legolas dev team - ## PR description - -> fill or delete + ## New features +**Legolas** +- -> fill or delete +**Pylbo** +- ## Bugfixes +**Legolas** +- -> fill or delete - +**Pylbo** +- ## Checklist - -**Source code stuff** -- [ ] everything is nicely formatted using [Black-like code style](https://black.readthedocs.io/en/stable/the_black_code_style.html) -- [ ] use-statements should use `use mod_xxx, only:` as much as possible -- [ ] things that can go in a separate module, should go in a separate module -- [ ] code additions|changes are also added|changed in the docstrings - -**Testing stuff** -- [ ] all tests pass -- [ ] in case of new features, new tests are added - -**Website stuff** -- [ ] relevant pages on the website are edited, if relevant -- [ ] on the edited pages, the `last_modified_at` frontmatter key is updated -- [ ] all links render as they should, check through a local `bundle exec jekyll serve` in the `docs` directory -- [ ] additions to the parfile are also added to `parameters_file.md` - -**Release stuff** -(only if this is a merge-in-master release) -- [ ] bump version number -- [ ] changelog pages on the website are updated +- [ ] all tests pass (and for new features, new tests are added) +- [ ] documentation has been updated (if applicable) From 067dab553d0bcb5917ef2c5f684bc271ad88e049 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Tue, 5 Oct 2021 17:15:46 +0200 Subject: [PATCH 5/9] pylbo: bugfix for derived eigenfunctions on a multispectrum plot --- post_processing/pylbo/visualisation/spectra.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/post_processing/pylbo/visualisation/spectra.py b/post_processing/pylbo/visualisation/spectra.py index 1bde09b4..e65ab48c 100644 --- a/post_processing/pylbo/visualisation/spectra.py +++ b/post_processing/pylbo/visualisation/spectra.py @@ -457,7 +457,7 @@ def add_derived_eigenfunctions(self): self._def_handler = DerivedEigenfunctionHandler( self.dataseries, self._def_ax, self.ax ) - super()._enable_interface(handle=self._pp_handler) + super()._enable_interface(handle=self._def_handler) class MergedSpectrumPlot(SpectrumFigure): From bd7816a12d26b3bb3f7563fc15164f623bbc27e9 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Tue, 5 Oct 2021 17:16:19 +0200 Subject: [PATCH 6/9] pylbo: fixes broken colourmapping for interactive legends on multispectrum plots --- post_processing/pylbo/visualisation/legend_interface.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/post_processing/pylbo/visualisation/legend_interface.py b/post_processing/pylbo/visualisation/legend_interface.py index f22053c2..fc8bac75 100644 --- a/post_processing/pylbo/visualisation/legend_interface.py +++ b/post_processing/pylbo/visualisation/legend_interface.py @@ -107,6 +107,9 @@ def make_legend_pickable(self): add_pickradius_to_item(item=legend_item, pickradius=self.pickradius) # add an attribute to this artist to tell it's from a legend setattr(legend_item, "is_legend_item", True) + # make sure colourmapping is done properly + legend_item.set_facecolor(drawn_item.get_facecolor()) + legend_item.set_edgecolor(drawn_item.get_edgecolor()) # we make the regions invisible until clicked, or set visible as default if self._make_visible_by_default: legend_item.set_alpha(self.alpha_point) From 71864c7bc3fe67cd0702842571b8669f372e7884 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Tue, 5 Oct 2021 17:16:47 +0200 Subject: [PATCH 7/9] minor tweaks to internal kink equilibrium --- .../smod_equil_internal_kink_instability.f08 | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/equilibria/smod_equil_internal_kink_instability.f08 b/src/equilibria/smod_equil_internal_kink_instability.f08 index f81c0018..688ea44e 100644 --- a/src/equilibria/smod_equil_internal_kink_instability.f08 +++ b/src/equilibria/smod_equil_internal_kink_instability.f08 @@ -24,7 +24,7 @@ module subroutine internal_kink_eq() use mod_equilibrium_params, only: cte_rho0, cte_v03, cte_p0, alpha - real(dp) :: r, x + real(dp) :: r, x, a0 real(dp) :: J0, J1, J2, DJ0, DJ1 integer :: i @@ -33,12 +33,13 @@ module subroutine internal_kink_eq() ) call initialise_grid() + a0 = x_end if (use_defaults) then ! LCOV_EXCL_START flow = .true. cte_rho0 = 1.0d0 cte_v03 = 1.0d0 cte_p0 = 9.0d0 - alpha = 5.0d0 / x_end + alpha = 5.0d0 / a0 k2 = 1.0d0 k3 = 0.16d0 * alpha @@ -46,7 +47,7 @@ module subroutine internal_kink_eq() do i = 1, gauss_gridpts r = grid_gauss(i) - x = r / x_end + x = r / a0 J0 = bessel_jn(0, alpha * x) J1 = bessel_jn(1, alpha * x) @@ -54,15 +55,17 @@ module subroutine internal_kink_eq() DJ0 = -alpha * J1 DJ1 = alpha * (0.5d0 * J0 - 0.5d0 * J2) - rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2) - v_field % v03(i) = cte_v03 * (1.0d0 - x**2) + rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2 / a0**2) + v_field % v03(i) = cte_v03 * (1.0d0 - x**2 / a0**2) B_field % B02(i) = J1 B_field % B03(i) = J0 B_field % B0(i) = sqrt((B_field % B02(i))**2 + (B_field % B03(i))**2) T_field % T0(i) = cte_p0 / (rho_field % rho0(i)) - rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x - T_field % d_T0_dr(i) = 2.0d0 * x * cte_p0 / (cte_rho0 * (1.0d0 - x**2)**2) + rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x / a0 + T_field % d_T0_dr(i) = 2.0d0 * x * cte_p0 / ( & + a0**2 * cte_rho0 * (1.0d0 - x**2)**2 & + ) v_field % d_v03_dr(i) = -2.0d0 * cte_v03 * x B_field % d_B02_dr(i) = DJ1 B_field % d_B03_dr(i) = DJ0 From 794bf3dfd048ce92c31d2615a8a774a731b277de Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Wed, 6 Oct 2021 10:57:16 +0200 Subject: [PATCH 8/9] pylbo: fix AttributeError for non-continua legend items --- .../pylbo/visualisation/legend_interface.py | 8 +++++--- tests/pylbo_tests/test_visualisations_spectra.py | 10 ++++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/post_processing/pylbo/visualisation/legend_interface.py b/post_processing/pylbo/visualisation/legend_interface.py index fc8bac75..c831aafd 100644 --- a/post_processing/pylbo/visualisation/legend_interface.py +++ b/post_processing/pylbo/visualisation/legend_interface.py @@ -1,5 +1,6 @@ import matplotlib.collections as mpl_collections import matplotlib.lines as mpl_lines +import matplotlib.patches as mpl_patches import numpy as np from pylbo.utilities.toolbox import add_pickradius_to_item @@ -107,9 +108,10 @@ def make_legend_pickable(self): add_pickradius_to_item(item=legend_item, pickradius=self.pickradius) # add an attribute to this artist to tell it's from a legend setattr(legend_item, "is_legend_item", True) - # make sure colourmapping is done properly - legend_item.set_facecolor(drawn_item.get_facecolor()) - legend_item.set_edgecolor(drawn_item.get_edgecolor()) + # make sure colourmapping is done properly, only for continua regions + if isinstance(legend_item, mpl_patches.Rectangle): + legend_item.set_facecolor(drawn_item.get_facecolor()) + legend_item.set_edgecolor(drawn_item.get_edgecolor()) # we make the regions invisible until clicked, or set visible as default if self._make_visible_by_default: legend_item.set_alpha(self.alpha_point) diff --git a/tests/pylbo_tests/test_visualisations_spectra.py b/tests/pylbo_tests/test_visualisations_spectra.py index 20b5990a..d9d29c64 100644 --- a/tests/pylbo_tests/test_visualisations_spectra.py +++ b/tests/pylbo_tests/test_visualisations_spectra.py @@ -33,6 +33,16 @@ def test_spectrum_plot_eigenfunctions(ds_v112): assert getattr(p, "ef_ax", None) is not None +def test_spectrum_plot_continua_and_eigenfunctions(ds_v112): + p = pylbo.plot_spectrum(ds_v112) + p.add_continua() + p.add_eigenfunctions() + labels = p.ax.get_legend_handles_labels()[1] + assert len(labels) > 2 + for attr in ("c_handler", "ef_handler", "ef_ax"): + assert getattr(p, attr, None) is not None + + def test_spectrum_plot_eigenfunctions_notpresent(ds_v100): p = pylbo.plot_spectrum(ds_v100) with pytest.raises(EigenfunctionsNotPresent): From 5d8b3a3be42acf5d72da5ba7bff67e01c7725e1e Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Wed, 6 Oct 2021 11:07:51 +0200 Subject: [PATCH 9/9] remove viscosity warning --- src/main.f08 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/main.f08 b/src/main.f08 index 96c33a46..e11287ae 100644 --- a/src/main.f08 +++ b/src/main.f08 @@ -62,7 +62,7 @@ subroutine initialisation() use mod_input, only: read_parfile, get_parfile use mod_equilibrium, only: initialise_equilibrium, set_equilibrium use mod_logging, only: print_logo - use mod_global_variables, only: viscosity, hall_mhd + use mod_global_variables, only: hall_mhd character(len=str_len) :: parfile integer :: nb_evs @@ -89,14 +89,10 @@ subroutine initialisation() call set_equilibrium() ! TODO: remove this warning when fully tested - if (viscosity) then - call log_message( & - "using viscous MHD, note that this is not yet fully tested!", level="warning" & - ) - end if if (hall_mhd) then call log_message( & - "using Hall MHD, note that this does not yet work properly!", level="warning" & + "using Hall MHD, note that some features are not yet **fully** tested!", & + level="warning" & ) end if