diff --git a/CMakeSources.in b/CMakeSources.in index 29030478..568c74ae 100644 --- a/CMakeSources.in +++ b/CMakeSources.in @@ -20,7 +20,13 @@ add_library(neo STATIC src/coordinates/vmec_coordinates.f90 src/coordinates/geoflux_coordinates.f90 src/coordinates/libneo_coordinate_conventions.f90 - src/coordinates/libneo_coordinates_merged.f90 + src/coordinates/libneo_coordinates_base.f90 + src/coordinates/libneo_coordinates_vmec.f90 + src/coordinates/libneo_coordinates_geoflux.f90 + src/coordinates/libneo_coordinates_chartmap.f90 + src/coordinates/libneo_coordinates_validator.f90 + src/coordinates/libneo_coordinates_file_detection.f90 + src/coordinates/libneo_coordinates.f90 src/coordinates/libneo_chartmap_vmec_generator.f90 src/coordinates/libneo_coordinates_mapping.f90 src/spline_vmec_data.f90 diff --git a/bind.toml b/bind.toml new file mode 100644 index 00000000..9f7572ac --- /dev/null +++ b/bind.toml @@ -0,0 +1,7 @@ +[python] +package-name = "pylibneo" +# Exclude files that cause issues with f90wrap +exclude = [ + "**/f2py_*.f90", + "**/local_rusage.c" +] diff --git a/fpm.toml b/fpm.toml index 3c8230f3..6689e292 100644 --- a/fpm.toml +++ b/fpm.toml @@ -6,13 +6,15 @@ maintainer = "albert@tugraz.at" copyright = "Copyright 2025, Plasma Physics at TU Graz" [dependencies] -openmp = "*" -mpi = "*" -hdf5 = "*" -netcdf = "*" +# System libraries (openmp, mpi, hdf5, netcdf) handled via compiler flags [build] -external-modules = ["fftw3.f03"] +external-modules = [ + "fftw3", "fftw3.f03", + "netcdf", "hdf5", "h5lt", + "omp_lib", "mpi", "mpi_f08", + "openacc" +] auto-executables = true auto-tests = true auto-examples = true diff --git a/src/coordinates/libneo_coordinates.f90 b/src/coordinates/libneo_coordinates.f90 index e240c66f..b5f42e71 100644 --- a/src/coordinates/libneo_coordinates.f90 +++ b/src/coordinates/libneo_coordinates.f90 @@ -1,265 +1,27 @@ module libneo_coordinates - use, intrinsic :: iso_fortran_env, only: dp => real64 - use interpolate, only: BatchSplineData3D - use libneo_coordinate_conventions + use libneo_coordinates_base + use libneo_coordinates_vmec + use libneo_coordinates_geoflux + use libneo_coordinates_chartmap + use libneo_coordinates_validator + use libneo_coordinates_file_detection implicit none - public :: UNKNOWN - public :: CYL, VMEC, BOOZER + public :: coordinate_system_t + public :: UNKNOWN, CYL, VMEC, BOOZER public :: RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM - character(len=*), parameter :: chartmap_netcdf_spec = & - "Chartmap NetCDF conventions:"//new_line('a')// & - "- Dimensions: rho, theta, zeta"//new_line('a')// & - "- Variables (required):"//new_line('a')// & - " rho(rho), theta(theta)"//new_line('a')// & - " zeta(zeta)"//new_line('a')// & - " x(zeta,theta,rho)"//new_line('a')// & - " y(zeta,theta,rho)"//new_line('a')// & - " z(zeta,theta,rho)"//new_line('a')// & - "- Optional: num_field_periods (integer >= 1)"// & - new_line('a')// & - "- Required: zeta_convention (global attribute)"// & - new_line('a')// & - " allowed: cyl, vmec"// & - new_line('a')// & - "- Ranges:"//new_line('a')// & - " rho in [0,1]"//new_line('a')// & - " theta in [0,2pi)"//new_line('a')// & - " zeta in [0,2pi/num_field_periods)"// & - new_line('a')// & - "periodic dims exclude endpoint"//new_line('a')// & - "- Storage order:"//new_line('a')// & - " file dims (zeta,theta,rho)"//new_line('a')// & - "Fortran reads x(rho,theta,zeta)"//new_line('a')// & - "- Units: x,y,z in cm" + public :: chartmap_netcdf_spec, chartmap_from_cyl_ierr_spec + public :: chartmap_from_cyl_ok, chartmap_from_cyl_err_max_iter + public :: chartmap_from_cyl_err_singular, chartmap_from_cyl_err_out_of_bounds + public :: chartmap_from_cyl_err_invalid + public :: refcoords_file_unknown, refcoords_file_chartmap, refcoords_file_vmec_wout - integer, parameter :: chartmap_from_cyl_ok = 0 - integer, parameter :: chartmap_from_cyl_err_max_iter = 1 - integer, parameter :: chartmap_from_cyl_err_singular = 2 - integer, parameter :: chartmap_from_cyl_err_out_of_bounds = 3 - integer, parameter :: chartmap_from_cyl_err_invalid = 4 + public :: vmec_coordinate_system_t, make_vmec_coordinate_system + public :: geoflux_coordinate_system_t, make_geoflux_coordinate_system + public :: chartmap_coordinate_system_t, make_chartmap_coordinate_system - integer, parameter :: refcoords_file_unknown = 0 - integer, parameter :: refcoords_file_chartmap = 1 - integer, parameter :: refcoords_file_vmec_wout = 2 - - character(len=*), parameter :: chartmap_from_cyl_ierr_spec = & - "chartmap_from_cyl ierr codes:"//new_line('a')// & - "0 ok"//new_line('a')// & - "1 max iterations or no progress"//new_line('a')// & - "2 singular normal equations"//new_line('a')// & - "3 step out of bounds for rho"//new_line('a')// & - "4 invalid mapping slice" - - type, abstract :: coordinate_system_t - contains - procedure(evaluate_cart_if), deferred :: evaluate_cart - procedure(evaluate_cyl_if), deferred :: evaluate_cyl - procedure(covariant_basis_if), deferred :: covariant_basis - procedure(metric_tensor_if), deferred :: metric_tensor - procedure(from_cyl_if), deferred :: from_cyl - procedure :: cov_to_cart => coordinate_system_cov_to_cart - procedure :: ctr_to_cart => coordinate_system_ctr_to_cart - end type coordinate_system_t - - abstract interface - subroutine evaluate_cart_if(self, u, x) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine - - subroutine evaluate_cyl_if(self, u, x) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine - - subroutine covariant_basis_if(self, u, e_cov) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3, 3) - end subroutine - - subroutine metric_tensor_if(self, u, g, ginv, sqrtg) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - end subroutine - - subroutine from_cyl_if(self, xcyl, u, ierr) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - end subroutine - end interface - - interface - module subroutine make_vmec_coordinate_system(cs) - class(coordinate_system_t), allocatable, intent(out) :: cs - end subroutine - - module subroutine make_geoflux_coordinate_system(cs) - class(coordinate_system_t), allocatable, intent(out) :: cs - end subroutine - - module subroutine make_chartmap_coordinate_system(cs, filename) - class(coordinate_system_t), allocatable, intent(out) :: cs - character(len=*), intent(in) :: filename - end subroutine - - module subroutine validate_chartmap_file(filename, ierr, message) - character(len=*), intent(in) :: filename - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - end subroutine validate_chartmap_file - - module subroutine detect_refcoords_file_type(filename, file_type, ierr, message) - character(len=*), intent(in) :: filename - integer, intent(out) :: file_type - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - end subroutine detect_refcoords_file_type - end interface - - type, extends(coordinate_system_t) :: vmec_coordinate_system_t - contains - procedure :: evaluate_cart => vmec_evaluate_cart - procedure :: evaluate_cyl => vmec_evaluate_cyl - procedure :: covariant_basis => vmec_covariant_basis - procedure :: metric_tensor => vmec_metric_tensor - procedure :: from_cyl => vmec_from_cyl - end type vmec_coordinate_system_t - - type, extends(coordinate_system_t) :: chartmap_coordinate_system_t - type(BatchSplineData3D) :: spl_cart - type(BatchSplineData3D) :: spl_rz - logical :: has_spl_rz = .false. - integer :: nrho = 0 - integer :: ntheta = 0 - integer :: nzeta = 0 - integer :: num_field_periods = 1 - integer :: zeta_convention = UNKNOWN - integer :: rho_convention = UNKNOWN - real(dp) :: tol_newton = 1.0e-12_dp - contains - procedure :: evaluate_cart => chartmap_evaluate_cart - procedure :: evaluate_cyl => chartmap_evaluate_cyl - procedure :: covariant_basis => chartmap_covariant_basis - procedure :: metric_tensor => chartmap_metric_tensor - procedure :: from_cyl => chartmap_from_cyl - procedure :: from_cart => chartmap_from_cart - end type chartmap_coordinate_system_t - - interface - module subroutine vmec_evaluate_cart(self, u, x) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine vmec_evaluate_cart - - module subroutine vmec_evaluate_cyl(self, u, x) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine vmec_evaluate_cyl - - module subroutine vmec_covariant_basis(self, u, e_cov) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3, 3) - end subroutine vmec_covariant_basis - - module subroutine vmec_metric_tensor(self, u, g, ginv, sqrtg) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - end subroutine vmec_metric_tensor - - module subroutine vmec_from_cyl(self, xcyl, u, ierr) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - end subroutine vmec_from_cyl - - module subroutine chartmap_evaluate_cart(self, u, x) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine chartmap_evaluate_cart - - module subroutine chartmap_evaluate_cyl(self, u, x) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine chartmap_evaluate_cyl - - module subroutine chartmap_covariant_basis(self, u, e_cov) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3, 3) - end subroutine chartmap_covariant_basis - - module subroutine chartmap_metric_tensor(self, u, g, ginv, sqrtg) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - end subroutine chartmap_metric_tensor - - module subroutine chartmap_from_cyl(self, xcyl, u, ierr) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - end subroutine chartmap_from_cyl - - module subroutine chartmap_from_cart(self, x, u, ierr) - class(chartmap_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: x(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - end subroutine chartmap_from_cart - end interface - -contains - - subroutine coordinate_system_cov_to_cart(self, u, v_cov, v_cart) - !> Convert covariant vector components to Cartesian. - !> v_cart_i = e_cov(i,j) * ginv(j,k) * v_cov(k) - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(in) :: v_cov(3) - real(dp), intent(out) :: v_cart(3) - - real(dp) :: e_cov(3, 3), g(3, 3), ginv(3, 3), sqrtg - real(dp) :: v_ctr(3) - - call self%metric_tensor(u, g, ginv, sqrtg) - call self%covariant_basis(u, e_cov) - v_ctr = matmul(ginv, v_cov) - v_cart = matmul(e_cov, v_ctr) - end subroutine coordinate_system_cov_to_cart - - subroutine coordinate_system_ctr_to_cart(self, u, v_ctr, v_cart) - !> Convert contravariant vector components to Cartesian. - !> v_cart_i = e_cov(i,j) * v_ctr(j) - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(in) :: v_ctr(3) - real(dp), intent(out) :: v_cart(3) - - real(dp) :: e_cov(3, 3) - - call self%covariant_basis(u, e_cov) - v_cart = matmul(e_cov, v_ctr) - end subroutine coordinate_system_ctr_to_cart + public :: validate_chartmap_file + public :: detect_refcoords_file_type end module libneo_coordinates diff --git a/src/coordinates/libneo_coordinates_base.f90 b/src/coordinates/libneo_coordinates_base.f90 new file mode 100644 index 00000000..9d6bb1ea --- /dev/null +++ b/src/coordinates/libneo_coordinates_base.f90 @@ -0,0 +1,135 @@ +module libneo_coordinates_base + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinate_conventions + implicit none + private + + public :: coordinate_system_t + public :: UNKNOWN, CYL, VMEC, BOOZER + public :: RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM + + public :: chartmap_netcdf_spec, chartmap_from_cyl_ierr_spec + public :: chartmap_from_cyl_ok, chartmap_from_cyl_err_max_iter + public :: chartmap_from_cyl_err_singular, chartmap_from_cyl_err_out_of_bounds + public :: chartmap_from_cyl_err_invalid + public :: refcoords_file_unknown, refcoords_file_chartmap, refcoords_file_vmec_wout + + character(len=*), parameter :: chartmap_netcdf_spec = & + "Chartmap NetCDF conventions:"//new_line('a')// & + "- Dimensions: rho, theta, zeta"//new_line('a')// & + "- Variables (required):"//new_line('a')// & + " rho(rho), theta(theta)"//new_line('a')// & + " zeta(zeta)"//new_line('a')// & + " x(zeta,theta,rho)"//new_line('a')// & + " y(zeta,theta,rho)"//new_line('a')// & + " z(zeta,theta,rho)"//new_line('a')// & + "- Optional: num_field_periods (integer >= 1)"//new_line('a')// & + "- Required: zeta_convention (global attribute)"//new_line('a')// & + " allowed: cyl, vmec"//new_line('a')// & + "- Ranges:"//new_line('a')// & + " rho in [0,1]"//new_line('a')// & + " theta in [0,2pi)"//new_line('a')// & + " zeta in [0,2pi/num_field_periods)"//new_line('a')// & + "periodic dims exclude endpoint"//new_line('a')// & + "- Storage order:"//new_line('a')// & + " file dims (zeta,theta,rho)"//new_line('a')// & + "Fortran reads x(rho,theta,zeta)"//new_line('a')// & + "- Units: x,y,z in cm" + + integer, parameter :: chartmap_from_cyl_ok = 0 + integer, parameter :: chartmap_from_cyl_err_max_iter = 1 + integer, parameter :: chartmap_from_cyl_err_singular = 2 + integer, parameter :: chartmap_from_cyl_err_out_of_bounds = 3 + integer, parameter :: chartmap_from_cyl_err_invalid = 4 + + integer, parameter :: refcoords_file_unknown = 0 + integer, parameter :: refcoords_file_chartmap = 1 + integer, parameter :: refcoords_file_vmec_wout = 2 + + character(len=*), parameter :: chartmap_from_cyl_ierr_spec = & + "chartmap_from_cyl ierr codes:"//new_line('a')// & + "0 ok"//new_line('a')// & + "1 max iterations or no progress"//new_line('a')// & + "2 singular normal equations"//new_line('a')// & + "3 step out of bounds for rho"//new_line('a')// & + "4 invalid mapping slice" + + type, abstract :: coordinate_system_t + contains + procedure(evaluate_cart_if), deferred :: evaluate_cart + procedure(evaluate_cyl_if), deferred :: evaluate_cyl + procedure(covariant_basis_if), deferred :: covariant_basis + procedure(metric_tensor_if), deferred :: metric_tensor + procedure(from_cyl_if), deferred :: from_cyl + procedure :: cov_to_cart => coordinate_system_cov_to_cart + procedure :: ctr_to_cart => coordinate_system_ctr_to_cart + end type coordinate_system_t + + abstract interface + subroutine evaluate_cart_if(self, u, x) + import :: coordinate_system_t, dp + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(out) :: x(3) + end subroutine + + subroutine evaluate_cyl_if(self, u, x) + import :: coordinate_system_t, dp + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(out) :: x(3) + end subroutine + + subroutine covariant_basis_if(self, u, e_cov) + import :: coordinate_system_t, dp + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(out) :: e_cov(3, 3) + end subroutine + + subroutine metric_tensor_if(self, u, g, ginv, sqrtg) + import :: coordinate_system_t, dp + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg + end subroutine + + subroutine from_cyl_if(self, xcyl, u, ierr) + import :: coordinate_system_t, dp + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: xcyl(3) + real(dp), intent(out) :: u(3) + integer, intent(out) :: ierr + end subroutine + end interface + +contains + + subroutine coordinate_system_cov_to_cart(self, u, v_cov, v_cart) + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(in) :: v_cov(3) + real(dp), intent(out) :: v_cart(3) + + real(dp) :: e_cov(3, 3), g(3, 3), ginv(3, 3), sqrtg + real(dp) :: v_ctr(3) + + call self%metric_tensor(u, g, ginv, sqrtg) + v_ctr = matmul(ginv, v_cov) + call self%covariant_basis(u, e_cov) + v_cart = matmul(e_cov, v_ctr) + end subroutine coordinate_system_cov_to_cart + + subroutine coordinate_system_ctr_to_cart(self, u, v_ctr, v_cart) + class(coordinate_system_t), intent(in) :: self + real(dp), intent(in) :: u(3) + real(dp), intent(in) :: v_ctr(3) + real(dp), intent(out) :: v_cart(3) + + real(dp) :: e_cov(3, 3) + + call self%covariant_basis(u, e_cov) + v_cart = matmul(e_cov, v_ctr) + end subroutine coordinate_system_ctr_to_cart + +end module libneo_coordinates_base diff --git a/src/coordinates/libneo_coordinates_chartmap.f90 b/src/coordinates/libneo_coordinates_chartmap.f90 index c85ca67f..d329f74e 100644 --- a/src/coordinates/libneo_coordinates_chartmap.f90 +++ b/src/coordinates/libneo_coordinates_chartmap.f90 @@ -1,11 +1,43 @@ -submodule(libneo_coordinates) libneo_coordinates_chartmap +module libneo_coordinates_chartmap + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates_base, only: coordinate_system_t, & + UNKNOWN, CYL, VMEC, BOOZER, & + RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM, & + chartmap_from_cyl_ok, chartmap_from_cyl_err_max_iter, & + chartmap_from_cyl_err_singular, & + chartmap_from_cyl_err_out_of_bounds, & + chartmap_from_cyl_err_invalid + use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & + evaluate_batch_splines_3d, evaluate_batch_splines_3d_der use nctools_module, only: nc_open, nc_close, nc_inq_dim, nc_get use math_constants, only: TWOPI - use interpolate, only: construct_batch_splines_3d, & - evaluate_batch_splines_3d, evaluate_batch_splines_3d_der use netcdf, only: NF90_CHAR, NF90_GLOBAL, NF90_NOERR, nf90_get_att, nf90_get_var, & nf90_inq_varid, nf90_inquire_attribute implicit none + private + + public :: chartmap_coordinate_system_t + public :: make_chartmap_coordinate_system + + type, extends(coordinate_system_t) :: chartmap_coordinate_system_t + type(BatchSplineData3D) :: spl_cart + type(BatchSplineData3D) :: spl_rz + logical :: has_spl_rz = .false. + integer :: nrho = 0 + integer :: ntheta = 0 + integer :: nzeta = 0 + integer :: num_field_periods = 1 + integer :: zeta_convention = UNKNOWN + integer :: rho_convention = UNKNOWN + real(dp) :: tol_newton = 1.0e-12_dp + contains + procedure :: evaluate_cart => chartmap_evaluate_cart + procedure :: evaluate_cyl => chartmap_evaluate_cyl + procedure :: covariant_basis => chartmap_covariant_basis + procedure :: metric_tensor => chartmap_metric_tensor + procedure :: from_cyl => chartmap_from_cyl + procedure :: from_cart => chartmap_from_cart + end type chartmap_coordinate_system_t contains @@ -287,30 +319,24 @@ subroutine chartmap_eval_cart_der(self, u, vals, dvals) vals(2) = rz(1)*sph vals(3) = rz(2) - ! dvals(cart_component, coord_index) = d(x_cart)/d(u_coord) - ! drz(derivative_dim, quantity) where quantity 1=R, 2=Z - ! X = R*cos(phi), Y = R*sin(phi), Z = Z + dvals(1, 1) = drz(1, 1)*cph + dvals(2, 1) = drz(1, 1)*sph + dvals(3, 1) = drz(1, 2) - ! d/drho derivatives: - dvals(1, 1) = drz(1, 1)*cph ! dX/drho = dR/drho * cos - dvals(2, 1) = drz(1, 1)*sph ! dY/drho = dR/drho * sin - dvals(3, 1) = drz(1, 2) ! dZ/drho = dZ/drho + dvals(1, 2) = drz(2, 1)*cph + dvals(2, 2) = drz(2, 1)*sph + dvals(3, 2) = drz(2, 2) - ! d/dtheta derivatives: - dvals(1, 2) = drz(2, 1)*cph ! dX/dtheta = dR/dtheta * cos - dvals(2, 2) = drz(2, 1)*sph ! dY/dtheta = dR/dtheta * sin - dvals(3, 2) = drz(2, 2) ! dZ/dtheta = dZ/dtheta - - ! d/dphi derivatives: - dvals(1, 3) = drz(3, 1)*cph - rz(1)*sph ! dX/dphi = dR/dphi*cos - R*sin - dvals(2, 3) = drz(3, 1)*sph + rz(1)*cph ! dY/dphi = dR/dphi*sin + R*cos - dvals(3, 3) = drz(3, 2) ! dZ/dphi = dZ/dphi + dvals(1, 3) = drz(3, 1)*cph - rz(1)*sph + dvals(2, 3) = drz(3, 1)*sph + rz(1)*cph + dvals(3, 3) = drz(3, 2) else call evaluate_batch_splines_3d_der(self%spl_cart, u_eval, vals, dvals) end if end subroutine chartmap_eval_cart_der - module subroutine make_chartmap_coordinate_system(cs, filename) + subroutine make_chartmap_coordinate_system(cs, filename) + use libneo_coordinates_validator, only: validate_chartmap_file class(coordinate_system_t), allocatable, intent(out) :: cs character(len=*), intent(in) :: filename @@ -325,6 +351,7 @@ module subroutine make_chartmap_coordinate_system(cs, filename) end subroutine make_chartmap_coordinate_system subroutine initialize_chartmap(ccs, filename) + use libneo_coordinates_validator, only: validate_chartmap_file type(chartmap_coordinate_system_t), intent(inout) :: ccs character(len=*), intent(in) :: filename @@ -522,9 +549,6 @@ subroutine chartmap_covariant_basis(self, u, e_cov) real(dp) :: dvals(3, 3) call chartmap_eval_cart_der(self, u, vals, dvals) - - ! e_cov(cart_component, coord_index) = d(x_cart)/d(u_coord) - ! dvals has the same indexing: dvals(cart_component, coord_index) e_cov = dvals end subroutine chartmap_covariant_basis @@ -788,8 +812,6 @@ subroutine chartmap_eval_residual_and_partials_cart(self, x_target, rho, & residual = x_target - vals res_norm = sqrt(sum(residual**2)) - ! dvals(cart_component, coord_index) = d(x_cart)/d(u_coord) - ! Extract columns: dx_drho = [dX/drho, dY/drho, dZ/drho] = dvals(:, 1) dx_drho = dvals(:, 1) dx_dtheta = dvals(:, 2) dx_dzeta = dvals(:, 3) @@ -1063,7 +1085,6 @@ subroutine chartmap_eval_residual_and_partials(self, x_target, rho, theta, zeta, residual = x_target - vals res_norm = sqrt(sum(residual**2)) - ! dvals(cart_component, coord_index) = d(x_cart)/d(u_coord) dx_drho = dvals(:, 1) dx_dtheta = dvals(:, 2) end subroutine chartmap_eval_residual_and_partials @@ -1188,4 +1209,4 @@ subroutine chartmap_solve_slice(self, x_target, zeta, rho, theta, ierr, trace) end if end subroutine chartmap_solve_slice -end submodule libneo_coordinates_chartmap +end module libneo_coordinates_chartmap diff --git a/src/coordinates/libneo_coordinates_file_detection.f90 b/src/coordinates/libneo_coordinates_file_detection.f90 index bca5a354..e780515f 100644 --- a/src/coordinates/libneo_coordinates_file_detection.f90 +++ b/src/coordinates/libneo_coordinates_file_detection.f90 @@ -1,11 +1,18 @@ -submodule(libneo_coordinates) libneo_coordinates_file_detection +module libneo_coordinates_file_detection + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates_base, only: refcoords_file_unknown, refcoords_file_chartmap, & + refcoords_file_vmec_wout + use libneo_coordinates_validator, only: validate_chartmap_file use netcdf, only: NF90_NOERR, NF90_NOWRITE, nf90_close, nf90_inq_dimid, & nf90_inq_varid, nf90_open, nf90_strerror implicit none + private + + public :: detect_refcoords_file_type contains - module subroutine detect_refcoords_file_type(filename, file_type, ierr, message) + subroutine detect_refcoords_file_type(filename, file_type, ierr, message) character(len=*), intent(in) :: filename integer, intent(out) :: file_type integer, intent(out) :: ierr @@ -56,4 +63,4 @@ module subroutine detect_refcoords_file_type(filename, file_type, ierr, message) status = nf90_close(ncid) end subroutine detect_refcoords_file_type -end submodule libneo_coordinates_file_detection +end module libneo_coordinates_file_detection diff --git a/src/coordinates/libneo_coordinates_geoflux.f90 b/src/coordinates/libneo_coordinates_geoflux.f90 index ca1b9c6d..1250cade 100644 --- a/src/coordinates/libneo_coordinates_geoflux.f90 +++ b/src/coordinates/libneo_coordinates_geoflux.f90 @@ -1,7 +1,12 @@ -submodule (libneo_coordinates) libneo_coordinates_geoflux +module libneo_coordinates_geoflux + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates_base, only: coordinate_system_t use geoflux_coordinates, only: geoflux_to_cyl, assign_geoflux_to_cyl_jacobian use cylindrical_cartesian, only: cyl_to_cart implicit none + private + + public :: geoflux_coordinate_system_t, make_geoflux_coordinate_system type, extends(coordinate_system_t) :: geoflux_coordinate_system_t contains @@ -14,7 +19,7 @@ contains - module subroutine make_geoflux_coordinate_system(cs) + subroutine make_geoflux_coordinate_system(cs) class(coordinate_system_t), allocatable, intent(out) :: cs allocate(geoflux_coordinate_system_t :: cs) end subroutine make_geoflux_coordinate_system @@ -24,6 +29,9 @@ subroutine geoflux_evaluate_cyl(self, u, x) real(dp), intent(in) :: u(3) real(dp), intent(out) :: x(3) + associate(dummy => self) + end associate + call geoflux_to_cyl(u, x) end subroutine geoflux_evaluate_cyl @@ -41,14 +49,17 @@ end subroutine geoflux_evaluate_cart subroutine geoflux_covariant_basis(self, u, e_cov) class(geoflux_coordinate_system_t), intent(in) :: self real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3,3) + real(dp), intent(out) :: e_cov(3, 3) real(dp) :: s, theta, phi real(dp) :: R, Z - real(dp) :: jac_cyl(3,3) + real(dp) :: jac_cyl(3, 3) real(dp) :: cos_phi, sin_phi real(dp) :: xgeo(3), xcyl(3) + associate(dummy => self) + end associate + s = u(1) theta = u(2) phi = u(3) @@ -63,25 +74,25 @@ subroutine geoflux_covariant_basis(self, u, e_cov) cos_phi = cos(phi) sin_phi = sin(phi) - e_cov(1, 1) = jac_cyl(1,1) * cos_phi - e_cov(2, 1) = jac_cyl(1,1) * sin_phi - e_cov(3, 1) = jac_cyl(3,1) + e_cov(1, 1) = jac_cyl(1, 1)*cos_phi + e_cov(2, 1) = jac_cyl(1, 1)*sin_phi + e_cov(3, 1) = jac_cyl(3, 1) - e_cov(1, 2) = jac_cyl(1,2) * cos_phi - e_cov(2, 2) = jac_cyl(1,2) * sin_phi - e_cov(3, 2) = jac_cyl(3,2) + e_cov(1, 2) = jac_cyl(1, 2)*cos_phi + e_cov(2, 2) = jac_cyl(1, 2)*sin_phi + e_cov(3, 2) = jac_cyl(3, 2) - e_cov(1, 3) = jac_cyl(1,3) * cos_phi - R * sin_phi - e_cov(2, 3) = jac_cyl(1,3) * sin_phi + R * cos_phi - e_cov(3, 3) = jac_cyl(3,3) + e_cov(1, 3) = jac_cyl(1, 3)*cos_phi - R*sin_phi + e_cov(2, 3) = jac_cyl(1, 3)*sin_phi + R*cos_phi + e_cov(3, 3) = jac_cyl(3, 3) end subroutine geoflux_covariant_basis subroutine geoflux_metric_tensor(self, u, g, ginv, sqrtg) class(geoflux_coordinate_system_t), intent(in) :: self real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg + real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - real(dp) :: e_cov(3,3) + real(dp) :: e_cov(3, 3) real(dp) :: det integer :: i, j @@ -116,7 +127,14 @@ subroutine geoflux_from_cyl(self, xcyl, u, ierr) real(dp), intent(out) :: u(3) integer, intent(out) :: ierr + associate(dummy => self) + end associate + associate(dummy2 => xcyl) + end associate + + u = 0.0_dp + ierr = 0 error stop "geoflux_from_cyl: not implemented" end subroutine geoflux_from_cyl -end submodule libneo_coordinates_geoflux +end module libneo_coordinates_geoflux diff --git a/src/coordinates/libneo_coordinates_merged.f90 b/src/coordinates/libneo_coordinates_merged.f90 deleted file mode 100644 index bd5194dc..00000000 --- a/src/coordinates/libneo_coordinates_merged.f90 +++ /dev/null @@ -1,1171 +0,0 @@ -module libneo_coordinates - use, intrinsic :: iso_fortran_env, only: dp => real64 - use interpolate, only: BatchSplineData3D, construct_batch_splines_3d, & - evaluate_batch_splines_3d, evaluate_batch_splines_3d_der - use libneo_coordinate_conventions - use nctools_module, only: nc_open, nc_close, nc_inq_dim, nc_get - use math_constants, only: TWOPI - use netcdf, only: NF90_BYTE, NF90_CHAR, NF90_DOUBLE, NF90_GLOBAL, NF90_INT, & - NF90_INT64, NF90_MAX_VAR_DIMS, NF90_NOERR, NF90_NOWRITE, & - NF90_SHORT, & - nf90_close, nf90_get_att, nf90_get_var, nf90_inq_dimid, & - nf90_inq_varid, nf90_inquire_attribute, & - nf90_inquire_dimension, nf90_inquire_variable, & - nf90_open, nf90_strerror - use spline_vmec_sub, only: splint_vmec_data - use cylindrical_cartesian, only: cyl_to_cart, cart_to_cyl - use geoflux_coordinates, only: geoflux_to_cyl, assign_geoflux_to_cyl_jacobian - implicit none - - public :: UNKNOWN - public :: CYL, VMEC, BOOZER - public :: RHO_TOR, RHO_POL, PSI_TOR_NORM, PSI_POL_NORM - - character(len=*), parameter :: chartmap_netcdf_spec = & - "Chartmap NetCDF conventions:"//new_line('a')// & - "- Dimensions: rho, theta, zeta"//new_line('a')// & - "- Variables (required):"//new_line('a')// & - " rho(rho), theta(theta)"//new_line('a')// & - " zeta(zeta)"//new_line('a')// & - " x(zeta,theta,rho)"//new_line('a')// & - " y(zeta,theta,rho)"//new_line('a')// & - " z(zeta,theta,rho)"//new_line('a')// & - "- Optional: num_field_periods (integer >= 1)"// & - new_line('a')// & - "- Required: zeta_convention (global attribute)"// & - new_line('a')// & - " allowed: cyl, vmec"// & - new_line('a')// & - "- Ranges:"//new_line('a')// & - " rho in [0,1]"//new_line('a')// & - " theta in [0,2pi)"//new_line('a')// & - " zeta in [0,2pi/num_field_periods)"// & - new_line('a')// & - "periodic dims exclude endpoint"//new_line('a')// & - "- Storage order:"//new_line('a')// & - " file dims (zeta,theta,rho)"//new_line('a')// & - "Fortran reads x(rho,theta,zeta)"//new_line('a')// & - "- Units: x,y,z in cm" - - integer, parameter :: chartmap_from_cyl_ok = 0 - integer, parameter :: chartmap_from_cyl_err_max_iter = 1 - integer, parameter :: chartmap_from_cyl_err_singular = 2 - integer, parameter :: chartmap_from_cyl_err_out_of_bounds = 3 - integer, parameter :: chartmap_from_cyl_err_invalid = 4 - - integer, parameter :: refcoords_file_unknown = 0 - integer, parameter :: refcoords_file_chartmap = 1 - integer, parameter :: refcoords_file_vmec_wout = 2 - - character(len=*), parameter :: chartmap_from_cyl_ierr_spec = & - "chartmap_from_cyl ierr codes:"//new_line('a')// & - "0 ok"//new_line('a')// & - "1 max iterations or no progress"//new_line('a')// & - "2 singular normal equations"//new_line('a')// & - "3 step out of bounds for rho"//new_line('a')// & - "4 invalid mapping slice" - - type, abstract :: coordinate_system_t - contains - procedure(evaluate_cart_if), deferred :: evaluate_cart - procedure(evaluate_cyl_if), deferred :: evaluate_cyl - procedure(covariant_basis_if), deferred :: covariant_basis - procedure(metric_tensor_if), deferred :: metric_tensor - procedure(from_cyl_if), deferred :: from_cyl - procedure :: cov_to_cart => coordinate_system_cov_to_cart - procedure :: ctr_to_cart => coordinate_system_ctr_to_cart - end type coordinate_system_t - - abstract interface - subroutine evaluate_cart_if(self, u, x) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine - - subroutine evaluate_cyl_if(self, u, x) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - end subroutine - - subroutine covariant_basis_if(self, u, e_cov) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3, 3) - end subroutine - - subroutine metric_tensor_if(self, u, g, ginv, sqrtg) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - end subroutine - - subroutine from_cyl_if(self, xcyl, u, ierr) - import :: coordinate_system_t, dp - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - end subroutine - end interface - - type, extends(coordinate_system_t) :: vmec_coordinate_system_t - contains - procedure :: evaluate_cart => vmec_evaluate_cart - procedure :: evaluate_cyl => vmec_evaluate_cyl - procedure :: covariant_basis => vmec_covariant_basis - procedure :: metric_tensor => vmec_metric_tensor - procedure :: from_cyl => vmec_from_cyl - end type vmec_coordinate_system_t - - type, extends(coordinate_system_t) :: chartmap_coordinate_system_t - type(BatchSplineData3D) :: spl_cart - type(BatchSplineData3D) :: spl_rz - logical :: has_spl_rz = .false. - integer :: nrho = 0 - integer :: ntheta = 0 - integer :: nzeta = 0 - integer :: num_field_periods = 1 - integer :: zeta_convention = UNKNOWN - integer :: rho_convention = UNKNOWN - real(dp) :: tol_newton = 1.0e-12_dp - contains - procedure :: evaluate_cart => chartmap_evaluate_cart - procedure :: evaluate_cyl => chartmap_evaluate_cyl - procedure :: covariant_basis => chartmap_covariant_basis - procedure :: metric_tensor => chartmap_metric_tensor - procedure :: from_cyl => chartmap_from_cyl - procedure :: from_cart => chartmap_from_cart - end type chartmap_coordinate_system_t - - type, extends(coordinate_system_t) :: geoflux_coordinate_system_t - contains - procedure :: evaluate_cart => geoflux_evaluate_cart - procedure :: evaluate_cyl => geoflux_evaluate_cyl - procedure :: covariant_basis => geoflux_covariant_basis - procedure :: metric_tensor => geoflux_metric_tensor - procedure :: from_cyl => geoflux_from_cyl - end type geoflux_coordinate_system_t - - public :: make_vmec_coordinate_system - public :: make_chartmap_coordinate_system - public :: make_geoflux_coordinate_system - public :: validate_chartmap_file - public :: detect_refcoords_file_type - -contains - - subroutine coordinate_system_cov_to_cart(self, u, v_cov, v_cart) - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(in) :: v_cov(3) - real(dp), intent(out) :: v_cart(3) - - real(dp) :: e_cov(3, 3), g(3, 3), ginv(3, 3), sqrtg - real(dp) :: v_ctr(3) - - call self%metric_tensor(u, g, ginv, sqrtg) - call self%covariant_basis(u, e_cov) - v_ctr = matmul(ginv, v_cov) - v_cart = matmul(e_cov, v_ctr) - end subroutine coordinate_system_cov_to_cart - - subroutine coordinate_system_ctr_to_cart(self, u, v_ctr, v_cart) - class(coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(in) :: v_ctr(3) - real(dp), intent(out) :: v_cart(3) - - real(dp) :: e_cov(3, 3) - - call self%covariant_basis(u, e_cov) - v_cart = matmul(e_cov, v_ctr) - end subroutine coordinate_system_ctr_to_cart - - ! ======================================================================== - ! VMEC coordinate system procedures - ! ======================================================================== - - subroutine make_vmec_coordinate_system(cs) - class(coordinate_system_t), allocatable, intent(out) :: cs - allocate(vmec_coordinate_system_t :: cs) - end subroutine make_vmec_coordinate_system - - subroutine vmec_evaluate_cyl(self, u, x) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - - real(dp) :: s, theta, varphi - real(dp) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota - real(dp) :: R, Z, alam, dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp - real(dp) :: dl_ds, dl_dt, dl_dp - - associate(dummy => self) - end associate - - s = u(1) - theta = u(2) - varphi = u(3) - - call splint_vmec_data(s, theta, varphi, A_phi, A_theta, dA_phi_ds, & - dA_theta_ds, aiota, R, Z, alam, dR_ds, dR_dt, & - dR_dp, dZ_ds, dZ_dt, dZ_dp, dl_ds, dl_dt, dl_dp) - - x(1) = R - x(2) = varphi - x(3) = Z - end subroutine vmec_evaluate_cyl - - subroutine vmec_evaluate_cart(self, u, x) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - - real(dp) :: xcyl(3) - - call self%evaluate_cyl(u, xcyl) - call cyl_to_cart(xcyl, x) - end subroutine vmec_evaluate_cart - - subroutine vmec_covariant_basis(self, u, e_cov) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3,3) - - real(dp) :: s, theta, varphi - real(dp) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota - real(dp) :: R, Z, alam, dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp - real(dp) :: dl_ds, dl_dt, dl_dp - real(dp) :: cos_phi, sin_phi - - associate(dummy => self) - end associate - - s = u(1) - theta = u(2) - varphi = u(3) - - call splint_vmec_data(s, theta, varphi, A_phi, A_theta, dA_phi_ds, & - dA_theta_ds, aiota, R, Z, alam, dR_ds, dR_dt, & - dR_dp, dZ_ds, dZ_dt, dZ_dp, dl_ds, dl_dt, dl_dp) - - cos_phi = cos(varphi) - sin_phi = sin(varphi) - - e_cov(1, 1) = dR_ds*cos_phi - e_cov(2, 1) = dR_ds*sin_phi - e_cov(3, 1) = dZ_ds - - e_cov(1, 2) = dR_dt*cos_phi - e_cov(2, 2) = dR_dt*sin_phi - e_cov(3, 2) = dZ_dt - - e_cov(1, 3) = dR_dp*cos_phi - R*sin_phi - e_cov(2, 3) = dR_dp*sin_phi + R*cos_phi - e_cov(3, 3) = dZ_dp - end subroutine vmec_covariant_basis - - subroutine vmec_metric_tensor(self, u, g, ginv, sqrtg) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg - - real(dp) :: e_cov(3,3) - real(dp) :: det - integer :: i, j - - call self%covariant_basis(u, e_cov) - - do i = 1, 3 - do j = 1, 3 - g(i, j) = dot_product(e_cov(:, i), e_cov(:, j)) - end do - end do - - det = g(1,1)*(g(2,2)*g(3,3) - g(2,3)*g(3,2)) & - - g(1,2)*(g(2,1)*g(3,3) - g(2,3)*g(3,1)) & - + g(1,3)*(g(2,1)*g(3,2) - g(2,2)*g(3,1)) - - sqrtg = sqrt(abs(det)) - - ginv(1,1) = (g(2,2)*g(3,3) - g(2,3)*g(3,2))/det - ginv(1,2) = (g(1,3)*g(3,2) - g(1,2)*g(3,3))/det - ginv(1,3) = (g(1,2)*g(2,3) - g(1,3)*g(2,2))/det - ginv(2,1) = (g(2,3)*g(3,1) - g(2,1)*g(3,3))/det - ginv(2,2) = (g(1,1)*g(3,3) - g(1,3)*g(3,1))/det - ginv(2,3) = (g(1,3)*g(2,1) - g(1,1)*g(2,3))/det - ginv(3,1) = (g(2,1)*g(3,2) - g(2,2)*g(3,1))/det - ginv(3,2) = (g(1,2)*g(3,1) - g(1,1)*g(3,2))/det - ginv(3,3) = (g(1,1)*g(2,2) - g(1,2)*g(2,1))/det - end subroutine vmec_metric_tensor - - subroutine vmec_from_cyl(self, xcyl, u, ierr) - class(vmec_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - - integer, parameter :: max_iter = 50 - real(dp), parameter :: tol_res = 1.0e-10_dp - real(dp), parameter :: tol_step = 1.0e-12_dp - - real(dp) :: s, theta, varphi - real(dp) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota - real(dp) :: R, Z, alam - real(dp) :: dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp - real(dp) :: dl_ds, dl_dt, dl_dp - real(dp) :: R_axis, Z_axis, R_bnd, Z_bnd - real(dp) :: rs, zs, dist_axis, dist_bnd - real(dp) :: f(2), J(2, 2), det, delta(2) - real(dp) :: res_norm, res_norm_try, alpha - real(dp) :: s_try, theta_try - integer :: iter, k - - associate(dummy => self) - end associate - - ierr = 0 - varphi = xcyl(2) - - call splint_vmec_data(0.0_dp, 0.0_dp, varphi, A_phi, A_theta, & - dA_phi_ds, dA_theta_ds, aiota, R_axis, Z_axis, & - alam, dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - - rs = xcyl(1) - R_axis - zs = xcyl(3) - Z_axis - theta = modulo(atan2(zs, rs), TWOPI) - - call splint_vmec_data(1.0_dp, theta, varphi, A_phi, A_theta, dA_phi_ds, & - dA_theta_ds, aiota, R_bnd, Z_bnd, alam, dR_ds, & - dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp, dl_ds, dl_dt, & - dl_dp) - - dist_axis = sqrt(rs**2 + zs**2) - dist_bnd = sqrt((R_bnd - R_axis)**2 + (Z_bnd - Z_axis)**2) - if (dist_bnd > 1.0e-12_dp) then - s = min(1.0_dp, max(0.0_dp, (dist_axis/dist_bnd)**2)) - else - s = 0.0_dp - end if - - do iter = 1, max_iter - call splint_vmec_data(s, theta, varphi, A_phi, A_theta, dA_phi_ds, & - dA_theta_ds, aiota, R, Z, alam, dR_ds, & - dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp, dl_ds, & - dl_dt, dl_dp) - - f(1) = R - xcyl(1) - f(2) = Z - xcyl(3) - res_norm = sqrt(f(1)**2 + f(2)**2) - if (res_norm < tol_res) exit - - J(1, 1) = dR_ds - J(1, 2) = dR_dt - J(2, 1) = dZ_ds - J(2, 2) = dZ_dt - det = J(1, 1)*J(2, 2) - J(1, 2)*J(2, 1) - if (abs(det) < 1.0e-14_dp) then - ierr = 2 - return - end if - - delta(1) = (-f(1)*J(2, 2) + f(2)*J(1, 2))/det - delta(2) = (-J(1, 1)*f(2) + J(2, 1)*f(1))/det - - if (maxval(abs(delta)) < tol_step) exit - - alpha = 1.0_dp - do k = 1, 10 - s_try = min(1.0_dp, max(0.0_dp, s + alpha*delta(1))) - theta_try = modulo(theta + alpha*delta(2), TWOPI) - call splint_vmec_data(s_try, theta_try, varphi, A_phi, A_theta, & - dA_phi_ds, dA_theta_ds, aiota, R, Z, alam, & - dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp, & - dl_ds, dl_dt, dl_dp) - res_norm_try = sqrt((R - xcyl(1))**2 + (Z - xcyl(3))**2) - if (res_norm_try < res_norm) then - s = s_try - theta = theta_try - exit - end if - alpha = 0.5_dp*alpha - end do - - if (k > 10) then - ierr = 1 - return - end if - end do - - if (iter > max_iter) then - ierr = 1 - return - end if - - u = [s, theta, varphi] - end subroutine vmec_from_cyl - - ! ======================================================================== - ! Geoflux coordinate system procedures - ! ======================================================================== - - subroutine make_geoflux_coordinate_system(cs) - class(coordinate_system_t), allocatable, intent(out) :: cs - allocate(geoflux_coordinate_system_t :: cs) - end subroutine make_geoflux_coordinate_system - - subroutine geoflux_evaluate_cyl(self, u, x) - class(geoflux_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - - associate(dummy => self) - end associate - - call geoflux_to_cyl(u, x) - end subroutine geoflux_evaluate_cyl - - subroutine geoflux_evaluate_cart(self, u, x) - class(geoflux_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: x(3) - - real(dp) :: xcyl(3) - - call self%evaluate_cyl(u, xcyl) - call cyl_to_cart(xcyl, x) - end subroutine geoflux_evaluate_cart - - subroutine geoflux_covariant_basis(self, u, e_cov) - class(geoflux_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3,3) - - real(dp) :: s, theta, phi - real(dp) :: R, Z - real(dp) :: jac_cyl(3,3) - real(dp) :: cos_phi, sin_phi - real(dp) :: xgeo(3), xcyl(3) - - associate(dummy => self) - end associate - - s = u(1) - theta = u(2) - phi = u(3) - - xgeo = [s, theta, phi] - call geoflux_to_cyl(xgeo, xcyl) - R = xcyl(1) - Z = xcyl(3) - - call assign_geoflux_to_cyl_jacobian(s, theta, phi, R, Z, jac_cyl) - - cos_phi = cos(phi) - sin_phi = sin(phi) - - e_cov(1, 1) = jac_cyl(1,1) * cos_phi - e_cov(2, 1) = jac_cyl(1,1) * sin_phi - e_cov(3, 1) = jac_cyl(3,1) - - e_cov(1, 2) = jac_cyl(1,2) * cos_phi - e_cov(2, 2) = jac_cyl(1,2) * sin_phi - e_cov(3, 2) = jac_cyl(3,2) - - e_cov(1, 3) = jac_cyl(1,3) * cos_phi - R * sin_phi - e_cov(2, 3) = jac_cyl(1,3) * sin_phi + R * cos_phi - e_cov(3, 3) = jac_cyl(3,3) - end subroutine geoflux_covariant_basis - - subroutine geoflux_metric_tensor(self, u, g, ginv, sqrtg) - class(geoflux_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg - - real(dp) :: e_cov(3,3) - real(dp) :: det - integer :: i, j - - call self%covariant_basis(u, e_cov) - - do i = 1, 3 - do j = 1, 3 - g(i, j) = dot_product(e_cov(:, i), e_cov(:, j)) - end do - end do - - det = g(1,1)*(g(2,2)*g(3,3) - g(2,3)*g(3,2)) & - - g(1,2)*(g(2,1)*g(3,3) - g(2,3)*g(3,1)) & - + g(1,3)*(g(2,1)*g(3,2) - g(2,2)*g(3,1)) - - sqrtg = sqrt(abs(det)) - - ginv(1,1) = (g(2,2)*g(3,3) - g(2,3)*g(3,2))/det - ginv(1,2) = (g(1,3)*g(3,2) - g(1,2)*g(3,3))/det - ginv(1,3) = (g(1,2)*g(2,3) - g(1,3)*g(2,2))/det - ginv(2,1) = (g(2,3)*g(3,1) - g(2,1)*g(3,3))/det - ginv(2,2) = (g(1,1)*g(3,3) - g(1,3)*g(3,1))/det - ginv(2,3) = (g(1,3)*g(2,1) - g(1,1)*g(2,3))/det - ginv(3,1) = (g(2,1)*g(3,2) - g(2,2)*g(3,1))/det - ginv(3,2) = (g(1,2)*g(3,1) - g(1,1)*g(3,2))/det - ginv(3,3) = (g(1,1)*g(2,2) - g(1,2)*g(2,1))/det - end subroutine geoflux_metric_tensor - - subroutine geoflux_from_cyl(self, xcyl, u, ierr) - class(geoflux_coordinate_system_t), intent(in) :: self - real(dp), intent(in) :: xcyl(3) - real(dp), intent(out) :: u(3) - integer, intent(out) :: ierr - - associate(dummy => self) - end associate - associate(dummy2 => xcyl) - end associate - u = 0.0_dp - ierr = 0 - - error stop "geoflux_from_cyl: not implemented" - end subroutine geoflux_from_cyl - - ! ======================================================================== - ! File detection procedure - ! ======================================================================== - - subroutine detect_refcoords_file_type(filename, file_type, ierr, message) - character(len=*), intent(in) :: filename - integer, intent(out) :: file_type - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: ncid - integer :: status - integer :: dim_rho, dim_theta, dim_zeta - integer :: var_x, var_y, var_z - integer :: var_rmnc - integer :: ierr_val - character(len=2048) :: message_val - - ierr = 0 - message = "" - file_type = refcoords_file_unknown - - status = nf90_open(trim(filename), NF90_NOWRITE, ncid) - if (status /= NF90_NOERR) then - ierr = 1 - message = "open failed: "//trim(nf90_strerror(status)) - return - end if - - if (nf90_inq_dimid(ncid, "rho", dim_rho) == NF90_NOERR .and. & - nf90_inq_dimid(ncid, "theta", dim_theta) == NF90_NOERR .and. & - nf90_inq_dimid(ncid, "zeta", dim_zeta) == NF90_NOERR .and. & - nf90_inq_varid(ncid, "x", var_x) == NF90_NOERR .and. & - nf90_inq_varid(ncid, "y", var_y) == NF90_NOERR .and. & - nf90_inq_varid(ncid, "z", var_z) == NF90_NOERR) then - status = nf90_close(ncid) - file_type = refcoords_file_chartmap - - call validate_chartmap_file(filename, ierr_val, message_val) - if (ierr_val /= 0) then - ierr = ierr_val - message = trim(message_val) - end if - return - end if - - if (nf90_inq_varid(ncid, "rmnc", var_rmnc) == NF90_NOERR) then - file_type = refcoords_file_vmec_wout - status = nf90_close(ncid) - return - end if - - status = nf90_close(ncid) - end subroutine detect_refcoords_file_type - - ! ======================================================================== - ! Chartmap validator procedures - ! ======================================================================== - - subroutine validate_chartmap_file(filename, ierr, message) - character(len=*), intent(in) :: filename - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer, parameter :: ok = 0 - integer, parameter :: err_open = 1 - integer, parameter :: err_missing_dim = 2 - integer, parameter :: err_invalid_dim = 3 - integer, parameter :: err_missing_var = 4 - integer, parameter :: err_bad_var = 5 - integer, parameter :: err_missing_units = 6 - integer, parameter :: err_bad_units = 7 - integer, parameter :: err_bad_grid = 8 - - integer :: status - integer :: ncid - integer :: dim_rho, dim_theta, dim_zeta - integer :: len_rho, len_theta, len_zeta - integer :: var_rho, var_theta, var_zeta, var_x, var_y, var_z - integer :: num_field_periods - real(dp), allocatable :: rho(:), theta(:), zeta(:) - real(dp) :: period - character(len=16) :: zeta_convention - - ierr = ok - message = "" - - status = nf90_open(trim(filename), NF90_NOWRITE, ncid) - if (status /= NF90_NOERR) then - ierr = err_open - message = "open failed: "//trim(nf90_strerror(status)) - return - end if - - do - call validator_require_dim(ncid, "rho", dim_rho, len_rho, ierr, message) - if (ierr /= ok) exit - call validator_require_dim(ncid, "theta", dim_theta, len_theta, ierr, & - message) - if (ierr /= ok) exit - call validator_require_dim(ncid, "zeta", dim_zeta, len_zeta, ierr, message) - if (ierr /= ok) exit - - if (len_rho < 2 .or. len_theta < 2 .or. len_zeta < 2) then - ierr = err_invalid_dim - message = "invalid dimension length(s)" - exit - end if - - call validator_require_1d_real64(ncid, "rho", dim_rho, var_rho, ierr, & - message) - if (ierr /= ok) exit - call validator_require_1d_real64(ncid, "theta", dim_theta, var_theta, & - ierr, message) - if (ierr /= ok) exit - call validator_require_1d_real64(ncid, "zeta", dim_zeta, var_zeta, ierr, & - message) - if (ierr /= ok) exit - - call validator_require_3d_real64(ncid, "x", dim_rho, dim_theta, dim_zeta, & - var_x, ierr, message) - if (ierr /= ok) exit - call validator_require_3d_real64(ncid, "y", dim_rho, dim_theta, dim_zeta, & - var_y, ierr, message) - if (ierr /= ok) exit - call validator_require_3d_real64(ncid, "z", dim_rho, dim_theta, dim_zeta, & - var_z, ierr, message) - if (ierr /= ok) exit - - call validator_require_units_cm(ncid, "x", var_x, ierr, message) - if (ierr /= ok) exit - call validator_require_units_cm(ncid, "y", var_y, ierr, message) - if (ierr /= ok) exit - call validator_require_units_cm(ncid, "z", var_z, ierr, message) - if (ierr /= ok) exit - - call validator_check_optional_zeta_convention(ncid, zeta_convention, & - ierr, message) - if (ierr /= ok) exit - - call validator_read_optional_num_field_periods(ncid, num_field_periods, & - ierr, message) - if (ierr /= ok) exit - - allocate (rho(len_rho), theta(len_theta), zeta(len_zeta)) - call validator_read_grid(ncid, var_rho, "rho", rho, ierr, message) - if (ierr /= ok) exit - call validator_read_grid(ncid, var_theta, "theta", theta, ierr, message) - if (ierr /= ok) exit - call validator_read_grid(ncid, var_zeta, "zeta", zeta, ierr, message) - if (ierr /= ok) exit - - call validator_check_rho_grid(rho, ierr, message) - if (ierr /= ok) exit - - call validator_check_periodic_grid("theta", theta, TWOPI, ierr, message) - if (ierr /= ok) exit - - period = TWOPI/real(num_field_periods, dp) - call validator_check_periodic_grid("zeta", zeta, period, ierr, message) - if (ierr /= ok) exit - - if (trim(zeta_convention) == "cyl" .or. trim(zeta_convention) == & - "vmec") then - call validator_check_cyl_phi_contract(ncid, var_x, var_y, len_rho, & - len_theta, zeta, & - period, ierr, message) - if (ierr /= ok) exit - end if - - exit - end do - - status = nf90_close(ncid) - end subroutine validate_chartmap_file - - subroutine validator_check_cyl_phi_contract(ncid, var_x, var_y, len_rho, & - len_theta, zeta, period, ierr, & - message) - integer, intent(in) :: ncid - integer, intent(in) :: var_x, var_y - integer, intent(in) :: len_rho, len_theta - real(dp), intent(in) :: zeta(:) - real(dp), intent(in) :: period - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer, parameter :: err_bad_grid = 8 - real(dp), parameter :: tol_phi = 1.0e-10_dp - integer :: ir, it, iz - integer :: i - integer, dimension(2) :: ir_list - integer, dimension(3) :: it_list - integer, dimension(4) :: iz_list - real(dp) :: xbuf(1, 1, 1), ybuf(1, 1, 1) - real(dp) :: phi, dphi, r_xy - integer :: status - - ierr = 0 - message = "" - - if (size(zeta) < 2) then - ierr = err_bad_grid - message = "zeta_convention=cyl requires at least 2 zeta points" - return - end if - - ir_list = [max(1, 1 + len_rho/2), len_rho] - it_list = [1, 1 + len_theta/3, 1 + (2*len_theta)/3] - iz_list = [1, 2, 1 + size(zeta)/2, size(zeta)] - - do ir = 1, size(ir_list) - do i = 1, size(iz_list) - iz = max(1, min(size(zeta), iz_list(i))) - do it = 1, size(it_list) - status = nf90_get_var(ncid, var_x, xbuf, & - start=[ir_list(ir), it_list(it), iz], & - count=[1, 1, 1]) - if (status /= NF90_NOERR) then - ierr = err_bad_grid - message = "could not read x for cyl contract check" - return - end if - status = nf90_get_var(ncid, var_y, ybuf, & - start=[ir_list(ir), it_list(it), iz], & - count=[1, 1, 1]) - if (status /= NF90_NOERR) then - ierr = err_bad_grid - message = "could not read y for cyl contract check" - return - end if - - r_xy = sqrt(xbuf(1, 1, 1)**2 + ybuf(1, 1, 1)**2) - if (r_xy < 1.0e-12_dp) cycle - - phi = atan2(ybuf(1, 1, 1), xbuf(1, 1, 1)) - dphi = abs(modulo(phi - zeta(iz) + 0.5_dp*period, period) - & - 0.5_dp*period) - - if (dphi > tol_phi) then - ierr = err_bad_grid - write (message, '(a,1x,es12.4)') & - "cyl requires atan2(y,x)=zeta; dphi=", dphi - return - end if - end do - end do - end do - end subroutine validator_check_cyl_phi_contract - - subroutine validator_check_optional_zeta_convention(ncid, value_out, ierr, & - message) - integer, intent(in) :: ncid - character(len=*), intent(out) :: value_out - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status - integer :: att_type - integer :: att_len - character(len=:), allocatable :: value - - ierr = 0 - message = "" - value_out = "" - - status = nf90_inquire_attribute(ncid, NF90_GLOBAL, "zeta_convention", & - xtype=att_type, len=att_len) - if (status /= NF90_NOERR) then - ierr = 5 - message = "missing zeta_convention global attribute" - return - end if - if (att_type /= NF90_CHAR .or. att_len < 1) then - ierr = 5 - message = "zeta_convention must be a global string attribute" - return - end if - - allocate (character(len=att_len) :: value) - status = nf90_get_att(ncid, NF90_GLOBAL, "zeta_convention", value) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not read zeta_convention" - return - end if - - select case (trim(value)) - case ("cyl", "vmec") - value_out = trim(value) - case default - ierr = 5 - message = "unsupported zeta_convention (must be cyl or vmec)" - end select - end subroutine validator_check_optional_zeta_convention - - subroutine validator_require_dim(ncid, name, dimid, dimlen, ierr, message) - integer, intent(in) :: ncid - character(len=*), intent(in) :: name - integer, intent(out) :: dimid - integer, intent(out) :: dimlen - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status - - ierr = 0 - message = "" - - status = nf90_inq_dimid(ncid, trim(name), dimid) - if (status /= NF90_NOERR) then - ierr = 2 - message = "missing dimension "//trim(name) - return - end if - - status = nf90_inquire_dimension(ncid, dimid, len=dimlen) - if (status /= NF90_NOERR) then - ierr = 3 - message = "invalid dimension "//trim(name) - return - end if - end subroutine validator_require_dim - - subroutine validator_require_1d_real64(ncid, name, dimid, varid, ierr, message) - integer, intent(in) :: ncid - character(len=*), intent(in) :: name - integer, intent(in) :: dimid - integer, intent(out) :: varid - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status, xtype, ndims - integer :: dimids(NF90_MAX_VAR_DIMS) - - ierr = 0 - message = "" - - status = nf90_inq_varid(ncid, trim(name), varid) - if (status /= NF90_NOERR) then - ierr = 4 - message = "missing variable "//trim(name) - return - end if - - status = nf90_inquire_variable(ncid, varid, xtype=xtype, ndims=ndims, & - dimids=dimids) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not inquire variable "//trim(name) - return - end if - - if (xtype /= NF90_DOUBLE .or. ndims /= 1 .or. dimids(1) /= dimid) then - ierr = 5 - message = "variable "//trim(name)//" must be float64 with correct dim" - return - end if - end subroutine validator_require_1d_real64 - - subroutine validator_require_3d_real64(ncid, name, dim_rho, dim_theta, dim_zeta, & - varid, ierr, message) - integer, intent(in) :: ncid - character(len=*), intent(in) :: name - integer, intent(in) :: dim_rho, dim_theta, dim_zeta - integer, intent(out) :: varid - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status, xtype, ndims - integer :: dimids(NF90_MAX_VAR_DIMS) - - ierr = 0 - message = "" - - status = nf90_inq_varid(ncid, trim(name), varid) - if (status /= NF90_NOERR) then - ierr = 4 - message = "missing variable "//trim(name) - return - end if - - status = nf90_inquire_variable(ncid, varid, xtype=xtype, ndims=ndims, & - dimids=dimids) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not inquire variable "//trim(name) - return - end if - - if (xtype /= NF90_DOUBLE .or. ndims /= 3) then - ierr = 5 - message = "variable "//trim(name)//" must be float64 with 3 dimensions" - return - end if - - if (dimids(1) /= dim_rho .or. dimids(2) /= dim_theta .or. dimids(3) /= & - dim_zeta) then - ierr = 5 - message = "variable "//trim(name)// & - " has wrong dimension order; expected (rho,theta,zeta)" - return - end if - end subroutine validator_require_3d_real64 - - subroutine validator_require_units_cm(ncid, name, varid, ierr, message) - integer, intent(in) :: ncid - character(len=*), intent(in) :: name - integer, intent(in) :: varid - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status, att_len, att_type - character(len=:), allocatable :: units - - ierr = 0 - message = "" - - status = nf90_inquire_attribute(ncid, varid, "units", xtype=att_type, & - len=att_len) - if (status /= NF90_NOERR) then - ierr = 6 - message = "missing units attribute for "//trim(name) - return - end if - - if (att_type /= NF90_CHAR .or. att_len < 1) then - ierr = 7 - message = "invalid units attribute for "//trim(name) - return - end if - - allocate (character(len=att_len) :: units) - status = nf90_get_att(ncid, varid, "units", units) - if (status /= NF90_NOERR) then - ierr = 7 - message = "could not read units attribute for "//trim(name) - return - end if - - if (trim(units) /= "cm") then - ierr = 7 - message = "units for "//trim(name)//" must be cm" - return - end if - end subroutine validator_require_units_cm - - subroutine validator_read_optional_num_field_periods(ncid, num_field_periods, & - ierr, message) - integer, intent(in) :: ncid - integer, intent(out) :: num_field_periods - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status, varid, xtype, ndims - character(len=*), parameter :: var_primary = "num_field_periods" - character(len=*), parameter :: var_legacy = "nfp" - - ierr = 0 - message = "" - num_field_periods = 1 - - status = nf90_inq_varid(ncid, var_primary, varid) - if (status /= NF90_NOERR) then - status = nf90_inq_varid(ncid, var_legacy, varid) - if (status /= NF90_NOERR) return - end if - - status = nf90_inquire_variable(ncid, varid, xtype=xtype, ndims=ndims) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not inquire variable num_field_periods" - return - end if - - if (ndims /= 0 .or. .not. validator_is_integer_xtype(xtype)) then - ierr = 5 - message = "num_field_periods must be a scalar integer variable" - return - end if - - status = nf90_get_var(ncid, varid, num_field_periods) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not read num_field_periods" - return - end if - - if (num_field_periods < 1) then - ierr = 5 - message = "num_field_periods must be >= 1" - return - end if - end subroutine validator_read_optional_num_field_periods - - subroutine validator_read_grid(ncid, varid, name, x, ierr, message) - integer, intent(in) :: ncid - integer, intent(in) :: varid - character(len=*), intent(in) :: name - real(dp), intent(out) :: x(:) - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - integer :: status - - ierr = 0 - message = "" - - status = nf90_get_var(ncid, varid, x) - if (status /= NF90_NOERR) then - ierr = 5 - message = "could not read "//trim(name) - end if - end subroutine validator_read_grid - - subroutine validator_check_rho_grid(rho, ierr, message) - real(dp), intent(in) :: rho(:) - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - real(dp), parameter :: tol = 1.0e-12_dp - real(dp), parameter :: rho_min_max = 0.01_dp - - ierr = 0 - message = "" - - if (.not. validator_is_strictly_increasing(rho)) then - ierr = 8 - message = "rho grid must be strictly increasing" - return - end if - - if (rho(1) < 0.0_dp .or. rho(1) > rho_min_max) then - ierr = 8 - message = "rho(1) must be in [0, 0.01] (axis or small positive for singularity avoidance)" - return - end if - - if (abs(rho(size(rho)) - 1.0_dp) > tol) then - ierr = 8 - message = "rho must end at 1.0" - return - end if - end subroutine validator_check_rho_grid - - subroutine validator_check_periodic_grid(name, x, period, ierr, message) - character(len=*), intent(in) :: name - real(dp), intent(in) :: x(:) - real(dp), intent(in) :: period - integer, intent(out) :: ierr - character(len=*), intent(out) :: message - - real(dp), parameter :: tol = 1.0e-12_dp - real(dp) :: step_expected - - ierr = 0 - message = "" - - if (size(x) == 1) then - if (abs(x(1) - 0.0_dp) > tol) then - ierr = 8 - message = trim(name)//" must start at 0" - end if - return - end if - - if (.not. validator_is_strictly_increasing(x)) then - ierr = 8 - message = trim(name)//" grid must be strictly increasing" - return - end if - - step_expected = period/real(size(x), dp) - if (abs(x(1) - 0.0_dp) > tol) then - ierr = 8 - message = trim(name)//" must start at 0" - return - end if - - if (.not. validator_is_uniform_step(x, step_expected, tol)) then - ierr = 8 - message = trim(name)//" must be uniform with endpoint excluded" - return - end if - - if (x(size(x)) > (period - step_expected + tol)) then - ierr = 8 - message = trim(name)//" must exclude the period endpoint" - return - end if - end subroutine validator_check_periodic_grid - - logical function validator_is_integer_xtype(xtype) - integer, intent(in) :: xtype - validator_is_integer_xtype = (xtype == NF90_BYTE) .or. & - (xtype == NF90_SHORT) .or. & - (xtype == NF90_INT) .or. (xtype == NF90_INT64) - end function validator_is_integer_xtype - - logical function validator_is_strictly_increasing(x) - real(dp), intent(in) :: x(:) - integer :: i - - validator_is_strictly_increasing = .true. - do i = 2, size(x) - if (x(i) <= x(i - 1)) then - validator_is_strictly_increasing = .false. - return - end if - end do - end function validator_is_strictly_increasing - - logical function validator_is_uniform_step(x, step, tol) - real(dp), intent(in) :: x(:) - real(dp), intent(in) :: step - real(dp), intent(in) :: tol - - real(dp) :: dx_max - - if (size(x) < 2) then - validator_is_uniform_step = .true. - return - end if - - dx_max = maxval(abs((x(2:) - x(:size(x) - 1)) - step)) - validator_is_uniform_step = (dx_max <= 10.0_dp*tol) - end function validator_is_uniform_step - - ! ======================================================================== - ! Chartmap coordinate system procedures - include from file - ! ======================================================================== - -#include "libneo_coordinates_chartmap_impl.inc" - -end module libneo_coordinates diff --git a/src/coordinates/libneo_coordinates_chartmap_validator.f90 b/src/coordinates/libneo_coordinates_validator.f90 similarity index 98% rename from src/coordinates/libneo_coordinates_chartmap_validator.f90 rename to src/coordinates/libneo_coordinates_validator.f90 index ab3685d8..980cfe3a 100644 --- a/src/coordinates/libneo_coordinates_chartmap_validator.f90 +++ b/src/coordinates/libneo_coordinates_validator.f90 @@ -1,4 +1,5 @@ -submodule(libneo_coordinates) libneo_coordinates_chartmap_validator +module libneo_coordinates_validator + use, intrinsic :: iso_fortran_env, only: dp => real64 use math_constants, only: TWOPI use netcdf, only: NF90_BYTE, NF90_CHAR, NF90_DOUBLE, NF90_GLOBAL, NF90_INT, & NF90_INT64, NF90_MAX_VAR_DIMS, NF90_NOERR, NF90_NOWRITE, & @@ -9,10 +10,13 @@ nf90_inquire_variable, & nf90_open, nf90_strerror implicit none + private + + public :: validate_chartmap_file contains - module subroutine validate_chartmap_file(filename, ierr, message) + subroutine validate_chartmap_file(filename, ierr, message) character(len=*), intent(in) :: filename integer, intent(out) :: ierr character(len=*), intent(out) :: message @@ -571,4 +575,4 @@ logical function is_uniform_step(x, step, tol) is_uniform_step = (dx_max <= 10.0_dp*tol) end function is_uniform_step -end submodule libneo_coordinates_chartmap_validator +end module libneo_coordinates_validator diff --git a/src/coordinates/libneo_coordinates_vmec.f90 b/src/coordinates/libneo_coordinates_vmec.f90 index d22772f9..65dd2125 100644 --- a/src/coordinates/libneo_coordinates_vmec.f90 +++ b/src/coordinates/libneo_coordinates_vmec.f90 @@ -1,12 +1,26 @@ -submodule (libneo_coordinates) libneo_coordinates_vmec +module libneo_coordinates_vmec + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates_base, only: coordinate_system_t use spline_vmec_sub, only: splint_vmec_data use cylindrical_cartesian, only: cyl_to_cart use math_constants, only: TWOPI implicit none + private + + public :: vmec_coordinate_system_t, make_vmec_coordinate_system + + type, extends(coordinate_system_t) :: vmec_coordinate_system_t + contains + procedure :: evaluate_cart => vmec_evaluate_cart + procedure :: evaluate_cyl => vmec_evaluate_cyl + procedure :: covariant_basis => vmec_covariant_basis + procedure :: metric_tensor => vmec_metric_tensor + procedure :: from_cyl => vmec_from_cyl + end type vmec_coordinate_system_t contains - module subroutine make_vmec_coordinate_system(cs) + subroutine make_vmec_coordinate_system(cs) class(coordinate_system_t), allocatable, intent(out) :: cs allocate(vmec_coordinate_system_t :: cs) end subroutine make_vmec_coordinate_system @@ -21,6 +35,9 @@ subroutine vmec_evaluate_cyl(self, u, x) real(dp) :: R, Z, alam, dR_ds, dR_dt, dR_dp, dZ_ds, dZ_dt, dZ_dp real(dp) :: dl_ds, dl_dt, dl_dp + associate(dummy => self) + end associate + s = u(1) theta = u(2) varphi = u(3) @@ -48,7 +65,7 @@ end subroutine vmec_evaluate_cart subroutine vmec_covariant_basis(self, u, e_cov) class(vmec_coordinate_system_t), intent(in) :: self real(dp), intent(in) :: u(3) - real(dp), intent(out) :: e_cov(3,3) + real(dp), intent(out) :: e_cov(3, 3) real(dp) :: s, theta, varphi real(dp) :: A_phi, A_theta, dA_phi_ds, dA_theta_ds, aiota @@ -56,6 +73,9 @@ subroutine vmec_covariant_basis(self, u, e_cov) real(dp) :: dl_ds, dl_dt, dl_dp real(dp) :: cos_phi, sin_phi + associate(dummy => self) + end associate + s = u(1) theta = u(2) varphi = u(3) @@ -83,9 +103,9 @@ end subroutine vmec_covariant_basis subroutine vmec_metric_tensor(self, u, g, ginv, sqrtg) class(vmec_coordinate_system_t), intent(in) :: self real(dp), intent(in) :: u(3) - real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg + real(dp), intent(out) :: g(3, 3), ginv(3, 3), sqrtg - real(dp) :: e_cov(3,3) + real(dp) :: e_cov(3, 3) real(dp) :: det integer :: i, j @@ -221,4 +241,4 @@ subroutine vmec_from_cyl(self, xcyl, u, ierr) u = [s, theta, varphi] end subroutine vmec_from_cyl -end submodule libneo_coordinates_vmec +end module libneo_coordinates_vmec diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 32d82149..f3dcf7cf 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -264,6 +264,16 @@ if(LIBNEO_TESTING_ENABLED) --basename geoflux ) set_tests_properties(test_geoflux_plot PROPERTIES DEPENDS test_geoflux) + + add_test( + NAME test_geoflux_coordinate_system + COMMAND ${Python_EXECUTABLE} ${CMAKE_SOURCE_DIR}/test/scripts/run_geoflux_test.py + --exe $ + --data-dir ${TEST_OUTPUT_DIR}/scripts + ) + set_tests_properties(test_geoflux_coordinate_system PROPERTIES + FAIL_REGULAR_EXPRESSION "STOP" + DEPENDS test_geoflux) endif() ## Custom tests @@ -444,6 +454,18 @@ set_tests_properties(test_vmec_from_cyl PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} DEPENDS setup_vmec_wout) +add_executable(test_vector_conversion.x source/test_vector_conversion.f90) +target_link_libraries(test_vector_conversion.x neo ${MAGFIE_LIB}) +add_test(NAME test_vector_conversion + COMMAND test_vector_conversion.x) +set_tests_properties(test_vector_conversion PROPERTIES + FAIL_REGULAR_EXPRESSION "STOP" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + DEPENDS setup_chartmap_volume) + +add_executable(test_geoflux_coordinate_system.x source/test_geoflux_coordinate_system.f90) +target_link_libraries(test_geoflux_coordinate_system.x neo ${MAGFIE_LIB}) + add_executable(test_chartmap_matches_vmec.x source/test_chartmap_matches_vmec.f90) target_link_libraries(test_chartmap_matches_vmec.x neo) add_test(NAME test_chartmap_matches_vmec diff --git a/test/source/test_geoflux_coordinate_system.f90 b/test/source/test_geoflux_coordinate_system.f90 new file mode 100644 index 00000000..b8499d71 --- /dev/null +++ b/test/source/test_geoflux_coordinate_system.f90 @@ -0,0 +1,310 @@ +program test_geoflux_coordinate_system + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates, only: coordinate_system_t, & + make_geoflux_coordinate_system, & + geoflux_coordinate_system_t + use geoflux_coordinates, only: geoflux_to_cyl, cyl_to_geoflux + use geoflux_field, only: spline_geoflux_data + use cylindrical_cartesian, only: cyl_to_cart + implicit none + + character(len=*), parameter :: fallback_geqdsk = '../../python/tests/test.geqdsk' + character(len=*), parameter :: env_geqdsk = 'LIBNEO_TEST_GEQDSK' + + character(len=512) :: geqdsk_file + character(len=512) :: arg_buffer + integer :: arg_status, arg_len + integer :: nerrors + class(coordinate_system_t), allocatable :: cs + + nerrors = 0 + + geqdsk_file = fallback_geqdsk + call get_command_argument(1, value=arg_buffer, length=arg_len, status=arg_status) + if (arg_status == 0 .and. arg_len > 0 .and. len_trim(arg_buffer) > 0) then + geqdsk_file = trim(arg_buffer) + else + call get_environment_variable(env_geqdsk, value=arg_buffer, & + length=arg_len, status=arg_status) + if (arg_status == 0 .and. arg_len > 0 .and. len_trim(arg_buffer) > 0) then + geqdsk_file = trim(arg_buffer) + end if + end if + + call spline_geoflux_data(trim(geqdsk_file), 64, 128) + + call make_geoflux_coordinate_system(cs) + if (.not. allocated(cs)) error stop "Failed to create geoflux coordinate system" + + select type (gcs => cs) + type is (geoflux_coordinate_system_t) + call test_evaluate_cyl_matches_module(gcs, nerrors) + call test_evaluate_cart_consistency(gcs, nerrors) + call test_coordinate_roundtrip(gcs, nerrors) + call test_covariant_basis_finite_difference(gcs, nerrors) + call test_metric_tensor_properties(gcs, nerrors) + class default + error stop "Unexpected coordinate system type" + end select + + if (nerrors > 0) then + print *, "FAILED:", nerrors, "error(s) in geoflux coordinate system tests" + error stop 1 + end if + print *, "All geoflux coordinate system tests passed!" + +contains + + subroutine test_evaluate_cyl_matches_module(cs, nerrors) + class(geoflux_coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), xcyl_oop(3), xcyl_mod(3) + real(dp) :: tol, err + integer :: i, n_passed + real(dp) :: test_s(3), test_theta(3), test_phi(2) + integer :: is, it, ip + + tol = 1.0e-14_dp + n_passed = 0 + + test_s = [0.1_dp, 0.5_dp, 0.9_dp] + test_theta = [0.3_dp, 1.5_dp, -0.8_dp] + test_phi = [0.0_dp, 1.2_dp] + + do ip = 1, 2 + do it = 1, 3 + do is = 1, 3 + u = [test_s(is), test_theta(it), test_phi(ip)] + + call cs%evaluate_cyl(u, xcyl_oop) + call geoflux_to_cyl(u, xcyl_mod) + + err = maxval(abs(xcyl_oop - xcyl_mod)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: evaluate_cyl mismatch at u=", u + print *, " OOP:", xcyl_oop + print *, " mod:", xcyl_mod + end if + end do + end do + end do + + if (n_passed == 18) then + print *, " PASS: evaluate_cyl matches geoflux_to_cyl (18 points)" + else + print *, " FAIL: evaluate_cyl mismatch", 18 - n_passed, "/18" + nerrors = nerrors + 1 + end if + end subroutine test_evaluate_cyl_matches_module + + subroutine test_evaluate_cart_consistency(cs, nerrors) + class(geoflux_coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), xcyl(3), xcart_oop(3), xcart_manual(3) + real(dp) :: tol, err + integer :: n_passed + real(dp) :: test_points(3, 4) + integer :: i + + tol = 1.0e-14_dp + n_passed = 0 + + test_points(:, 1) = [0.2_dp, 0.5_dp, 0.0_dp] + test_points(:, 2) = [0.6_dp, -1.0_dp, 0.8_dp] + test_points(:, 3) = [0.4_dp, 2.0_dp, 1.5_dp] + test_points(:, 4) = [0.8_dp, 0.1_dp, 2.3_dp] + + do i = 1, 4 + u = test_points(:, i) + + call cs%evaluate_cart(u, xcart_oop) + call cs%evaluate_cyl(u, xcyl) + call cyl_to_cart(xcyl, xcart_manual) + + err = maxval(abs(xcart_oop - xcart_manual)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: evaluate_cart inconsistent at u=", u + print *, " OOP:", xcart_oop + print *, " manual:", xcart_manual + end if + end do + + if (n_passed == 4) then + print *, " PASS: evaluate_cart = cyl_to_cart(evaluate_cyl)" + else + print *, " FAIL: evaluate_cart consistency", 4 - n_passed, "/4" + nerrors = nerrors + 1 + end if + end subroutine test_evaluate_cart_consistency + + subroutine test_coordinate_roundtrip(cs, nerrors) + class(geoflux_coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), xcyl(3), u_back(3), xcyl_back(3) + real(dp), parameter :: tol = 1.0e-8_dp + real(dp) :: err + integer :: n_passed + real(dp) :: test_points(3, 5) + integer :: i + + n_passed = 0 + + test_points(:, 1) = [0.2_dp, 0.3_dp, 0.0_dp] + test_points(:, 2) = [0.5_dp, 1.0_dp, 0.5_dp] + test_points(:, 3) = [0.7_dp, -0.5_dp, 1.2_dp] + test_points(:, 4) = [0.3_dp, 2.0_dp, 0.8_dp] + test_points(:, 5) = [0.9_dp, 0.1_dp, 0.0_dp] + + do i = 1, 5 + u = test_points(:, i) + + call cs%evaluate_cyl(u, xcyl) + call cyl_to_geoflux(xcyl, u_back) + call cs%evaluate_cyl(u_back, xcyl_back) + + err = maxval(abs(xcyl_back - xcyl)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: coordinate roundtrip at u=", u + print *, " u_back=", u_back, " err=", err + end if + end do + + if (n_passed == 5) then + print *, " PASS: u->cyl->u->cyl roundtrip (5 points)" + else + print *, " FAIL: coordinate roundtrip", 5 - n_passed, "/5" + nerrors = nerrors + 1 + end if + end subroutine test_coordinate_roundtrip + + subroutine test_covariant_basis_finite_difference(cs, nerrors) + class(geoflux_coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), e_cov(3, 3), e_cov_fd(3, 3) + real(dp) :: u_plus(3), u_minus(3), x_plus(3), x_minus(3) + real(dp), parameter :: h = 1.0e-6_dp + real(dp), parameter :: tol = 1.0e-4_dp + real(dp) :: rel_err, max_err + integer :: k, n_tests, n_passed + real(dp) :: test_points(3, 4) + integer :: i + + n_tests = 0 + n_passed = 0 + max_err = 0.0_dp + + test_points(:, 1) = [0.3_dp, 0.5_dp, 0.2_dp] + test_points(:, 2) = [0.5_dp, 1.2_dp, 0.0_dp] + test_points(:, 3) = [0.7_dp, -0.3_dp, 1.0_dp] + test_points(:, 4) = [0.4_dp, 2.5_dp, 0.7_dp] + + do i = 1, 4 + u = test_points(:, i) + n_tests = n_tests + 1 + + call cs%covariant_basis(u, e_cov) + + do k = 1, 3 + u_plus = u + u_minus = u + u_plus(k) = u(k) + h + u_minus(k) = u(k) - h + call cs%evaluate_cart(u_plus, x_plus) + call cs%evaluate_cart(u_minus, x_minus) + e_cov_fd(:, k) = (x_plus - x_minus) / (2.0_dp * h) + end do + + rel_err = maxval(abs(e_cov - e_cov_fd)) / max(maxval(abs(e_cov_fd)), 1.0e-10_dp) + if (rel_err > max_err) max_err = rel_err + + if (rel_err < tol) then + n_passed = n_passed + 1 + else + print *, " covariant_basis FD mismatch at u=", u, " rel_err=", rel_err + end if + end do + + if (n_passed == n_tests) then + print *, " PASS: covariant_basis matches FD (", n_tests, " points, max_err=", max_err, ")" + else + print *, " FAIL: covariant_basis FD check", n_tests - n_passed, "/", n_tests + nerrors = nerrors + 1 + end if + end subroutine test_covariant_basis_finite_difference + + subroutine test_metric_tensor_properties(cs, nerrors) + class(geoflux_coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), g(3, 3), ginv(3, 3), sqrtg + real(dp) :: identity(3, 3), prod(3, 3) + real(dp), parameter :: tol = 1.0e-10_dp + integer :: i, j, k + logical :: sym_ok, inv_ok, pos_ok + real(dp) :: test_points(3, 3) + integer :: ip, n_passed + + n_passed = 0 + + test_points(:, 1) = [0.3_dp, 0.8_dp, 0.1_dp] + test_points(:, 2) = [0.6_dp, 1.5_dp, 0.5_dp] + test_points(:, 3) = [0.5_dp, -0.5_dp, 1.2_dp] + + do ip = 1, 3 + u = test_points(:, ip) + call cs%metric_tensor(u, g, ginv, sqrtg) + + sym_ok = .true. + do i = 1, 3 + do j = i + 1, 3 + if (abs(g(i, j) - g(j, i)) > tol) sym_ok = .false. + if (abs(ginv(i, j) - ginv(j, i)) > tol) sym_ok = .false. + end do + end do + + identity = 0.0_dp + do i = 1, 3 + identity(i, i) = 1.0_dp + end do + + prod = 0.0_dp + do i = 1, 3 + do j = 1, 3 + do k = 1, 3 + prod(i, j) = prod(i, j) + g(i, k) * ginv(k, j) + end do + end do + end do + + inv_ok = maxval(abs(prod - identity)) < tol + pos_ok = sqrtg > 0.0_dp + + if (sym_ok .and. inv_ok .and. pos_ok) then + n_passed = n_passed + 1 + else + print *, " metric_tensor fail at u=", u + if (.not. sym_ok) print *, " not symmetric" + if (.not. inv_ok) print *, " g*ginv /= I" + if (.not. pos_ok) print *, " sqrtg <= 0" + end if + end do + + if (n_passed == 3) then + print *, " PASS: metric_tensor properties verified (3 points)" + else + print *, " FAIL: metric_tensor properties", 3 - n_passed, "/3" + nerrors = nerrors + 1 + end if + end subroutine test_metric_tensor_properties + +end program test_geoflux_coordinate_system diff --git a/test/source/test_vector_conversion.f90 b/test/source/test_vector_conversion.f90 new file mode 100644 index 00000000..9aa9fdc7 --- /dev/null +++ b/test/source/test_vector_conversion.f90 @@ -0,0 +1,332 @@ +program test_vector_conversion + use, intrinsic :: iso_fortran_env, only: dp => real64 + use libneo_coordinates, only: coordinate_system_t, & + make_chartmap_coordinate_system, & + chartmap_coordinate_system_t + use math_constants, only: TWOPI + implicit none + + integer :: nerrors + class(coordinate_system_t), allocatable :: cs + character(len=*), parameter :: chartmap_file = "chartmap.nc" + + nerrors = 0 + + call make_chartmap_coordinate_system(cs, chartmap_file) + if (.not. allocated(cs)) error stop "Failed to load chartmap" + + select type (ccs => cs) + type is (chartmap_coordinate_system_t) + call test_ctr_to_cart_basis_vectors(ccs, nerrors) + call test_ctr_to_cart_linearity(ccs, nerrors) + call test_cov_to_cart_norm_preservation(ccs, nerrors) + call test_cov_ctr_cart_consistency(ccs, nerrors) + call test_cart_ctr_cart_roundtrip(ccs, nerrors) + call test_cart_cov_cart_roundtrip(ccs, nerrors) + class default + error stop "Unexpected coordinate system type" + end select + + if (nerrors > 0) then + print *, "FAILED:", nerrors, "error(s) in vector conversion tests" + error stop 1 + end if + print *, "All vector conversion tests passed!" + +contains + + subroutine test_ctr_to_cart_basis_vectors(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), e_cov(3, 3) + real(dp) :: v_ctr(3), v_cart(3) + real(dp) :: tol, err + integer :: k + + tol = 1.0e-12_dp + u = [0.4_dp, 1.2_dp, 0.8_dp] + + call cs%covariant_basis(u, e_cov) + + do k = 1, 3 + v_ctr = 0.0_dp + v_ctr(k) = 1.0_dp + call cs%ctr_to_cart(u, v_ctr, v_cart) + err = maxval(abs(v_cart - e_cov(:, k))) + if (err > tol) then + print *, " FAIL: ctr_to_cart unit vector", k, "err=", err + nerrors = nerrors + 1 + end if + end do + + if (nerrors == 0) then + print *, " PASS: ctr_to_cart maps unit vectors to basis vectors" + end if + end subroutine test_ctr_to_cart_basis_vectors + + subroutine test_ctr_to_cart_linearity(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3) + real(dp) :: v1_ctr(3), v2_ctr(3), v_combined_ctr(3) + real(dp) :: v1_cart(3), v2_cart(3), v_combined_cart(3) + real(dp) :: v_expected(3) + real(dp) :: a, b, tol, err + + tol = 1.0e-12_dp + u = [0.5_dp, 2.5_dp, 1.0_dp] + v1_ctr = [1.3_dp, -0.7_dp, 2.1_dp] + v2_ctr = [-0.5_dp, 1.8_dp, 0.3_dp] + a = 2.5_dp + b = -1.3_dp + + call cs%ctr_to_cart(u, v1_ctr, v1_cart) + call cs%ctr_to_cart(u, v2_ctr, v2_cart) + + v_combined_ctr = a * v1_ctr + b * v2_ctr + call cs%ctr_to_cart(u, v_combined_ctr, v_combined_cart) + + v_expected = a * v1_cart + b * v2_cart + err = maxval(abs(v_combined_cart - v_expected)) + + if (err > tol) then + print *, " FAIL: ctr_to_cart linearity err=", err + nerrors = nerrors + 1 + else + print *, " PASS: ctr_to_cart is linear" + end if + end subroutine test_ctr_to_cart_linearity + + subroutine test_cov_to_cart_norm_preservation(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3) + real(dp) :: g(3, 3), ginv(3, 3), sqrtg + real(dp) :: v_cov(3), v_cart(3) + real(dp) :: norm_sq_curv, norm_sq_cart + real(dp) :: tol, rel_err + integer :: i, j, n_tests, n_passed + real(dp) :: test_points(3, 4) + + tol = 1.0e-10_dp + n_tests = 0 + n_passed = 0 + + test_points(:, 1) = [0.2_dp, 0.5_dp, 0.3_dp] + test_points(:, 2) = [0.7_dp, 3.0_dp, 1.5_dp] + test_points(:, 3) = [0.5_dp, 5.0_dp, 0.1_dp] + test_points(:, 4) = [0.9_dp, 1.0_dp, 2.0_dp] + + do i = 1, 4 + u = test_points(:, i) + call cs%metric_tensor(u, g, ginv, sqrtg) + + v_cov = [1.5_dp, -0.8_dp, 2.3_dp] + + norm_sq_curv = 0.0_dp + do j = 1, 3 + norm_sq_curv = norm_sq_curv + v_cov(j) * sum(ginv(j, :) * v_cov) + end do + + call cs%cov_to_cart(u, v_cov, v_cart) + norm_sq_cart = sum(v_cart**2) + + n_tests = n_tests + 1 + rel_err = abs(norm_sq_cart - norm_sq_curv) / max(norm_sq_curv, 1.0e-15_dp) + + if (rel_err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: norm mismatch at u=", u + print *, " |v|²_curv=", norm_sq_curv, " |v|²_cart=", norm_sq_cart + end if + end do + + if (n_passed == n_tests) then + print *, " PASS: cov_to_cart preserves vector norm (", n_tests, " points)" + else + print *, " FAIL: norm preservation failed", n_tests - n_passed, "/", n_tests + nerrors = nerrors + 1 + end if + end subroutine test_cov_to_cart_norm_preservation + + subroutine test_cov_ctr_cart_consistency(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3) + real(dp) :: g(3, 3), ginv(3, 3), sqrtg + real(dp) :: v_ctr(3), v_cov(3) + real(dp) :: v_cart_from_ctr(3), v_cart_from_cov(3) + real(dp) :: tol, err + integer :: i, n_tests, n_passed + real(dp) :: test_points(3, 3) + + tol = 1.0e-11_dp + n_tests = 0 + n_passed = 0 + + test_points(:, 1) = [0.3_dp, 1.0_dp, 0.5_dp] + test_points(:, 2) = [0.6_dp, 4.0_dp, 1.2_dp] + test_points(:, 3) = [0.8_dp, 2.5_dp, 0.8_dp] + + do i = 1, 3 + u = test_points(:, i) + + v_ctr = [0.7_dp, -1.2_dp, 0.9_dp] + + call cs%ctr_to_cart(u, v_ctr, v_cart_from_ctr) + + call cs%metric_tensor(u, g, ginv, sqrtg) + v_cov = matmul(g, v_ctr) + + call cs%cov_to_cart(u, v_cov, v_cart_from_cov) + + n_tests = n_tests + 1 + err = maxval(abs(v_cart_from_ctr - v_cart_from_cov)) + + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: cov/ctr inconsistency at u=", u, " err=", err + end if + end do + + if (n_passed == n_tests) then + print *, " PASS: cov_to_cart and ctr_to_cart are consistent (", n_tests, " points)" + else + print *, " FAIL: cov/ctr consistency failed", n_tests - n_passed, "/", n_tests + nerrors = nerrors + 1 + end if + end subroutine test_cov_ctr_cart_consistency + + subroutine test_cart_ctr_cart_roundtrip(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), e_cov(3, 3) + real(dp) :: v_cart(3), v_ctr(3), v_cart_back(3) + real(dp) :: tol, err + integer :: i, n_passed + real(dp) :: test_points(3, 4) + real(dp) :: test_vectors(3, 3) + integer :: iv + + tol = 1.0e-10_dp + n_passed = 0 + + test_points(:, 1) = [0.3_dp, 0.8_dp, 0.4_dp] + test_points(:, 2) = [0.6_dp, 2.5_dp, 1.0_dp] + test_points(:, 3) = [0.5_dp, 4.0_dp, 0.2_dp] + test_points(:, 4) = [0.8_dp, 1.5_dp, 1.8_dp] + + test_vectors(:, 1) = [1.0_dp, 0.0_dp, 0.0_dp] + test_vectors(:, 2) = [0.5_dp, -0.3_dp, 0.8_dp] + test_vectors(:, 3) = [-0.2_dp, 1.5_dp, 0.7_dp] + + do i = 1, 4 + u = test_points(:, i) + call cs%covariant_basis(u, e_cov) + + do iv = 1, 3 + v_cart = test_vectors(:, iv) + call solve_3x3(e_cov, v_cart, v_ctr) + call cs%ctr_to_cart(u, v_ctr, v_cart_back) + + err = maxval(abs(v_cart_back - v_cart)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: cart->ctr->cart roundtrip at u=", u, " err=", err + end if + end do + end do + + if (n_passed == 12) then + print *, " PASS: cart->ctr->cart roundtrip (12 tests)" + else + print *, " FAIL: cart->ctr->cart roundtrip", 12 - n_passed, "/12" + nerrors = nerrors + 1 + end if + end subroutine test_cart_ctr_cart_roundtrip + + subroutine test_cart_cov_cart_roundtrip(cs, nerrors) + class(coordinate_system_t), intent(in) :: cs + integer, intent(inout) :: nerrors + + real(dp) :: u(3), e_cov(3, 3) + real(dp) :: g(3, 3), ginv(3, 3), sqrtg + real(dp) :: v_cart(3), v_ctr(3), v_cov(3), v_cart_back(3) + real(dp) :: tol, err + integer :: i, n_passed + real(dp) :: test_points(3, 4) + real(dp) :: test_vectors(3, 3) + integer :: iv + + tol = 1.0e-10_dp + n_passed = 0 + + test_points(:, 1) = [0.25_dp, 1.0_dp, 0.5_dp] + test_points(:, 2) = [0.55_dp, 3.0_dp, 0.8_dp] + test_points(:, 3) = [0.45_dp, 5.5_dp, 0.3_dp] + test_points(:, 4) = [0.75_dp, 2.0_dp, 1.5_dp] + + test_vectors(:, 1) = [0.0_dp, 1.0_dp, 0.0_dp] + test_vectors(:, 2) = [0.8_dp, 0.2_dp, -0.5_dp] + test_vectors(:, 3) = [-0.4_dp, 0.9_dp, 1.2_dp] + + do i = 1, 4 + u = test_points(:, i) + call cs%covariant_basis(u, e_cov) + call cs%metric_tensor(u, g, ginv, sqrtg) + + do iv = 1, 3 + v_cart = test_vectors(:, iv) + call solve_3x3(e_cov, v_cart, v_ctr) + v_cov = matmul(g, v_ctr) + call cs%cov_to_cart(u, v_cov, v_cart_back) + + err = maxval(abs(v_cart_back - v_cart)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: cart->cov->cart roundtrip at u=", u, " err=", err + end if + end do + end do + + if (n_passed == 12) then + print *, " PASS: cart->cov->cart roundtrip (12 tests)" + else + print *, " FAIL: cart->cov->cart roundtrip", 12 - n_passed, "/12" + nerrors = nerrors + 1 + end if + end subroutine test_cart_cov_cart_roundtrip + + subroutine solve_3x3(A, b, x) + real(dp), intent(in) :: A(3, 3), b(3) + real(dp), intent(out) :: x(3) + + real(dp) :: det, Ainv(3, 3) + + det = A(1,1)*(A(2,2)*A(3,3) - A(2,3)*A(3,2)) & + - A(1,2)*(A(2,1)*A(3,3) - A(2,3)*A(3,1)) & + + A(1,3)*(A(2,1)*A(3,2) - A(2,2)*A(3,1)) + + Ainv(1,1) = (A(2,2)*A(3,3) - A(2,3)*A(3,2))/det + Ainv(1,2) = (A(1,3)*A(3,2) - A(1,2)*A(3,3))/det + Ainv(1,3) = (A(1,2)*A(2,3) - A(1,3)*A(2,2))/det + Ainv(2,1) = (A(2,3)*A(3,1) - A(2,1)*A(3,3))/det + Ainv(2,2) = (A(1,1)*A(3,3) - A(1,3)*A(3,1))/det + Ainv(2,3) = (A(1,3)*A(2,1) - A(1,1)*A(2,3))/det + Ainv(3,1) = (A(2,1)*A(3,2) - A(2,2)*A(3,1))/det + Ainv(3,2) = (A(1,2)*A(3,1) - A(1,1)*A(3,2))/det + Ainv(3,3) = (A(1,1)*A(2,2) - A(1,2)*A(2,1))/det + + x = matmul(Ainv, b) + end subroutine solve_3x3 + +end program test_vector_conversion diff --git a/test/source/test_vmec_coordinate_system.f90 b/test/source/test_vmec_coordinate_system.f90 index 10e9997d..01f69d67 100644 --- a/test/source/test_vmec_coordinate_system.f90 +++ b/test/source/test_vmec_coordinate_system.f90 @@ -1,15 +1,13 @@ program test_vmec_coordinate_system use, intrinsic :: iso_fortran_env, only: dp => real64 use libneo_coordinates, only: vmec_coordinate_system_t + use cylindrical_cartesian, only: cyl_to_cart use math_constants, only: TWOPI use new_vmec_stuff_mod, only: netcdffile, nper use spline_vmec_sub, only: spline_vmec_data implicit none type(vmec_coordinate_system_t) :: vmec - real(dp) :: u(3), xcyl(3) - real(dp) :: e_cov(3, 3) - real(dp) :: g(3, 3), ginv(3, 3), sqrtg real(dp) :: zeta_period integer :: nerrors @@ -23,43 +21,207 @@ program test_vmec_coordinate_system end if zeta_period = TWOPI/real(max(1, nper), dp) - u = [0.25_dp, 0.13_dp*TWOPI, 0.2_dp*zeta_period] - call vmec%evaluate_cyl(u, xcyl) - if (.not. (xcyl(1) > 100.0_dp .and. xcyl(1) < 2.0e5_dp)) then - print *, " FAIL: VMEC R not in cm range, R=", xcyl(1) - nerrors = nerrors + 1 - end if - if (abs(xcyl(2) - u(3)) > 1.0e-14_dp) then - print *, " FAIL: VMEC phi mismatch" - nerrors = nerrors + 1 - end if - if (.not. (abs(xcyl(3)) < 2.0e5_dp)) then - print *, " FAIL: VMEC Z not finite/plausible, Z=", xcyl(3) - nerrors = nerrors + 1 - end if - - call vmec%covariant_basis(u, e_cov) - if (any(.not. (abs(e_cov) < huge(1.0_dp)))) then - print *, " FAIL: VMEC covariant basis contains invalid values" - nerrors = nerrors + 1 - end if - - call vmec%metric_tensor(u, g, ginv, sqrtg) - if (.not. (sqrtg > 0.0_dp .and. sqrtg < huge(1.0_dp))) then - print *, " FAIL: VMEC sqrtg invalid, sqrtg=", sqrtg - nerrors = nerrors + 1 - end if - if (maxval(abs(g - transpose(g))) > 1.0e-10_dp) then - print *, " FAIL: VMEC metric tensor not symmetric" - nerrors = nerrors + 1 - end if + call test_basic_evaluate_cyl(vmec, zeta_period, nerrors) + call test_evaluate_cart_consistency(vmec, zeta_period, nerrors) + call test_covariant_basis_fd(vmec, zeta_period, nerrors) + call test_metric_tensor_properties(vmec, zeta_period, nerrors) if (nerrors > 0) then print *, "FAILED: ", nerrors, " error(s) in VMEC coordinate system test" error stop 1 end if - print *, " PASS: VMEC coordinate system basic checks" + print *, "All VMEC coordinate system tests passed!" + +contains + + subroutine test_basic_evaluate_cyl(vmec, zeta_period, nerrors) + type(vmec_coordinate_system_t), intent(in) :: vmec + real(dp), intent(in) :: zeta_period + integer, intent(inout) :: nerrors + + real(dp) :: u(3), xcyl(3) + + u = [0.25_dp, 0.13_dp*TWOPI, 0.2_dp*zeta_period] + call vmec%evaluate_cyl(u, xcyl) + + if (.not. (xcyl(1) > 100.0_dp .and. xcyl(1) < 2.0e5_dp)) then + print *, " FAIL: VMEC R not in cm range, R=", xcyl(1) + nerrors = nerrors + 1 + end if + if (abs(xcyl(2) - u(3)) > 1.0e-14_dp) then + print *, " FAIL: VMEC phi mismatch" + nerrors = nerrors + 1 + end if + if (.not. (abs(xcyl(3)) < 2.0e5_dp)) then + print *, " FAIL: VMEC Z not finite/plausible, Z=", xcyl(3) + nerrors = nerrors + 1 + else + print *, " PASS: evaluate_cyl basic sanity" + end if + end subroutine test_basic_evaluate_cyl + + subroutine test_evaluate_cart_consistency(vmec, zeta_period, nerrors) + type(vmec_coordinate_system_t), intent(in) :: vmec + real(dp), intent(in) :: zeta_period + integer, intent(inout) :: nerrors + + real(dp) :: u(3), xcyl(3), xcart_oop(3), xcart_manual(3) + real(dp), parameter :: tol = 1.0e-14_dp + real(dp) :: err + integer :: i, n_passed + real(dp) :: test_points(3, 4) + + n_passed = 0 + + test_points(:, 1) = [0.2_dp, 0.5_dp, 0.1_dp*zeta_period] + test_points(:, 2) = [0.5_dp, 2.0_dp, 0.3_dp*zeta_period] + test_points(:, 3) = [0.8_dp, 4.5_dp, 0.7_dp*zeta_period] + test_points(:, 4) = [0.3_dp, 1.0_dp, 0.0_dp] + + do i = 1, 4 + u = test_points(:, i) + + call vmec%evaluate_cart(u, xcart_oop) + call vmec%evaluate_cyl(u, xcyl) + call cyl_to_cart(xcyl, xcart_manual) + + err = maxval(abs(xcart_oop - xcart_manual)) + if (err < tol) then + n_passed = n_passed + 1 + else + print *, " FAIL: evaluate_cart inconsistent at u=", u + print *, " err=", err + end if + end do + + if (n_passed == 4) then + print *, " PASS: evaluate_cart = cyl_to_cart(evaluate_cyl)" + else + print *, " FAIL: evaluate_cart consistency", 4 - n_passed, "/4" + nerrors = nerrors + 1 + end if + end subroutine test_evaluate_cart_consistency + + subroutine test_covariant_basis_fd(vmec, zeta_period, nerrors) + type(vmec_coordinate_system_t), intent(in) :: vmec + real(dp), intent(in) :: zeta_period + integer, intent(inout) :: nerrors + + real(dp) :: u(3), e_cov(3, 3), e_cov_fd(3, 3) + real(dp) :: u_plus(3), u_minus(3), x_plus(3), x_minus(3) + real(dp), parameter :: h = 1.0e-6_dp + real(dp), parameter :: tol = 1.0e-4_dp + real(dp) :: rel_err, max_err + integer :: k, i, n_passed + real(dp) :: test_points(3, 4) + + n_passed = 0 + max_err = 0.0_dp + + test_points(:, 1) = [0.3_dp, 1.0_dp, 0.2_dp*zeta_period] + test_points(:, 2) = [0.5_dp, 3.0_dp, 0.5_dp*zeta_period] + test_points(:, 3) = [0.7_dp, 5.0_dp, 0.8_dp*zeta_period] + test_points(:, 4) = [0.4_dp, 0.5_dp, 0.0_dp] + + do i = 1, 4 + u = test_points(:, i) + + call vmec%covariant_basis(u, e_cov) + + do k = 1, 3 + u_plus = u + u_minus = u + u_plus(k) = u(k) + h + u_minus(k) = u(k) - h + call vmec%evaluate_cart(u_plus, x_plus) + call vmec%evaluate_cart(u_minus, x_minus) + e_cov_fd(:, k) = (x_plus - x_minus) / (2.0_dp * h) + end do + + rel_err = maxval(abs(e_cov - e_cov_fd)) / max(maxval(abs(e_cov_fd)), 1.0e-10_dp) + if (rel_err > max_err) max_err = rel_err + + if (rel_err < tol) then + n_passed = n_passed + 1 + else + print *, " covariant_basis FD mismatch at u=", u, " rel_err=", rel_err + end if + end do + + if (n_passed == 4) then + print *, " PASS: covariant_basis matches FD (4 points, max_err=", max_err, ")" + else + print *, " FAIL: covariant_basis FD", 4 - n_passed, "/4" + nerrors = nerrors + 1 + end if + end subroutine test_covariant_basis_fd + + subroutine test_metric_tensor_properties(vmec, zeta_period, nerrors) + type(vmec_coordinate_system_t), intent(in) :: vmec + real(dp), intent(in) :: zeta_period + integer, intent(inout) :: nerrors + + real(dp) :: u(3), g(3, 3), ginv(3, 3), sqrtg + real(dp) :: identity(3, 3), prod(3, 3) + real(dp), parameter :: tol = 1.0e-10_dp + integer :: i, j, k, ip, n_passed + logical :: sym_ok, inv_ok, pos_ok + real(dp) :: test_points(3, 3) + + n_passed = 0 + + test_points(:, 1) = [0.25_dp, 0.13_dp*TWOPI, 0.2_dp*zeta_period] + test_points(:, 2) = [0.6_dp, 2.5_dp, 0.4_dp*zeta_period] + test_points(:, 3) = [0.8_dp, 4.0_dp, 0.9_dp*zeta_period] + + do ip = 1, 3 + u = test_points(:, ip) + call vmec%metric_tensor(u, g, ginv, sqrtg) + + sym_ok = .true. + do i = 1, 3 + do j = i + 1, 3 + if (abs(g(i, j) - g(j, i)) > tol) sym_ok = .false. + if (abs(ginv(i, j) - ginv(j, i)) > tol) sym_ok = .false. + end do + end do + + identity = 0.0_dp + do i = 1, 3 + identity(i, i) = 1.0_dp + end do + + prod = 0.0_dp + do i = 1, 3 + do j = 1, 3 + do k = 1, 3 + prod(i, j) = prod(i, j) + g(i, k) * ginv(k, j) + end do + end do + end do + + inv_ok = maxval(abs(prod - identity)) < tol + pos_ok = sqrtg > 0.0_dp .and. sqrtg < huge(1.0_dp) + + if (sym_ok .and. inv_ok .and. pos_ok) then + n_passed = n_passed + 1 + else + print *, " metric_tensor fail at u=", u + if (.not. sym_ok) print *, " not symmetric" + if (.not. inv_ok) print *, " g*ginv /= I, err=", maxval(abs(prod - identity)) + if (.not. pos_ok) print *, " sqrtg invalid:", sqrtg + end if + end do + + if (n_passed == 3) then + print *, " PASS: metric_tensor properties (3 points)" + else + print *, " FAIL: metric_tensor properties", 3 - n_passed, "/3" + nerrors = nerrors + 1 + end if + end subroutine test_metric_tensor_properties + end program test_vmec_coordinate_system