From 95f08794d4f0ac4aab9dc901ac1c072650a5968e Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 17 Apr 2025 19:20:38 +0200 Subject: [PATCH 1/6] remove .gitkeep in conflict with symlink --- run/.gitkeep | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 run/.gitkeep diff --git a/run/.gitkeep b/run/.gitkeep deleted file mode 100644 index e69de29..0000000 From e6169497c97875388e87477f365b2f46ff8e1a82 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 17 Apr 2025 21:14:27 +0200 Subject: [PATCH 2/6] use explicit temporary arrays --- src/mephit_iter.F90 | 20 +- src/mephit_mesh.F90 | 873 ++++++++++++++++++++++++++++---------------- 2 files changed, 563 insertions(+), 330 deletions(-) diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index 7ad4e19..2e94b2d 100644 --- a/src/mephit_iter.F90 +++ b/src/mephit_iter.F90 @@ -1285,7 +1285,8 @@ subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_ complex(dp), intent(out) :: Phi_mn(0:, mesh%m_res_min:) complex(dp), intent(out) :: Phi_aligned_mn(0:, mesh%m_res_min:) integer :: m, m_res, kf_min - complex(dp) :: Bmnpsi_over_B0phi(0:mesh%nflux) + real(dp) :: Bmnrho_over_B0theta(1:mesh%nflux), rtemp(1:mesh%nflux) + complex(dp) :: Bmnpsi_over_B0phi(1:mesh%nflux) type(vec_polmodes_t) :: Bmn kf_min = 0 ! max(1, mesh%res_ind(mesh%m_res_min) / 4) @@ -1297,13 +1298,18 @@ subroutine compute_flr2_current(Bn, jmnpar_over_Bmod, flr2, Phi_mn, Phi_aligned_ call RT0_poloidal_modes(Bn, Bmn) do m = mesh%m_res_min, mesh%m_res_max m_res = -equil%cocos%sgn_q * m - call resample1d(fs_half%psi, Bmn%coeff_rad(m_res, :)%Re / fs_half%q * & - equil%cocos%sgn_dpsi, fs%psi(1:), Bmnpsi_over_B0phi(1:)%Re, 3) - call resample1d(fs_half%psi, Bmn%coeff_rad(m_res, :)%Im / fs_half%q * & - equil%cocos%sgn_dpsi, fs%psi(1:), Bmnpsi_over_B0phi(1:)%Im, 3) + Bmnrho_over_B0theta(:) = Bmn%coeff_rad(m_res, :)%Re + call resample1d(fs_half%psi, Bmnrho_over_B0theta / fs_half%q * & + equil%cocos%sgn_dpsi, fs%psi(1:), rtemp, 3) + Bmnpsi_over_B0phi%Re = rtemp + Bmnrho_over_B0theta(:) = Bmn%coeff_rad(m_res, :)%Im + call resample1d(fs_half%psi, Bmnrho_over_B0theta / fs_half%q * & + equil%cocos%sgn_dpsi, fs%psi(1:), rtemp, 3) + Bmnpsi_over_B0phi%Im = rtemp call flr2_response_current(flr2, 1, m_res, mesh%nflux, & - Bmnpsi_over_B0phi(1:), jmnpar_over_Bmod%coeff(m_res, 1:), Phi_mn(1:, m)) - Phi_aligned_mn(:, m) = imun * Bmnpsi_over_B0phi * fs%q * dPhi0_dpsi%y / (m_res + mesh%n * fs%q) + Bmnpsi_over_B0phi, jmnpar_over_Bmod%coeff(m_res, 1:), Phi_mn(1:, m)) + Phi_aligned_mn(1:, m) = imun * Bmnpsi_over_B0phi * & + fs%q(1:) * dPhi0_dpsi%y(1:) / (m_res + mesh%n * fs%q(1:)) jmnpar_over_Bmod%coeff(m_res, :kf_min) = (0d0, 0d0) ! suppress spurious current near axis end do call vec_polmodes_deinit(Bmn) diff --git a/src/mephit_mesh.F90 b/src/mephit_mesh.F90 index e958993..53738ea 100644 --- a/src/mephit_mesh.F90 +++ b/src/mephit_mesh.F90 @@ -299,11 +299,31 @@ end subroutine rtree_query real(dp) :: R, Z, psi, theta, sqrt_g, B0_R, B0_phi, B0_Z, dR_dtheta, dZ_dtheta end type coord_cache_t + interface coord_cache_write + module procedure coord_cache_write_1 + module procedure coord_cache_write_2 + end interface coord_cache_write + + interface coord_cache_read + module procedure coord_cache_read_1 + module procedure coord_cache_read_2 + end interface coord_cache_read + type :: field_cache_t real(dp) :: psi, theta, B0(3), j0(3), Bmod, dBmod_dR, dBmod_dZ, & dB0_dR(3), dB0_dZ(3), dj0_dR(3), dj0_dZ(3) end type field_cache_t + interface field_cache_write + module procedure field_cache_write_1 + module procedure field_cache_write_2 + end interface field_cache_write + + interface field_cache_read + module procedure field_cache_read_1 + module procedure field_cache_read_2 + end interface field_cache_read + type :: shielding_t real(dp), allocatable :: GL_weights(:) type(coord_cache_t), allocatable :: sample_Ires(:, :) @@ -454,354 +474,561 @@ subroutine flux_func_cache_read(cache, file, dataset) call h5_close(h5id_root) end subroutine flux_func_cache_read - subroutine coord_cache_write(cache, file, dataset, comment) + subroutine coord_cache_write_1(cache, file, group, comment) use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close - class(coord_cache_t), dimension(..), intent(in) :: cache + type(coord_cache_t), dimension(:), intent(in) :: cache character(len = *), intent(in) :: file - character(len = *), intent(in) :: dataset + character(len = *), intent(in) :: group character(len = *), intent(in) :: comment + character(len = len_trim(group)) :: grp + character(len = len_trim(comment)) :: cmnt integer(HID_T) :: h5id_root + integer, dimension(:), allocatable :: itemp + real(dp), dimension(:), allocatable :: rtemp - select rank (cache) - rank (2) - call h5_open_rw(file, h5id_root) - call h5_create_parent_groups(h5id_root, trim(adjustl(dataset)) // '/') - call h5_add(h5id_root, trim(adjustl(dataset)) // '/ktri', cache(:, :)%ktri, & - lbound(cache), ubound(cache), & - comment = 'triangle index of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/R', cache(:, :)%R, & - lbound(cache), ubound(cache), unit = 'cm', & - comment = 'R coordinate of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/Z', cache(:, :)%Z, & - lbound(cache), ubound(cache), unit = 'cm', & - comment = 'Z coordinate of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:, :)%psi, & - lbound(cache), ubound(cache), unit = 'Mx', & - comment = 'poloidal flux at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:, :)%theta, & - lbound(cache), ubound(cache), unit = 'rad', & - comment = 'flux poloidal angle at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/sqrt_g', cache(:, :)%sqrt_g, & - lbound(cache), ubound(cache), unit = 'cm G^-1', & - comment = 'Jacobian at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:, :)%B0_R, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:, :)%B0_phi, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:, :)%B0_Z, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dR_dtheta', cache(:, :)%dR_dtheta, & - lbound(cache), ubound(cache), unit = 'cm rad^-1', & - comment = 'Jacobian element (R, theta) at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dZ_dtheta', cache(:, :)%dZ_dtheta, & - lbound(cache), ubound(cache), unit = 'cm rad^-1', & - comment = 'Jacobian element (Z, theta) at ' // trim(adjustl(comment))) - call h5_close(h5id_root) - rank (1) - call h5_open_rw(file, h5id_root) - call h5_create_parent_groups(h5id_root, trim(adjustl(dataset)) // '/') - call h5_add(h5id_root, trim(adjustl(dataset)) // '/ktri', cache(:)%ktri, & - lbound(cache), ubound(cache), & - comment = 'triangle index of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/R', cache(:)%R, & - lbound(cache), ubound(cache), unit = 'cm', & - comment = 'R coordinate of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/Z', cache(:)%Z, & - lbound(cache), ubound(cache), unit = 'cm', & - comment = 'Z coordinate of ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:)%psi, & - lbound(cache), ubound(cache), unit = 'Mx', & - comment = 'poloidal flux at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:)%theta, & - lbound(cache), ubound(cache), unit = 'rad', & - comment = 'flux poloidal angle at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/sqrt_g', cache(:)%sqrt_g, & - lbound(cache), ubound(cache), unit = 'cm G^-1', & - comment = 'Jacobian at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:)%B0_R, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:)%B0_phi, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:)%B0_Z, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dR_dtheta', cache(:)%dR_dtheta, & - lbound(cache), ubound(cache), unit = 'cm rad^-1', & - comment = 'Jacobian element (R, theta) at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dZ_dtheta', cache(:)%dZ_dtheta, & - lbound(cache), ubound(cache), unit = 'cm rad^-1', & - comment = 'Jacobian element (Z, theta) at ' // trim(adjustl(comment))) - call h5_close(h5id_root) - rank default - error stop 'coord_cache_write: rank-1 or rank-2 array expected for argument ''cache''' - end select - end subroutine coord_cache_write + grp = trim(group) + cmnt = trim(comment) + allocate(itemp(size(cache)), rtemp(size(cache))) + call h5_open_rw(file, h5id_root) + call h5_create_parent_groups(h5id_root, grp // '/') + itemp(:) = cache%ktri + call h5_add(h5id_root, grp // '/ktri', itemp, & + lbound(cache), ubound(cache), & + comment = 'triangle index of ' // cmnt) + rtemp(:) = cache%R + call h5_add(h5id_root, grp // '/R', rtemp, & + lbound(cache), ubound(cache), unit = 'cm', & + comment = 'R coordinate of ' // cmnt) + rtemp(:) = cache%Z + call h5_add(h5id_root, grp // '/Z', rtemp, & + lbound(cache), ubound(cache), unit = 'cm', & + comment = 'Z coordinate of ' // cmnt) + rtemp(:) = cache%psi + call h5_add(h5id_root, grp // '/psi', rtemp, & + lbound(cache), ubound(cache), unit = 'Mx', & + comment = 'poloidal flux at ' // cmnt) + rtemp(:) = cache%theta + call h5_add(h5id_root, grp // '/theta', rtemp, & + lbound(cache), ubound(cache), unit = 'rad', & + comment = 'flux poloidal angle at ' // cmnt) + rtemp(:) = cache%sqrt_g + call h5_add(h5id_root, grp // '/sqrt_g', rtemp, & + lbound(cache), ubound(cache), unit = 'cm G^-1', & + comment = 'Jacobian at ' // cmnt) + rtemp(:) = cache%B0_R + call h5_add(h5id_root, grp // '/B0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'R component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%B0_phi + call h5_add(h5id_root, grp // '/B0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%B0_Z + call h5_add(h5id_root, grp // '/B0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dR_dtheta + call h5_add(h5id_root, grp // '/dR_dtheta', rtemp, & + lbound(cache), ubound(cache), unit = 'cm rad^-1', & + comment = 'Jacobian element (R, theta) at ' // cmnt) + rtemp(:) = cache%dZ_dtheta + call h5_add(h5id_root, grp // '/dZ_dtheta', rtemp, & + lbound(cache), ubound(cache), unit = 'cm rad^-1', & + comment = 'Jacobian element (Z, theta) at ' // cmnt) + call h5_close(h5id_root) + deallocate(itemp, rtemp) + end subroutine coord_cache_write_1 + + subroutine coord_cache_write_2(cache, file, group, comment) + use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close + type(coord_cache_t), dimension(:, :), intent(in) :: cache + character(len = *), intent(in) :: file + character(len = *), intent(in) :: group + character(len = *), intent(in) :: comment + character(len = len_trim(group)) :: grp + character(len = len_trim(comment)) :: cmnt + integer(HID_T) :: h5id_root + integer, dimension(:, :), allocatable :: itemp + real(dp), dimension(:, :), allocatable :: rtemp - subroutine coord_cache_read(cache, file, dataset) + grp = trim(group) + cmnt = trim(comment) + allocate(itemp(size(cache, 1), size(cache, 2)), rtemp(size(cache, 1), size(cache, 2))) + call h5_open_rw(file, h5id_root) + call h5_create_parent_groups(h5id_root, grp // '/') + itemp(:, :) = cache%ktri + call h5_add(h5id_root, grp // '/ktri', itemp, & + lbound(cache), ubound(cache), & + comment = 'triangle index of ' // cmnt) + rtemp(:, :) = cache%R + call h5_add(h5id_root, grp // '/R', rtemp, & + lbound(cache), ubound(cache), unit = 'cm', & + comment = 'R coordinate of ' // cmnt) + rtemp(:, :) = cache%Z + call h5_add(h5id_root, grp // '/Z', rtemp, & + lbound(cache), ubound(cache), unit = 'cm', & + comment = 'Z coordinate of ' // cmnt) + rtemp(:, :) = cache%psi + call h5_add(h5id_root, grp // '/psi', rtemp, & + lbound(cache), ubound(cache), unit = 'Mx', & + comment = 'poloidal flux at ' // cmnt) + rtemp(:, :) = cache%theta + call h5_add(h5id_root, grp // '/theta', rtemp, & + lbound(cache), ubound(cache), unit = 'rad', & + comment = 'flux poloidal angle at ' // cmnt) + rtemp(:, :) = cache%sqrt_g + call h5_add(h5id_root, grp // '/sqrt_g', rtemp, & + lbound(cache), ubound(cache), unit = 'cm G^-1', & + comment = 'Jacobian at ' // cmnt) + rtemp(:, :) = cache%B0_R + call h5_add(h5id_root, grp // '/B0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'R component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%B0_phi + call h5_add(h5id_root, grp // '/B0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%B0_Z + call h5_add(h5id_root, grp // '/B0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dR_dtheta + call h5_add(h5id_root, grp // '/dR_dtheta', rtemp, & + lbound(cache), ubound(cache), unit = 'cm rad^-1', & + comment = 'Jacobian element (R, theta) at ' // cmnt) + rtemp(:, :) = cache%dZ_dtheta + call h5_add(h5id_root, grp // '/dZ_dtheta', rtemp, & + lbound(cache), ubound(cache), unit = 'cm rad^-1', & + comment = 'Jacobian element (Z, theta) at ' // cmnt) + call h5_close(h5id_root) + deallocate(itemp, rtemp) + end subroutine coord_cache_write_2 + + subroutine coord_cache_read_1(cache, file, group) use hdf5_tools, only: HID_T, h5_open, h5_get, h5_close - class(coord_cache_t), dimension(..), intent(inout) :: cache + type(coord_cache_t), dimension(:), intent(inout) :: cache character(len = *), intent(in) :: file - character(len = *), intent(in) :: dataset + character(len = *), intent(in) :: group + character(len = len_trim(group)) :: grp integer(HID_T) :: h5id_root + integer, dimension(:), allocatable :: itemp + real(dp), dimension(:), allocatable :: rtemp - select rank (cache) - rank (2) - call h5_open(file, h5id_root) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/ktri', cache(:, :)%ktri) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/R', cache(:, :)%R) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/Z', cache(:, :)%Z) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:, :)%psi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:, :)%theta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/sqrt_g', cache(:, :)%sqrt_g) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:, :)%B0_R) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:, :)%B0_phi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:, :)%B0_Z) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dR_dtheta', cache(:, :)%dR_dtheta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dZ_dtheta', cache(:, :)%dZ_dtheta) - call h5_close(h5id_root) - rank (1) - call h5_open(file, h5id_root) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/ktri', cache(:)%ktri) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/R', cache(:)%R) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/Z', cache(:)%Z) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:)%psi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:)%theta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/sqrt_g', cache(:)%sqrt_g) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:)%B0_R) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:)%B0_phi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:)%B0_Z) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dR_dtheta', cache(:)%dR_dtheta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dZ_dtheta', cache(:)%dZ_dtheta) - call h5_close(h5id_root) - rank default - error stop 'coord_cache_read: rank-1 or rank-2 array expected for argument ''cache''' - end select - end subroutine coord_cache_read + grp = trim(group) + allocate(itemp(size(cache)), rtemp(size(cache))) + call h5_open(file, h5id_root) + call h5_get(h5id_root, grp // '/ktri', itemp) + cache%ktri = itemp + call h5_get(h5id_root, grp // '/R', rtemp) + cache%R = rtemp + call h5_get(h5id_root, grp // '/Z', rtemp) + cache%Z = rtemp + call h5_get(h5id_root, grp // '/psi', rtemp) + cache%psi = rtemp + call h5_get(h5id_root, grp // '/theta', rtemp) + cache%theta = rtemp + call h5_get(h5id_root, grp // '/sqrt_g', rtemp) + cache%sqrt_g = rtemp + call h5_get(h5id_root, grp // '/B0_R', rtemp) + cache%B0_R = rtemp + call h5_get(h5id_root, grp // '/B0_phi', rtemp) + cache%B0_phi = rtemp + call h5_get(h5id_root, grp // '/B0_Z', rtemp) + cache%B0_Z = rtemp + call h5_get(h5id_root, grp // '/dR_dtheta', rtemp) + cache%dR_dtheta = rtemp + call h5_get(h5id_root, grp // '/dZ_dtheta', rtemp) + cache%dZ_dtheta = rtemp + call h5_close(h5id_root) + deallocate(itemp, rtemp) + end subroutine coord_cache_read_1 - subroutine field_cache_write(cache, file, dataset, comment) + subroutine coord_cache_read_2(cache, file, group) + use hdf5_tools, only: HID_T, h5_open, h5_get, h5_close + type(coord_cache_t), dimension(:, :), intent(inout) :: cache + character(len = *), intent(in) :: file + character(len = *), intent(in) :: group + character(len = len_trim(group)) :: grp + integer(HID_T) :: h5id_root + integer, dimension(:, :), allocatable :: itemp + real(dp), dimension(:, :), allocatable :: rtemp + + grp = trim(group) + allocate(itemp(size(cache, 1), size(cache, 2)), rtemp(size(cache, 1), size(cache, 2))) + call h5_open(file, h5id_root) + call h5_get(h5id_root, grp // '/ktri', itemp) + cache%ktri = itemp + call h5_get(h5id_root, grp // '/R', rtemp) + cache%R = rtemp + call h5_get(h5id_root, grp // '/Z', rtemp) + cache%Z = rtemp + call h5_get(h5id_root, grp // '/psi', rtemp) + cache%psi = rtemp + call h5_get(h5id_root, grp // '/theta', rtemp) + cache%theta = rtemp + call h5_get(h5id_root, grp // '/sqrt_g', rtemp) + cache%sqrt_g = rtemp + call h5_get(h5id_root, grp // '/B0_R', rtemp) + cache%B0_R = rtemp + call h5_get(h5id_root, grp // '/B0_phi', rtemp) + cache%B0_phi = rtemp + call h5_get(h5id_root, grp // '/B0_Z', rtemp) + cache%B0_Z = rtemp + call h5_get(h5id_root, grp // '/dR_dtheta', rtemp) + cache%dR_dtheta = rtemp + call h5_get(h5id_root, grp // '/dZ_dtheta', rtemp) + cache%dZ_dtheta = rtemp + call h5_close(h5id_root) + deallocate(itemp, rtemp) + end subroutine coord_cache_read_2 + + subroutine field_cache_write_1(cache, file, group, comment) use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close - class(field_cache_t), dimension(..), intent(in) :: cache + type(field_cache_t), dimension(:), intent(in) :: cache character(len = *), intent(in) :: file - character(len = *), intent(in) :: dataset + character(len = *), intent(in) :: group character(len = *), intent(in) :: comment + character(len = len_trim(group)) :: grp + character(len = len_trim(comment)) :: cmnt integer(HID_T) :: h5id_root + real(dp), dimension(:), allocatable :: rtemp + grp = trim(group) + cmnt = trim(comment) + allocate(rtemp(size(cache))) call h5_open_rw(file, h5id_root) - call h5_create_parent_groups(h5id_root, trim(adjustl(dataset)) // '/') - select rank (cache) - rank (1) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:)%psi, & - lbound(cache), ubound(cache), unit = 'Mx', & - comment = 'poloidal flux at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:)%theta, & - lbound(cache), ubound(cache), unit = 'rad', & - comment = 'flux poloidal angle at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:)%B0(1), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:)%B0(3), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:)%B0(2), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_R', cache(:)%j0(1), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_Z', cache(:)%j0(3), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'Z component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_phi', cache(:)%j0(2), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0', cache(:)%Bmod, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0_dR', cache(:)%dBmod_dR, & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0_dZ', cache(:)%dBmod_dZ, & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0R_dR', cache(:)%dB0_dR(1), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0R_dZ', cache(:)%dB0_dZ(1), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dR', cache(:)%dB0_dR(2), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dZ', cache(:)%dB0_dZ(2), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dR', cache(:)%dB0_dR(3), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dZ', cache(:)%dB0_dZ(3), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0R_dR', cache(:)%dj0_dR(1), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0R_dZ', cache(:)%dj0_dZ(1), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dR', cache(:)%dj0_dR(2), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dZ', cache(:)%dj0_dZ(2), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dR', cache(:)%dj0_dR(3), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of Z component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dZ', cache(:)%dj0_dZ(3), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of Z component of equilibrium current density at ' // trim(adjustl(comment))) - rank (2) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:, :)%psi, & - lbound(cache), ubound(cache), unit = 'Mx', & - comment = 'poloidal flux at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:, :)%theta, & - lbound(cache), ubound(cache), unit = 'rad', & - comment = 'flux poloidal angle at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:, :)%B0(1), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:, :)%B0(3), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:, :)%B0(2), & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_R', cache(:, :)%j0(1), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_Z', cache(:, :)%j0(3), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'Z component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/j0_phi', cache(:, :)%j0(2), & - lbound(cache), ubound(cache), unit = 'statA', & - comment = 'physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/B0', cache(:, :)%Bmod, & - lbound(cache), ubound(cache), unit = 'G', & - comment = 'magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0_dR', cache(:, :)%dBmod_dR, & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0_dZ', cache(:, :)%dBmod_dZ, & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of magnitude of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0R_dR', cache(:, :)%dB0_dR(1), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0R_dZ', cache(:, :)%dB0_dZ(1), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of R component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dR', cache(:, :)%dB0_dR(2), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dZ', cache(:, :)%dB0_dZ(2), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of physical phi component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dR', cache(:, :)%dB0_dR(3), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'R derivative of Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dZ', cache(:, :)%dB0_dZ(3), & - lbound(cache), ubound(cache), unit = 'G cm^-1', & - comment = 'Z derivative of Z component of equilibrium magnetic field at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0R_dR', cache(:, :)%dj0_dR(1), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0R_dZ', cache(:, :)%dj0_dZ(1), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of R component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dR', cache(:, :)%dj0_dR(2), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dZ', cache(:, :)%dj0_dZ(2), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of physical phi component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dR', cache(:, :)%dj0_dR(3), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'R derivative of Z component of equilibrium current density at ' // trim(adjustl(comment))) - call h5_add(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dZ', cache(:, :)%dj0_dZ(3), & - lbound(cache), ubound(cache), unit = 'statA cm^-1', & - comment = 'Z derivative of Z component of equilibrium current density at ' // trim(adjustl(comment))) - rank default - error stop 'field_cache_write: rank-1 or rank-2 array expected for argument ''cache''' - end select + call h5_create_parent_groups(h5id_root, grp // '/') + rtemp(:) = cache%psi + call h5_add(h5id_root, grp // '/psi', rtemp, & + lbound(cache), ubound(cache), unit = 'Mx', & + comment = 'poloidal flux at ' // cmnt) + rtemp(:) = cache%theta + call h5_add(h5id_root, grp // '/theta', rtemp, & + lbound(cache), ubound(cache), unit = 'rad', & + comment = 'flux poloidal angle at ' // cmnt) + rtemp(:) = cache%B0(1) + call h5_add(h5id_root, grp // '/B0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'R component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%B0(3) + call h5_add(h5id_root, grp // '/B0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%B0(2) + call h5_add(h5id_root, grp // '/B0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%j0(1) + call h5_add(h5id_root, grp // '/j0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'R component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%j0(3) + call h5_add(h5id_root, grp // '/j0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'Z component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%j0(2) + call h5_add(h5id_root, grp // '/j0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'physical phi component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%Bmod + call h5_add(h5id_root, grp // '/B0', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dBmod_dR + call h5_add(h5id_root, grp // '/dB0_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dBmod_dZ + call h5_add(h5id_root, grp // '/dB0_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dR(1) + call h5_add(h5id_root, grp // '/dB0R_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of R component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dZ(1) + call h5_add(h5id_root, grp // '/dB0R_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of R component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dR(2) + call h5_add(h5id_root, grp // '/dB0phi_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dZ(2) + call h5_add(h5id_root, grp // '/dB0phi_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dR(3) + call h5_add(h5id_root, grp // '/dB0Z_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dB0_dZ(3) + call h5_add(h5id_root, grp // '/dB0Z_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:) = cache%dj0_dR(1) + call h5_add(h5id_root, grp // '/dj0R_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of R component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%dj0_dZ(1) + call h5_add(h5id_root, grp // '/dj0R_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of R component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%dj0_dR(2) + call h5_add(h5id_root, grp // '/dj0phi_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of physical phi component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%dj0_dZ(2) + call h5_add(h5id_root, grp // '/dj0phi_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of physical phi component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%dj0_dR(3) + call h5_add(h5id_root, grp // '/dj0Z_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of Z component of equilibrium current density at ' // cmnt) + rtemp(:) = cache%dj0_dZ(3) + call h5_add(h5id_root, grp // '/dj0Z_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of Z component of equilibrium current density at ' // cmnt) call h5_close(h5id_root) - end subroutine field_cache_write + deallocate(rtemp) + end subroutine field_cache_write_1 - subroutine field_cache_read(cache, file, dataset) + subroutine field_cache_write_2(cache, file, group, comment) + use hdf5_tools, only: HID_T, h5_open_rw, h5_create_parent_groups, h5_add, h5_close + type(field_cache_t), dimension(:, :), intent(in) :: cache + character(len = *), intent(in) :: file + character(len = *), intent(in) :: group + character(len = *), intent(in) :: comment + character(len = len_trim(group)) :: grp + character(len = len_trim(comment)) :: cmnt + integer(HID_T) :: h5id_root + real(dp), dimension(:, :), allocatable :: rtemp + + grp = trim(group) + cmnt = trim(comment) + allocate(rtemp(size(cache, 1), size(cache, 2))) + call h5_open_rw(file, h5id_root) + call h5_create_parent_groups(h5id_root, grp // '/') + rtemp(:, :) = cache%psi + call h5_add(h5id_root, grp // '/psi', rtemp, & + lbound(cache), ubound(cache), unit = 'Mx', & + comment = 'poloidal flux at ' // cmnt) + rtemp(:, :) = cache%theta + call h5_add(h5id_root, grp // '/theta', rtemp, & + lbound(cache), ubound(cache), unit = 'rad', & + comment = 'flux poloidal angle at ' // cmnt) + rtemp(:, :) = cache%B0(1) + call h5_add(h5id_root, grp // '/B0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'R component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%B0(3) + call h5_add(h5id_root, grp // '/B0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%B0(2) + call h5_add(h5id_root, grp // '/B0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%j0(1) + call h5_add(h5id_root, grp // '/j0_R', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'R component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%j0(3) + call h5_add(h5id_root, grp // '/j0_Z', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'Z component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%j0(2) + call h5_add(h5id_root, grp // '/j0_phi', rtemp, & + lbound(cache), ubound(cache), unit = 'statA', & + comment = 'physical phi component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%Bmod + call h5_add(h5id_root, grp // '/B0', rtemp, & + lbound(cache), ubound(cache), unit = 'G', & + comment = 'magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dBmod_dR + call h5_add(h5id_root, grp // '/dB0_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dBmod_dZ + call h5_add(h5id_root, grp // '/dB0_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of magnitude of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dR(1) + call h5_add(h5id_root, grp // '/dB0R_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of R component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dZ(1) + call h5_add(h5id_root, grp // '/dB0R_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of R component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dR(2) + call h5_add(h5id_root, grp // '/dB0phi_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dZ(2) + call h5_add(h5id_root, grp // '/dB0phi_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of physical phi component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dR(3) + call h5_add(h5id_root, grp // '/dB0Z_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'R derivative of Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dB0_dZ(3) + call h5_add(h5id_root, grp // '/dB0Z_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'G cm^-1', & + comment = 'Z derivative of Z component of equilibrium magnetic field at ' // cmnt) + rtemp(:, :) = cache%dj0_dR(1) + call h5_add(h5id_root, grp // '/dj0R_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of R component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%dj0_dZ(1) + call h5_add(h5id_root, grp // '/dj0R_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of R component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%dj0_dR(2) + call h5_add(h5id_root, grp // '/dj0phi_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of physical phi component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%dj0_dZ(2) + call h5_add(h5id_root, grp // '/dj0phi_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of physical phi component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%dj0_dR(3) + call h5_add(h5id_root, grp // '/dj0Z_dR', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'R derivative of Z component of equilibrium current density at ' // cmnt) + rtemp(:, :) = cache%dj0_dZ(3) + call h5_add(h5id_root, grp // '/dj0Z_dZ', rtemp, & + lbound(cache), ubound(cache), unit = 'statA cm^-1', & + comment = 'Z derivative of Z component of equilibrium current density at ' // cmnt) + call h5_close(h5id_root) + deallocate(rtemp) + end subroutine field_cache_write_2 + + subroutine field_cache_read_1(cache, file, group) use hdf5_tools, only: HID_T, h5_open, h5_get, h5_close - class(field_cache_t), dimension(..), intent(inout) :: cache + type(field_cache_t), dimension(:), intent(inout) :: cache character(len = *), intent(in) :: file - character(len = *), intent(in) :: dataset + character(len = *), intent(in) :: group + character(len = len_trim(group)) :: grp integer(HID_T) :: h5id_root + real(dp), dimension(:), allocatable :: rtemp + grp = trim(group) + allocate(rtemp(size(cache))) call h5_open(file, h5id_root) - select rank (cache) - rank (1) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:)%psi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:)%theta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:)%B0(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:)%B0(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:)%B0(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_R', cache(:)%j0(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_Z', cache(:)%j0(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_phi', cache(:)%j0(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0', cache(:)%Bmod) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0_dR', cache(:)%dBmod_dR) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0_dZ', cache(:)%dBmod_dZ) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0R_dR', cache(:)%dB0_dR(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0R_dZ', cache(:)%dB0_dZ(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dR', cache(:)%dB0_dR(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dZ', cache(:)%dB0_dZ(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dR', cache(:)%dB0_dR(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dZ', cache(:)%dB0_dZ(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0R_dR', cache(:)%dj0_dR(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0R_dZ', cache(:)%dj0_dZ(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dR', cache(:)%dj0_dR(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dZ', cache(:)%dj0_dZ(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dR', cache(:)%dj0_dR(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dZ', cache(:)%dj0_dZ(3)) - rank (2) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/psi', cache(:, :)%psi) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/theta', cache(:, :)%theta) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_R', cache(:, :)%B0(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_Z', cache(:, :)%B0(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0_phi', cache(:, :)%B0(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_R', cache(:, :)%j0(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_Z', cache(:, :)%j0(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/j0_phi', cache(:, :)%j0(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/B0', cache(:, :)%Bmod) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0_dR', cache(:, :)%dBmod_dR) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0_dZ', cache(:, :)%dBmod_dZ) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0R_dR', cache(:, :)%dB0_dR(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0R_dZ', cache(:, :)%dB0_dZ(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dR', cache(:, :)%dB0_dR(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0phi_dZ', cache(:, :)%dB0_dZ(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dR', cache(:, :)%dB0_dR(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dB0Z_dZ', cache(:, :)%dB0_dZ(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0R_dR', cache(:, :)%dj0_dR(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0R_dZ', cache(:, :)%dj0_dZ(1)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dR', cache(:, :)%dj0_dR(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0phi_dZ', cache(:, :)%dj0_dZ(2)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dR', cache(:, :)%dj0_dR(3)) - call h5_get(h5id_root, trim(adjustl(dataset)) // '/dj0Z_dZ', cache(:, :)%dj0_dZ(3)) - rank default - error stop 'field_cache_read: rank-1 or rank-2 array expected for argument ''cache''' - end select + call h5_get(h5id_root, grp // '/psi', rtemp) + cache%psi = rtemp + call h5_get(h5id_root, grp // '/theta', rtemp) + cache%theta = rtemp + call h5_get(h5id_root, grp // '/B0_R', rtemp) + cache%B0(1) = rtemp + call h5_get(h5id_root, grp // '/B0_Z', rtemp) + cache%B0(3) = rtemp + call h5_get(h5id_root, grp // '/B0_phi', rtemp) + cache%B0(2) = rtemp + call h5_get(h5id_root, grp // '/j0_R', rtemp) + cache%j0(1) = rtemp + call h5_get(h5id_root, grp // '/j0_Z', rtemp) + cache%j0(3) = rtemp + call h5_get(h5id_root, grp // '/j0_phi', rtemp) + cache%j0(2) = rtemp + call h5_get(h5id_root, grp // '/B0', rtemp) + cache%Bmod = rtemp + call h5_get(h5id_root, grp // '/dB0_dR', rtemp) + cache%dBmod_dR = rtemp + call h5_get(h5id_root, grp // '/dB0_dZ', rtemp) + cache%dBmod_dZ = rtemp + call h5_get(h5id_root, grp // '/dB0R_dR', rtemp) + cache%dB0_dR(1) = rtemp + call h5_get(h5id_root, grp // '/dB0R_dZ', rtemp) + cache%dB0_dZ(1) = rtemp + call h5_get(h5id_root, grp // '/dB0phi_dR', rtemp) + cache%dB0_dR(2) = rtemp + call h5_get(h5id_root, grp // '/dB0phi_dZ', rtemp) + cache%dB0_dZ(2) = rtemp + call h5_get(h5id_root, grp // '/dB0Z_dR', rtemp) + cache%dB0_dR(3) = rtemp + call h5_get(h5id_root, grp // '/dB0Z_dZ', rtemp) + cache%dB0_dZ(3) = rtemp + call h5_get(h5id_root, grp // '/dj0R_dR', rtemp) + cache%dj0_dR(1) = rtemp + call h5_get(h5id_root, grp // '/dj0R_dZ', rtemp) + cache%dj0_dZ(1) = rtemp + call h5_get(h5id_root, grp // '/dj0phi_dR', rtemp) + cache%dj0_dR(2) = rtemp + call h5_get(h5id_root, grp // '/dj0phi_dZ', rtemp) + cache%dj0_dZ(2) = rtemp + call h5_get(h5id_root, grp // '/dj0Z_dR', rtemp) + cache%dj0_dR(3) = rtemp + call h5_get(h5id_root, grp // '/dj0Z_dZ', rtemp) + cache%dj0_dZ(3) = rtemp + call h5_close(h5id_root) + deallocate(rtemp) + end subroutine field_cache_read_1 + + subroutine field_cache_read_2(cache, file, group) + use hdf5_tools, only: HID_T, h5_open, h5_get, h5_close + type(field_cache_t), dimension(:, :), intent(inout) :: cache + character(len = *), intent(in) :: file + character(len = *), intent(in) :: group + character(len = len_trim(group)) :: grp + integer(HID_T) :: h5id_root + real(dp), dimension(:, :), allocatable :: rtemp + + grp = trim(group) + allocate(rtemp(size(cache, 1), size(cache, 2))) + call h5_open(file, h5id_root) + call h5_get(h5id_root, grp // '/psi', rtemp) + cache%psi = rtemp + call h5_get(h5id_root, grp // '/theta', rtemp) + cache%theta = rtemp + call h5_get(h5id_root, grp // '/B0_R', rtemp) + cache%B0(1) = rtemp + call h5_get(h5id_root, grp // '/B0_Z', rtemp) + cache%B0(3) = rtemp + call h5_get(h5id_root, grp // '/B0_phi', rtemp) + cache%B0(2) = rtemp + call h5_get(h5id_root, grp // '/j0_R', rtemp) + cache%j0(1) = rtemp + call h5_get(h5id_root, grp // '/j0_Z', rtemp) + cache%j0(3) = rtemp + call h5_get(h5id_root, grp // '/j0_phi', rtemp) + cache%j0(2) = rtemp + call h5_get(h5id_root, grp // '/B0', rtemp) + cache%Bmod = rtemp + call h5_get(h5id_root, grp // '/dB0_dR', rtemp) + cache%dBmod_dR = rtemp + call h5_get(h5id_root, grp // '/dB0_dZ', rtemp) + cache%dBmod_dZ = rtemp + call h5_get(h5id_root, grp // '/dB0R_dR', rtemp) + cache%dB0_dR(1) = rtemp + call h5_get(h5id_root, grp // '/dB0R_dZ', rtemp) + cache%dB0_dZ(1) = rtemp + call h5_get(h5id_root, grp // '/dB0phi_dR', rtemp) + cache%dB0_dR(2) = rtemp + call h5_get(h5id_root, grp // '/dB0phi_dZ', rtemp) + cache%dB0_dZ(2) = rtemp + call h5_get(h5id_root, grp // '/dB0Z_dR', rtemp) + cache%dB0_dR(3) = rtemp + call h5_get(h5id_root, grp // '/dB0Z_dZ', rtemp) + cache%dB0_dZ(3) = rtemp + call h5_get(h5id_root, grp // '/dj0R_dR', rtemp) + cache%dj0_dR(1) = rtemp + call h5_get(h5id_root, grp // '/dj0R_dZ', rtemp) + cache%dj0_dZ(1) = rtemp + call h5_get(h5id_root, grp // '/dj0phi_dR', rtemp) + cache%dj0_dR(2) = rtemp + call h5_get(h5id_root, grp // '/dj0phi_dZ', rtemp) + cache%dj0_dZ(2) = rtemp + call h5_get(h5id_root, grp // '/dj0Z_dR', rtemp) + cache%dj0_dR(3) = rtemp + call h5_get(h5id_root, grp // '/dj0Z_dZ', rtemp) + cache%dj0_dZ(3) = rtemp call h5_close(h5id_root) - end subroutine field_cache_read + deallocate(rtemp) + end subroutine field_cache_read_2 subroutine shielding_init(s, m) type(shielding_t), intent(inout) :: s From 07f492e2359bb4326325e60b418562f3bbb561bf Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 24 Apr 2025 14:28:10 +0200 Subject: [PATCH 3/6] remove obsolete shielding projection auxiliaries --- src/mephit_mesh.F90 | 90 +++------------------------------------------ 1 file changed, 5 insertions(+), 85 deletions(-) diff --git a/src/mephit_mesh.F90 b/src/mephit_mesh.F90 index 53738ea..39f4a06 100644 --- a/src/mephit_mesh.F90 +++ b/src/mephit_mesh.F90 @@ -330,11 +330,6 @@ end subroutine rtree_query !> Free parameter in the compensated scheme for shielding real(dp) :: coeff - - !> pressure perturbation projection - integer :: inner_kp_max, outer_kp_max - integer, allocatable :: kpois(:, :) - real(dp), allocatable :: weights(:, :) end type shielding_t type :: cache_t @@ -1030,29 +1025,12 @@ subroutine field_cache_read_2(cache, file, group) deallocate(rtemp) end subroutine field_cache_read_2 - subroutine shielding_init(s, m) - type(shielding_t), intent(inout) :: s - integer, intent(in) :: m - integer :: kf - - call shielding_deinit(s) - kf = mesh%res_ind(m) - s%inner_kp_max = mesh%kp_max(kf - 1) - s%outer_kp_max = mesh%kp_max(kf) - allocate(s%kpois(2, s%inner_kp_max + s%outer_kp_max)) - allocate(s%weights(2, s%inner_kp_max + s%outer_kp_max)) - end subroutine shielding_init - - subroutine shielding_deinit(s) + elemental subroutine shielding_deinit(s) type(shielding_t), intent(inout) :: s s%coeff = 0d0 - s%inner_kp_max = 0 - s%outer_kp_max = 0 if (allocated(s%GL_weights)) deallocate(s%GL_weights) if (allocated(s%sample_Ires)) deallocate(s%sample_Ires) - if (allocated(s%kpois)) deallocate(s%kpois) - if (allocated(s%weights)) deallocate(s%weights) end subroutine shielding_deinit subroutine shielding_write(shielding, file, group) @@ -1072,16 +1050,6 @@ subroutine shielding_write(shielding, file, group) unit = 'Mx', comment = 'G-L quadrature weights including psi interval') call h5_add(h5id_root, grp // '/coeff', shielding%coeff, unit = 'g^-1 cm s', & comment = 'coefficient in the compensated scheme for shielding') - call h5_add(h5id_root, grp // '/inner_kp_max', shielding%inner_kp_max, & - comment = 'number of points on inner flux surface') - call h5_add(h5id_root, grp // '/outer_kp_max', shielding%outer_kp_max, & - comment = 'number of points on outer flux surface') - call h5_add(h5id_root, grp // '/kpois', shielding%kpois, & - lbound(shielding%kpois), ubound(shielding%kpois), & - comment = 'point indices of surrounding points of theta projection') - call h5_add(h5id_root, grp // '/weights', shielding%weights, & - lbound(shielding%weights), ubound(shielding%weights), & - unit = '1', comment = 'relative weighting of surrounding points of theta projection') call h5_close(h5id_root) end subroutine shielding_write @@ -1097,10 +1065,6 @@ subroutine shielding_read(shielding, file, group) grp = trim(group) call h5_open(file, h5id_root) call h5_get(h5id_root, grp // '/coeff', shielding%coeff) - call h5_get(h5id_root, grp // '/inner_kp_max', shielding%inner_kp_max) - call h5_get(h5id_root, grp // '/outer_kp_max', shielding%outer_kp_max) - call h5_get(h5id_root, grp // '/kpois', shielding%kpois) - call h5_get(h5id_root, grp // '/weights', shielding%weights) call h5_get_bounds(h5id_root, grp // '/sample_Ires/R', lb1, lb2, ub1, ub2) if (allocated(shielding%GL_weights)) deallocate(shielding%GL_weights) allocate(shielding%GL_weights(lb2:ub2)) @@ -1114,7 +1078,6 @@ end subroutine shielding_read subroutine cache_init(cache, GL_order) type(cache_t), intent(inout) :: cache integer, intent(in) :: GL_order - integer :: m call cache_deinit(cache) cache%GL_order = GL_order @@ -1122,25 +1085,19 @@ subroutine cache_init(cache, GL_order) allocate(cache%sample_polmodes_half(mesh%ntri)) allocate(cache%sample_jnperp(mesh%npoint - 1)) allocate(cache%shielding(mesh%m_res_min:mesh%m_res_max)) - do m = mesh%m_res_min, mesh%m_res_max - call shielding_init(cache%shielding(m), m) - end do allocate(cache%edge_fields(mesh%GL_order, mesh%nedge), cache%area_fields(mesh%GL2_order, mesh%ntri)) allocate(cache%mid_fields(mesh%nedge), cache%cntr_fields(mesh%ntri)) end subroutine cache_init subroutine cache_deinit(cache) type(cache_t), intent(inout) :: cache - integer :: m cache%GL_order = 0 if (allocated(cache%sample_polmodes)) deallocate(cache%sample_polmodes) if (allocated(cache%sample_polmodes_half)) deallocate(cache%sample_polmodes_half) if (allocated(cache%sample_jnperp)) deallocate(cache%sample_jnperp) if (allocated(cache%shielding)) then - do m = lbound(cache%shielding, 1), ubound(cache%shielding, 1) - call shielding_deinit(cache%shielding(m)) - end do + call shielding_deinit(cache%shielding) deallocate(cache%shielding) end if if (allocated(cache%edge_fields)) deallocate(cache%edge_fields) @@ -1241,9 +1198,9 @@ subroutine generate_mesh call compute_sample_polmodes(cache%sample_polmodes, .false.) call compute_sample_jnperp(cache%sample_jnperp) do m = mesh%m_res_min, mesh%m_res_max + call compute_shielding_auxiliaries(cache%shielding(m), m) call compute_sample_Ires(cache%shielding(m)%sample_Ires, & cache%shielding(m)%GL_weights, cache%GL_order, m) - call compute_shielding_auxiliaries(cache%shielding(m), m) end do call compute_kilca_auxiliaries call compute_gpec_jacfac @@ -2397,54 +2354,17 @@ subroutine check_mesh end subroutine check_mesh subroutine compute_shielding_auxiliaries(s, m) - use mephit_util, only: pi, clight, resample1d, interp1d, binsearch + use mephit_util, only: pi, clight, interp1d type(shielding_t), intent(inout) :: s integer, intent(in) :: m - integer :: kf, inner_kp_low, outer_kp_low + integer :: kf real(dp) :: dq_dpsi kf = mesh%res_ind(m) - inner_kp_low = mesh%kp_low(kf - 1) - outer_kp_low = mesh%kp_low(kf) - s%kpois(:, :) = 0 - s%weights(:, :) = 0d0 - call project_theta(mesh%node_theta_flux, mesh%node_theta_geom, & - s%outer_kp_max, outer_kp_low, s%inner_kp_max, inner_kp_low, & - s%kpois(:, :s%inner_kp_max), s%weights(:, :s%inner_kp_max)) - call project_theta(mesh%node_theta_flux, mesh%node_theta_geom, & - s%inner_kp_max, inner_kp_low, s%outer_kp_max, outer_kp_low, & - s%kpois(:, s%inner_kp_max + 1:), s%weights(:, s%inner_kp_max + 1:)) dq_dpsi = interp1d(equil%psi_eqd, equil%qpsi, fs_half%psi(kf), 3, .true.) s%coeff = clight * mesh%n / (4d0 * pi * mesh%R_O) * & abs(dq_dpsi / (fs_half%q(kf) * fs_half%dp_dpsi(kf))) / & (mesh%n * abs(fs%q(kf) - fs%q(kf - 1))) - - contains - subroutine project_theta(theta_flux, theta_geom, dest_kp_max, dest_kpoi_low, & - src_kp_max, src_kpoi_low, kpois, weights) - real(dp), dimension(:), intent(in) :: theta_flux, theta_geom - integer, intent(in) :: dest_kp_max, dest_kpoi_low, src_kp_max, src_kpoi_low - integer, dimension(:, :), intent(out) :: kpois - real(dp), dimension(:, :), intent(out) :: weights - real(dp), dimension(-1:dest_kp_max+2) :: theta_flux_ext, theta_geom_ext - real(dp), dimension(src_kp_max) :: src_theta, dest_theta - integer :: kp, inf_kp - - call extend_over_branch(theta_flux(dest_kpoi_low + 1:& - dest_kpoi_low + dest_kp_max), 2, theta_flux_ext) - call extend_over_branch(theta_geom(dest_kpoi_low + 1:& - dest_kpoi_low + dest_kp_max), 2, theta_geom_ext) - src_theta = theta_flux(src_kpoi_low + 1:src_kpoi_low + src_kp_max) - call resample1d(theta_flux_ext, theta_geom_ext, src_theta, dest_theta, 3) - kpois(:, 1) = [dest_kpoi_low + 1, dest_kpoi_low + 2] - weights(:, 1) = [1d0, 0d0] - do kp = 2, src_kp_max - call binsearch(theta_geom_ext(1:dest_kp_max + 1), 0, dest_theta(kp), inf_kp) - kpois(:, kp) = [dest_kpoi_low + inf_kp, dest_kpoi_low + mod(inf_kp, dest_kp_max) + 1] - weights(:, kp) = abs(dest_theta(kp) - theta_geom_ext(inf_kp:inf_kp + 1)) / & - (theta_geom_ext(inf_kp + 1) - theta_geom_ext(inf_kp)) - end do - end subroutine project_theta end subroutine compute_shielding_auxiliaries subroutine compute_gpec_jacfac From 8a734bf0fcd82a83cebc72ad1616463257e1e6c1 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 24 Apr 2025 16:21:20 +0200 Subject: [PATCH 4/6] specify FLR2 coefficient units --- src/mephit_flr2.f90 | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/src/mephit_flr2.f90 b/src/mephit_flr2.f90 index 0f6f8f0..2a9c385 100644 --- a/src/mephit_flr2.f90 +++ b/src/mephit_flr2.f90 @@ -102,44 +102,42 @@ subroutine flr2_write(flr2, file, dataset) call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', & flr2%psi, lbound(flr2%psi), ubound(flr2%psi), unit = 'Mx', & comment = 'poloidal flux') - call h5_add(h5id_root, trim(adjustl(dataset)) // '/psi', & - flr2%psi, lbound(flr2%psi), ubound(flr2%psi), unit = '?', & - comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/a2_in', & - flr2%a2_in, lbound(flr2%a2_in), ubound(flr2%a2_in), unit = '?', & + flr2%a2_in, lbound(flr2%a2_in), ubound(flr2%a2_in), unit = '1', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/a2_out', & - flr2%a2_out, lbound(flr2%a2_out), ubound(flr2%a2_out), unit = '?', & + flr2%a2_out, lbound(flr2%a2_out), ubound(flr2%a2_out), unit = '1', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/a0', & - flr2%a0, lbound(flr2%a0), ubound(flr2%a0), unit = '?', & + flr2%a0, lbound(flr2%a0), ubound(flr2%a0), unit = 'Mx^-2', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/b2_in', & - flr2%b2_in, lbound(flr2%b2_in), ubound(flr2%b2_in), unit = '?', & + flr2%b2_in, lbound(flr2%b2_in), ubound(flr2%b2_in), unit = 'statV Mx^-1', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/b2_out', & - flr2%b2_out, lbound(flr2%b2_out), ubound(flr2%b2_out), unit = '?', & + flr2%b2_out, lbound(flr2%b2_out), ubound(flr2%b2_out), unit = 'statV Mx^-1', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/b0', & - flr2%b0, lbound(flr2%b0), ubound(flr2%b0), unit = '?', & + flr2%b0, lbound(flr2%b0), ubound(flr2%b0), unit = 'statV Mx^-3', & comment = 'Poisson equation coefficient') + ! note: statA with c_0 = 1 call h5_add(h5id_root, trim(adjustl(dataset)) // '/c2_in', & - flr2%c2_in, lbound(flr2%c2_in), ubound(flr2%c2_in), unit = '?', & + flr2%c2_in, lbound(flr2%c2_in), ubound(flr2%c2_in), unit = 'statA statV^-1 Mx', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/c2_out', & - flr2%c2_out, lbound(flr2%c2_out), ubound(flr2%c2_out), unit = '?', & + flr2%c2_out, lbound(flr2%c2_out), ubound(flr2%c2_out), unit = 'statA statV^-1 Mx', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/c0', & - flr2%c0, lbound(flr2%c0), ubound(flr2%c0), unit = '?', & + flr2%c0, lbound(flr2%c0), ubound(flr2%c0), unit = 'statA statV^-1 M^-1', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/d2_in', & - flr2%d2_in, lbound(flr2%d2_in), ubound(flr2%d2_in), unit = '?', & + flr2%d2_in, lbound(flr2%d2_in), ubound(flr2%d2_in), unit = 'statA', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/d2_out', & - flr2%d2_out, lbound(flr2%d2_out), ubound(flr2%d2_out), unit = '?', & + flr2%d2_out, lbound(flr2%d2_out), ubound(flr2%d2_out), unit = 'statA', & comment = 'Poisson equation coefficient') call h5_add(h5id_root, trim(adjustl(dataset)) // '/d0', & - flr2%d0, lbound(flr2%d0), ubound(flr2%d0), unit = '?', & + flr2%d0, lbound(flr2%d0), ubound(flr2%d0), unit = 'statA Mx^-2', & comment = 'Poisson equation coefficient') call h5_close(h5id_root) end subroutine flr2_write From 769df29b1eebbe530d9bf096da4c5535d26b7ff0 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Thu, 24 Apr 2025 16:39:05 +0200 Subject: [PATCH 5/6] default shielding factor to 1, apply only if refined --- data/mephit_g000001.0001_TCFP_hip.in | 27 +++++++++++---- data/mephit_g33353_2900_EQH.in | 9 ----- data/mephit_illustration.in | 50 +++++++++++++--------------- src/mephit_conf.f90 | 6 ++-- src/mephit_iter.F90 | 2 +- 5 files changed, 49 insertions(+), 45 deletions(-) diff --git a/data/mephit_g000001.0001_TCFP_hip.in b/data/mephit_g000001.0001_TCFP_hip.in index fa36e78..f38e00c 100644 --- a/data/mephit_g000001.0001_TCFP_hip.in +++ b/data/mephit_g000001.0001_TCFP_hip.in @@ -13,11 +13,29 @@ config%refinement_scheme = 0 ! 0: iMHD, 1: FLRE 2nd order config%currn_model = 0 + ! 0: Viktor Nemov's code, 1: GPEC, 2: pre-computed Fourier transform + config%vac_src = 2 + ! prefactor in Biot-Savart law, e.g. 5 windings and conversion from A to statA / c_0 + config%Biot_Savart_prefactor = 5.0d-1 + !> file containing coil field data + config%coil_file = '/proj/plasma/CODE/lainer_p/MEPHIT/run/AUG_B_coils.h5' + !> file containing coil currents + config%currents_file = '/proj/plasma/DATA/BALANCE/COIL/33353/33353.2900_coil_markl.dat' + !> file containing equilibrium density profile + config%dens_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_ne_PED_MMARKL_rho_pol.dat' + !> file containing equilibrium electron temperature profile + config%temp_e_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Te_PED_MMARKL_rho_pol.dat' + !> file containing equilibrium ion temperature profile + config%temp_i_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Ti_PED_MMARKL_rho_pol.dat' + !> file containing equilibrium radial electric field profile + config%E_r_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353_2900_Er.dat' + ! mass number of ions + config%m_i = 2.0 ! toroidal mode number of perturbation config%n = 2 - ! maximum number of iteration steps + ! maximum number of (preconditioned) fixed-point iteration steps config%niter = 50 - ! maximum relative error in fixed-point iteration + ! maximum relative error in (preconditioned) fixed-point iteration config%iter_rel_err = 1d-12 ! maximum number of iterations in Arnoldi method config%nkrylov = 100 @@ -49,8 +67,5 @@ &arrays ! relative refinement around resonances refinement(4) = 1.414213562373095d0 - ! deleted flux surfaces around resonances - deletions(4) = 1 - ! coefficients for sheet currents - sheet_current_factor(4) = 1d4 + Delta_rad_res(4) = 0.1125d0 / diff --git a/data/mephit_g33353_2900_EQH.in b/data/mephit_g33353_2900_EQH.in index 71bcb3b..a2fb6ce 100644 --- a/data/mephit_g33353_2900_EQH.in +++ b/data/mephit_g33353_2900_EQH.in @@ -90,13 +90,4 @@ Delta_rad_res(9) = 5.347772320733589e-02 Delta_rad_res(10) = 6.066328245075899e-02 Delta_rad_res(11) = 3.464590203985631e-02 - sheet_current_factor(3) = 1d5 - sheet_current_factor(4) = 1d5 - sheet_current_factor(5) = 1d5 - sheet_current_factor(6) = 1d5 - sheet_current_factor(7) = 1d5 - sheet_current_factor(8) = 1d5 - sheet_current_factor(9) = 1d5 - sheet_current_factor(10) = 0d0 - sheet_current_factor(11) = 0d0 / diff --git a/data/mephit_illustration.in b/data/mephit_illustration.in index 0867431..b7813e8 100644 --- a/data/mephit_illustration.in +++ b/data/mephit_illustration.in @@ -3,28 +3,24 @@ config%log_level = 4 ! 0: single iteration, 1: direct iterations, 2: preconditioned iterations config%runmode = 2 - ! use non-resonant test case - config%nonres = .false. - ! average over quadrilaterals for non-resonant test case - config%quad_avg = .true. - ! 0: EPS paper profile, 1: parabolic profile, 2: EQDSK profile + ! 0: EPS 2018 paper profile, 1: parabolic profile, 2: EQDSK profile config%pres_prof = 2 ! 0: only PS current, 1: rot B_0, 2: EQDSK profile config%curr_prof = 2 ! 0: surface integral, 1: field line integral, 2: EQDSK profile config%q_prof = 1 - ! 0: iMHD, 1: KiLCA interface + ! 0: geometric series around resonance, 1: Gaussian profile + config%refinement_scheme = 0 + ! 0: iMHD, 1: FLRE 2nd order config%currn_model = 1 ! 0: Viktor Nemov's code, 1: GPEC, 2: pre-computed Fourier transform config%vac_src = 2 - ! mass number of ions - config%m_i = 2.0 ! prefactor in Biot-Savart law, e.g. 5 windings and conversion from A to statA / c_0 config%Biot_Savart_prefactor = 5.0d-1 + !> file containing coil field data + config%coil_file = '/proj/plasma/CODE/lainer_p/MEPHIT/run/AUG_B_coils.h5' !> file containing coil currents config%currents_file = '/proj/plasma/DATA/BALANCE/COIL/33353/33353.2900_coil_markl.dat' - !> file containing coil geometries - config%coil_file = '/temp/lainer_p/MEPHIT-run/AUG_B_coils.h5' !> file containing equilibrium density profile config%dens_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_ne_PED_MMARKL_rho_pol.dat' !> file containing equilibrium electron temperature profile @@ -33,11 +29,13 @@ config%temp_i_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353.2900_Ti_PED_MMARKL_rho_pol.dat' !> file containing equilibrium radial electric field profile config%E_r_file = '/proj/plasma/DATA/BALANCE/PROF/33353/33353_2900_Er.dat' + ! mass number of ions + config%m_i = 2.0 ! toroidal mode number of perturbation config%n = 2 - ! maximum number of iteration steps + ! maximum number of (preconditioned) fixed-point iteration steps config%niter = 75 - ! maximum relative error in fixed-point iteration + ! maximum relative error in (preconditioned) fixed-point iteration config%iter_rel_err = 1d-9 ! maximum number of iterations in Arnoldi method config%nkrylov = 100 @@ -47,13 +45,22 @@ config%ritz_rel_err = 1d-9 ! maximum distance between flux surfaces along theta = 0 config%max_Delta_rad = 1.75d0 - ! enable damping of Pfirsch-Schlueter current + ! maximum number of points per flux surface (0: no maximum) + config%pol_max = 0 + ! minimum number of points per flux surface (0: no minimum) + config%pol_min = 0 + ! ignore resonance position where q = 1, which is usually spurious + config%ignore_q1_res = .true. + ! enable damping of Pfirsch-Schlueter current near resonances config%damp = .true. - ! number of points in resonance sweeep + ! maximal number of eigenvectors to be saved + config%max_eig_out = 10 + ! maximum poloidal mode number for Fourier transform of results + config%m_max = 24 + ! number of points in sweep over electron fluid resonance config%resonance_sweep = 0 - - config%debug_mfem = .false. - config%pol_max = 0 + ! offset added to radial electric field + config%offset_E_r = 0.0d0 / &arrays refinement(3) = 1.618d0 @@ -83,13 +90,4 @@ Delta_rad_res(9) = 5.347772320733589e-02 Delta_rad_res(10) = 6.066328245075899e-02 Delta_rad_res(11) = 3.464590203985631e-02 - sheet_current_factor(3) = 1d5 - sheet_current_factor(4) = 1d5 - sheet_current_factor(5) = 1d5 - sheet_current_factor(6) = 1d5 - sheet_current_factor(7) = 1d5 - sheet_current_factor(8) = 1d5 - sheet_current_factor(9) = 1d5 - sheet_current_factor(10) = 0d0 - sheet_current_factor(11) = 0d0 / diff --git a/src/mephit_conf.f90 b/src/mephit_conf.f90 index 0983c7e..6a18197 100644 --- a/src/mephit_conf.f90 +++ b/src/mephit_conf.f90 @@ -228,7 +228,7 @@ module mephit_conf integer :: m_min, m_max - !> Width of refined flux surfaces around resonances in cm. + !> Width of refined flux surfaces around resonances in cm. Defaults to 0. real(dp), dimension(:), allocatable :: Delta_rad_res !> Number of additional fine flux surfaces. Defaults to 0. @@ -237,7 +237,7 @@ module mephit_conf !> Width ratio of neighbouring refined flux surfaces. Defaults to 0 (no refinement). real(dp), dimension(:), allocatable :: refinement - !> Free parameters setting the magnitudes of sheet currents. Defaults to 0. + !> Free parameters setting the magnitudes of sheet currents. Only applied with refinement. Defaults to 1. real(dp), dimension(:), allocatable :: sheet_current_factor contains @@ -372,7 +372,7 @@ subroutine config_delayed_read(config, filename, m_min, m_max) Delta_rad_res = 0d0 add_fine = 0 refinement = 0d0 - sheet_current_factor = 0d0 + sheet_current_factor = 1d0 open(newunit = fid, file = filename) read(fid, nml = arrays) close(fid) diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index 2e94b2d..04c5849 100644 --- a/src/mephit_iter.F90 +++ b/src/mephit_iter.F90 @@ -1329,7 +1329,7 @@ subroutine compute_shielding_current(pn, jmnpar_over_Bmod) call L1_poloidal_modes(pn, pmn) do m = mesh%m_res_min, mesh%m_res_max m_res = -equil%cocos%sgn_q * m - if (abs(conf_arr%sheet_current_factor(m)) > 0d0) then + if (abs(conf_arr%refinement(m)) > 1d0) then kf = mesh%res_ind(m) jmnpar_over_Bmod%coeff(m_res, kf-1:kf) = & (pmn%coeff(m_res, kf) - pmn%coeff(m_res, kf-1)) * & From 64951ea9d651e1f86f4b601b2845977172b8aa99 Mon Sep 17 00:00:00 2001 From: Patrick Lainer Date: Fri, 25 Apr 2025 11:47:55 +0200 Subject: [PATCH 6/6] inline finite element projections --- src/mephit_iter.F90 | 3 +-- src/mephit_pert.f90 | 39 ++++++++++++++++++--------------------- 2 files changed, 19 insertions(+), 23 deletions(-) diff --git a/src/mephit_iter.F90 b/src/mephit_iter.F90 index 04c5849..d1e423f 100644 --- a/src/mephit_iter.F90 +++ b/src/mephit_iter.F90 @@ -1063,8 +1063,7 @@ subroutine compute_currn(perteq, fdm, flr2, apply_damping, debug_initial) use mephit_conf, only: conf, currn_model_kilca, currn_model_mhd, logger, datafile use mephit_mesh, only: mesh use mephit_pert, only: polmodes_t, polmodes_init, polmodes_write, polmodes_deinit, & - L1_t, L1_init, L1_interp, L1_deinit, RT0_t, RT0_init, RT0_deinit, & - RT0_interp, RT0_project_pol_comp, RT0_project_tor_comp + L1_t, L1_init, L1_deinit, RT0_t, RT0_init, RT0_deinit use mephit_flr2, only: flr2_t type(perteq_t), intent(inout) :: perteq type(fdm_t), intent(in) :: fdm diff --git a/src/mephit_pert.f90 b/src/mephit_pert.f90 index 5368ee4..474736f 100644 --- a/src/mephit_pert.f90 +++ b/src/mephit_pert.f90 @@ -336,14 +336,14 @@ subroutine RT0_project_pol_comp(elem, proj) type(RT0_t), intent(inout) :: elem procedure(vector_element_projection) :: proj integer :: kedge, k, ktri - real(dp) :: n_f(3) + real(dp) :: edge_perp(3) do kedge = 1, mesh%nedge ktri = mesh%edge_tri(1, kedge) - n_f = [mesh%edge_Z(kedge), 0d0, -mesh%edge_R(kedge)] + edge_perp = [mesh%edge_Z(kedge), 0d0, -mesh%edge_R(kedge)] do k = 1, mesh%GL_order elem%DOF(kedge) = elem%DOF(kedge) + proj(ktri, mesh%GL_weights(k) * mesh%GL_R(k, kedge), & - mesh%GL_R(k, kedge), mesh%GL_Z(k, kedge), n_f, cache%edge_fields(k, kedge)) + mesh%GL_R(k, kedge), mesh%GL_Z(k, kedge), edge_perp, cache%edge_fields(k, kedge)) end do end do end subroutine RT0_project_pol_comp @@ -353,12 +353,12 @@ subroutine RT0_project_tor_comp(elem, proj) type(RT0_t), intent(inout) :: elem procedure(vector_element_projection) :: proj integer :: ktri, k - real(dp), parameter :: n_f(3) = [0d0, 1d0, 0d0] + real(dp), parameter :: area_perp(3) = [0d0, 1d0, 0d0] do ktri = 1, mesh%ntri do k = 1, mesh%GL2_order elem%comp_phi(ktri) = elem%comp_phi(ktri) + proj(ktri, mesh%GL2_weights(k), & - mesh%GL2_R(k, ktri), mesh%GL2_Z(k, ktri), n_f, cache%area_fields(k, ktri)) + mesh%GL2_R(k, ktri), mesh%GL2_Z(k, ktri), area_perp, cache%area_fields(k, ktri)) end do end do end subroutine RT0_project_tor_comp @@ -1036,11 +1036,13 @@ end subroutine read_Bnvac_GPEC subroutine compute_Bnvac(Bn) use coil_tools, only: read_currents, read_Bnvac_Fourier use mephit_conf, only: conf, logger, vac_src_nemov, vac_src_gpec, vac_src_fourier + use mephit_mesh, only: mesh type(RT0_t), intent(inout) :: Bn - integer :: nR, nZ - real(dp) :: Rmin, Rmax, Zmin, Zmax + integer :: nR, nZ, kedge, k + real(dp) :: Rmin, Rmax, Zmin, Zmax, edge_perp(3) real(dp), dimension(:), allocatable :: Ic complex(dp), dimension(:, :), allocatable :: Bn_R, Bn_Z + complex(dp) :: Bn_splined(3) ! initialize vacuum field select case (conf%vac_src) @@ -1062,21 +1064,16 @@ subroutine compute_Bnvac(Bn) deallocate(Bn_R, Bn_Z) ! project to finite elements Bn%DOF = (0d0, 0d0) - call RT0_project_pol_comp(Bn, project_splined) + do kedge = 1, mesh%nedge + edge_perp = [mesh%edge_Z(kedge), 0d0, -mesh%edge_R(kedge)] + do k = 1, mesh%GL_order + call spline_bn(conf%n, mesh%GL_R(k, kedge), mesh%GL_Z(k, kedge), & + Bn_splined(1), Bn_splined(2), Bn_splined(3)) + Bn%DOF(kedge) = Bn%DOF(kedge) + mesh%GL_weights(k) * mesh%GL_R(k, kedge) * & + sum(Bn_splined * edge_perp) + end do + end do call RT0_tor_comp_from_zero_div(Bn) - - contains - function project_splined(ktri, weight, R, Z, n_f, f) - use mephit_mesh, only: field_cache_t - complex(dp) :: project_splined - integer, intent(in) :: ktri - real(dp), intent(in) :: weight, R, Z, n_f(:) - type(field_cache_t), intent(in) :: f - complex(dp) :: B(3) - - call spline_bn(conf%n, R, Z, B(1), B(2), B(3)) - project_splined = weight * sum(B * n_f) - end function project_splined end subroutine compute_Bnvac subroutine generate_vacfield(vac)