diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4e7c145c..b71dafcf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -38,7 +38,13 @@ orbit_symplectic_base.f90 orbit_symplectic_quasi.f90 orbit_symplectic_euler1.f90 + orbit_rk_core.f90 orbit_symplectic.f90 + orbit_cpp.f90 + orbit_full_provider.f90 + orbit_full_mock_cart.f90 + orbit_full_mock_cyl.f90 + orbit_full.f90 util.F90 samplers.f90 cut_detector.f90 diff --git a/src/orbit_cpp.f90 b/src/orbit_cpp.f90 new file mode 100644 index 00000000..1ea8b04d --- /dev/null +++ b/src/orbit_cpp.f90 @@ -0,0 +1,175 @@ +module orbit_cpp + ! Classical Pauli Particle (CPP) pusher, flux-canonical realization + ! (Xiao & Qin, CPC 265 (2021) 107981). The guiding-center motion is the slow + ! manifold of the Pauli particle; a structure-preserving (variational / + ! Gauss-collocation) integrator stays on it and reproduces GC at GC-sized + ! (bounce-scale) steps. + ! + ! State and field machinery are shared with the GC integrator verbatim: + ! z(4) = (r, theta, phi, p_phi) in symplectic_integrator_t, + ! field_can_t carries mu (fixed parameter), ro0, and the canonical field + ! quantities Ath, Aph, hth, hph, Bmod with 1st and 2nd derivatives. + ! + ! The discrete scheme is the degenerate-Lagrangian Euler-Lagrange system + ! (implicit Gauss collocation) on field_can_t with mu held fixed. Because the + ! GC canonical equations ARE the slow-manifold equations of the Pauli + ! particle, the CPP Gauss residual coincides with the GC Gauss residual at + ! fixed mu; test_cpp_equals_gc_largestep verifies the trajectories agree to + ! Newton tolerance. + ! + ! GPU portability: residual, Jacobian, the dense LU solve, and the Newton + ! shell are pure, fixed-size, !$acc routine seq-able. No procedure pointers, + ! no class() polymorphism in the per-step loop; the only runtime indirection + ! is the field-evaluation pointer in field_can_mod (shared with GC and resolved + ! once at init). The stage count s (<= S_CPP_MAX) parameterizes GAUSS1..4. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: pi + use field_can_mod, only: field_can_t, get_derivatives2, eval_field => evaluate + use orbit_symplectic_base, only: symplectic_integrator_t, coeff_rk_gauss, & + GAUSS1, GAUSS2, GAUSS3, GAUSS4, extrap_field + use orbit_rk_core, only: rk_gauss_newton, MODEL_CPP + use diag_counters, only: count_event, EVT_RK_GAUSS_MAXIT, EVT_R_NEGATIVE + + implicit none + private + + integer, parameter, public :: S_CPP_MAX = 4 + + public :: orbit_cpp_init, orbit_timestep_cpp, cpp_stages_from_mode + +contains + + ! Map an integmode-style mode (GAUSS1..GAUSS4) to a Gauss stage count. + pure function cpp_stages_from_mode(mode) result(s) + integer, intent(in) :: mode + integer :: s + select case (mode) + case (GAUSS1) + s = 1 + case (GAUSS2) + s = 2 + case (GAUSS3) + s = 3 + case (GAUSS4) + s = 4 + case default + s = 1 + end select + end function cpp_stages_from_mode + + ! Initialize a CPP step from an already-prepared (si, f). The caller fills + ! si%z, si%dt, si%ntau, si%atol, si%rtol and f%mu, f%ro0 (e.g. via + ! simple%init_sympl, which sets mu from the initial pitch identically to GC). + ! This seeds f%pth at the current state so the first residual is consistent. + subroutine orbit_cpp_init(si, f) + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + + call eval_field(f, si%z(1), si%z(2), si%z(3), 0) + call get_derivatives2(f, si%z(4)) + end subroutine orbit_cpp_init + + ! One CPP macro-step (si%ntau micro-steps of si%dt). Advances si%z and f + ! exactly as orbit_timestep_sympl_rk_gauss does for GC, on the CPP residual. + recursive subroutine orbit_timestep_cpp(si, f, s, ierr) + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(in) :: s + integer, intent(out) :: ierr + + integer, parameter :: maxit = 32 + real(dp) :: x(4*s), xlast(4*s) + real(dp) :: a(s,s), b(s), c(s), Hprime(s), dz(4) + type(field_can_t) :: fs(s) + integer :: k, l, ktau + logical :: hit_maxit + + do k = 1, s + fs(k) = f + end do + + ierr = 0 + ktau = 0 + do while (ktau < si%ntau) + si%pthold = f%pth + + do k = 1, s + x((4*k-3):(4*k)) = si%z + end do + + call rk_gauss_newton(MODEL_CPP, si, fs, s, x, si%atol, si%rtol, maxit, & + xlast, hit_maxit) + if (hit_maxit) call count_event(EVT_RK_GAUSS_MAXIT) + + if (x(1) > 1.0d0) then + ierr = 1 + return + end if + + if (x(1) < 0.0d0) then + x(1) = -x(1) + x(2) = x(2) + pi + call count_event(EVT_R_NEGATIVE) + if (x(1) > 1.0d0) then + ierr = 1 + return + end if + end if + + call coeff_rk_gauss(s, a, b, c) + + if (extrap_field) then + do k = 1, s + dz(1) = x(4*k-3) - xlast(4*k-3) + dz(2) = x(4*k-2) - xlast(4*k-2) + dz(3) = x(4*k-1) - xlast(4*k-1) + dz(4) = x(4*k) - xlast(4*k) + + fs(k)%pth = fs(k)%pth + fs(k)%dpth(1)*dz(1) + fs(k)%dpth(2)*dz(2) & + + fs(k)%dpth(3)*dz(3) + fs(k)%dpth(4)*dz(4) + fs(k)%dpth(1) = fs(k)%dpth(1) + fs(k)%d2pth(1)*dz(1) + fs(k)%d2pth(2)*dz(2) & + + fs(k)%d2pth(3)*dz(3) + fs(k)%d2pth(7)*dz(4) + fs(k)%dpth(2) = fs(k)%dpth(2) + fs(k)%d2pth(2)*dz(1) + fs(k)%d2pth(4)*dz(2) & + + fs(k)%d2pth(5)*dz(3) + fs(k)%d2pth(8)*dz(4) + fs(k)%dpth(3) = fs(k)%dpth(3) + fs(k)%d2pth(2)*dz(1) + fs(k)%d2pth(5)*dz(2) & + + fs(k)%d2pth(6)*dz(3) + fs(k)%d2pth(9)*dz(4) + fs(k)%dH(1) = fs(k)%dH(1) + fs(k)%d2H(1)*dz(1) + fs(k)%d2H(2)*dz(2) & + + fs(k)%d2H(3)*dz(3) + fs(k)%d2H(7)*dz(4) + fs(k)%dH(2) = fs(k)%dH(2) + fs(k)%d2H(2)*dz(1) + fs(k)%d2H(4)*dz(2) & + + fs(k)%d2H(5)*dz(3) + fs(k)%d2H(8)*dz(4) + fs(k)%dH(3) = fs(k)%dH(3) + fs(k)%d2H(3)*dz(1) + fs(k)%d2H(5)*dz(2) & + + fs(k)%d2H(6)*dz(3) + fs(k)%d2H(9)*dz(4) + fs(k)%vpar = fs(k)%vpar + fs(k)%dvpar(1)*dz(1) + fs(k)%dvpar(2)*dz(2) & + + fs(k)%dvpar(3)*dz(3) + fs(k)%hth = fs(k)%hth + fs(k)%dhth(1)*dz(1) + fs(k)%dhth(2)*dz(2) & + + fs(k)%dhth(3)*dz(3) + fs(k)%hph = fs(k)%hph + fs(k)%dhph(1)*dz(1) + fs(k)%dhph(2)*dz(2) & + + fs(k)%dhph(3)*dz(3) + end do + else + do k = 1, s + call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 2) + call get_derivatives2(fs(k), x(4*k)) + end do + end if + + do k = 1, s + Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) + end do + + f = fs(s) + f%pth = si%pthold + si%z(1) = x(4*s-3) + + do l = 1, s + f%pth = f%pth - si%dt*b(l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) + si%z(2) = si%z(2) + si%dt*b(l)*Hprime(l) + si%z(3) = si%z(3) + si%dt*b(l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph + si%z(4) = si%z(4) - si%dt*b(l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) + end do + + ktau = ktau + 1 + end do + end subroutine orbit_timestep_cpp + +end module orbit_cpp diff --git a/src/orbit_full.f90 b/src/orbit_full.f90 new file mode 100644 index 00000000..b0e85b95 --- /dev/null +++ b/src/orbit_full.f90 @@ -0,0 +1,449 @@ +module orbit_full + ! Full-orbit (gyro-resolved Lorentz) pushers for SIMPLE, alongside the + ! existing symplectic guiding-center tracer. CGS Gaussian units throughout + ! (see src/util.F90): the Lorentz force carries the explicit 1/c, + ! m dv/dt = (q/c) v x B. + ! + ! ORBIT_BORIS is the explicit gyro-resolved pusher (this module). + ! ORBIT_FOSYMPL is the implicit-midpoint curvilinear full orbit (foimpl_step): + ! structure-preserving large-step counterpart of Boris, sharing the provider + ! seam, FullOrbitState, and helpers. + ! ORBIT_PAULI (CPP, 4D flux-canonical) lives in orbit_cpp, not here; it reuses + ! field_can_t and the symplectic GC state, dispatched from simple_main. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use orbit_full_provider, only: field_metric_provider_t, & + FO_OK, FO_ERR_FIELD, FO_ERR_NO_CONVERGE, FO_ERR_OUT_OF_DOMAIN + use orbit_rk_core, only: rk_solve + use util, only: c + implicit none + private + + ! orbit models (0 reserved for the existing symplectic guiding-center path) + integer, parameter, public :: ORBIT_GC = 0 + integer, parameter, public :: ORBIT_PAULI = 1 ! CPP variational, big dt, implicit + integer, parameter, public :: ORBIT_BORIS = 2 ! gyro-resolved Lorentz, explicit + integer, parameter, public :: ORBIT_FOSYMPL = 3 ! implicit-midpoint full orbit (B1) + + ! coordinate kinds (3..5 reserved for the libneo PR: VMEC, Boozer, chartmap) + integer, parameter, public :: COORD_CART = 1 + integer, parameter, public :: COORD_CYL = 2 + + ! error codes re-exported from the provider module + public :: FO_OK, FO_ERR_FIELD, FO_ERR_NO_CONVERGE, FO_ERR_OUT_OF_DOMAIN + + ! Name kept as FullOrbitState (no _t) for API compatibility with + ! test_full_orbit.f90. + type, public :: FullOrbitState + integer :: orbit_model = ORBIT_BORIS + integer :: ncoord = COORD_CART + ! Boris/Cartesian: z = (x1,x2,x3, v1,v2,v3) physical position and velocity. + ! Boris/curvilinear: z = (u1,u2,u3, v1,v2,v3) with v contravariant. + ! Pauli (CPP): z = (u1,u2,u3, vpar,0,0); z(1:4) carry state. + real(dp) :: z(6) = 0.0_dp + real(dp) :: mu = 0.0_dp ! m*vperp^2/(2 B), adiabatic invariant + real(dp) :: dt = 0.0_dp + real(dp) :: qm = 0.0_dp ! charge/mass (esu/g) + real(dp) :: mass = 0.0_dp ! particle mass (g) + real(dp) :: charge = 0.0_dp ! particle charge (esu) + class(field_metric_provider_t), pointer :: prov => null() + end type FullOrbitState + + public :: init_full_orbit_state, timestep_full_orbit, & + convert_full_to_gc, compute_energy + + interface init_full_orbit_state + module procedure init_full_orbit_state_prov + end interface init_full_orbit_state + +contains + + subroutine init_full_orbit_state_prov(state, x0, v0, orbit_model, ncoord, & + mass, charge, dt, prov) + type(FullOrbitState), intent(out) :: state + real(dp), intent(in) :: x0(3), v0(3) + integer, intent(in) :: orbit_model, ncoord + real(dp), intent(in) :: mass, charge, dt + class(field_metric_provider_t), intent(in), target :: prov + + state%orbit_model = orbit_model + state%ncoord = ncoord + state%mass = mass + state%charge = charge + state%qm = charge / mass + state%dt = dt + state%prov => prov + state%z(1:3) = x0 + state%z(4:6) = v0 + call set_mu_from_state(state) + end subroutine init_full_orbit_state_prov + + subroutine set_mu_from_state(state) + type(FullOrbitState), intent(inout) :: state + real(dp) :: Bvec(3), Bmod, hcov(3), vperp2, vpar + integer :: ierr + + call state%prov%eval_field(state%z(1:3), Bvec, Bmod, hcov, ierr) + if (ierr /= FO_OK .or. Bmod <= 0.0_dp) then + state%mu = 0.0_dp + return + end if + select case (state%ncoord) + case (COORD_CYL) + vperp2 = vperp_sq_cyl(state%z(1), state%z(4:6), hcov) + case default + vpar = dot_product(state%z(4:6), hcov) + vperp2 = max(dot_product(state%z(4:6), state%z(4:6)) - vpar*vpar, 0.0_dp) + end select + state%mu = state%mass * vperp2 / (2.0_dp * Bmod) + end subroutine set_mu_from_state + + subroutine timestep_full_orbit(state, ierr) + type(FullOrbitState), intent(inout) :: state + integer, intent(out) :: ierr + + select case (state%orbit_model) + case (ORBIT_BORIS) + select case (state%ncoord) + case (COORD_CART) + call boris_step_cart(state, ierr) + case (COORD_CYL) + call boris_step_cyl(state, ierr) + case default + ierr = FO_ERR_OUT_OF_DOMAIN + end select + case (ORBIT_FOSYMPL) + ! Implicit-midpoint curvilinear full orbit (VENUS-LEVIS geometry, B1): + ! m(dv^i/dt + Gamma^i_mn v^m v^n) = (q/c)(v x B)^i. The 6D residual is + ! fed to the shared Newton/LU core; structure-preserving large-step + ! counterpart of the explicit Boris pusher. + call foimpl_step(state, ierr) + case (ORBIT_PAULI) + ! CPP (4D flux-canonical) runs through the field_can / orbit_cpp path in + ! simple_main, not the 6D full-orbit state. Reaching it here means a + ! caller mis-routed a 4D model into the 6D pusher. + ierr = FO_ERR_OUT_OF_DOMAIN + case default + ierr = FO_ERR_OUT_OF_DOMAIN + end select + end subroutine timestep_full_orbit + + ! Contravariant velocity RHS for the curvilinear full orbit: + ! acc^i = -Gamma^i_{mn} v^m v^n + (q/(m c)) (v x B)^i. + ! The Lorentz term is built in the local orthonormal frame (mirroring the + ! Boris cyl rotation) and converted back to contravariant components. + subroutine fo_accel(state, u, v, acc, ierr) + type(FullOrbitState), intent(in) :: state + real(dp), intent(in) :: u(3), v(3) + real(dp), intent(out) :: acc(3) + integer, intent(out) :: ierr + real(dp) :: Bvec(3), Bmod, hcov(3), Gamma(3,3,3) + real(dp) :: vorth(3), forth(3), fcon(3), geo(3), qmc, R + integer :: i, m, n + + qmc = state%qm / c + call state%prov%eval_field(u, Bvec, Bmod, hcov, ierr) + if (ierr /= FO_OK) then + acc = 0.0_dp + return + end if + + call state%prov%christoffel(u, Gamma) + geo = 0.0_dp + do i = 1, 3 + do m = 1, 3 + do n = 1, 3 + geo(i) = geo(i) - Gamma(i,m,n) * v(m) * v(n) + end do + end do + end do + + select case (state%ncoord) + case (COORD_CYL) + R = u(1) + vorth = contravar_to_orth(v, R) + forth = qmc * cross(vorth, Bvec) + fcon = orth_to_contravar(forth, R) + case default + ! flat metric: orthonormal == contravariant + fcon = qmc * cross(v, Bvec) + end select + + acc = geo + fcon + end subroutine fo_accel + + ! Implicit-midpoint step for the 6D curvilinear full orbit. Residual + ! F = znew - zold - dt * rhs((zold+znew)/2), + ! rhs = (v, acc(u,v)). Solved by Newton with a finite-difference Jacobian and + ! the dense LU core; symplectic for the canonical lift, energy-bounded. + subroutine foimpl_step(state, ierr) + type(FullOrbitState), intent(inout) :: state + integer, intent(out) :: ierr + integer, parameter :: maxit = 32 + real(dp), parameter :: atol = 1.0e-13_dp, rtol = 1.0e-12_dp + real(dp) :: zold(6), z(6), fvec(6), fjac(6,6), zpert(6), fpert(6) + real(dp) :: dz(6), reltol(6), h + integer :: kit, i, j, info + logical :: conv + + zold = state%z + z = zold ! initial guess: explicit-Euler-free, start at old state + + do kit = 1, maxit + call foimpl_residual(state, zold, z, fvec, ierr) + if (ierr /= FO_OK) return + + do j = 1, 6 + h = max(1.0e-8_dp, 1.0e-7_dp*abs(z(j))) + zpert = z; zpert(j) = z(j) + h + call foimpl_residual(state, zold, zpert, fpert, ierr) + if (ierr /= FO_OK) return + do i = 1, 6 + fjac(i,j) = (fpert(i) - fvec(i)) / h + end do + end do + + dz = fvec + call rk_solve(6, fjac, dz, info) + if (info /= 0) then + ierr = FO_ERR_NO_CONVERGE + return + end if + z = z - dz + + do i = 1, 3 + reltol(i) = max(abs(z(i)), 1.0_dp) + reltol(i+3) = max(abs(z(i+3)), 1.0_dp) + end do + conv = .true. + do i = 1, 6 + if (abs(dz(i)) >= rtol*reltol(i) .and. abs(fvec(i)) >= atol) conv = .false. + end do + if (conv) exit + end do + + if (state%ncoord == COORD_CYL .and. z(1) <= 0.0_dp) then + ierr = FO_ERR_OUT_OF_DOMAIN + return + end if + + state%z = z + ierr = FO_OK + end subroutine foimpl_step + + subroutine foimpl_residual(state, zold, z, fvec, ierr) + type(FullOrbitState), intent(in) :: state + real(dp), intent(in) :: zold(6), z(6) + real(dp), intent(out) :: fvec(6) + integer, intent(out) :: ierr + real(dp) :: zmid(6), umid(3), vmid(3), acc(3), dt + + dt = state%dt + zmid = 0.5_dp * (zold + z) + umid = zmid(1:3) + vmid = zmid(4:6) + call fo_accel(state, umid, vmid, acc, ierr) + if (ierr /= FO_OK) then + fvec = 0.0_dp + return + end if + fvec(1:3) = z(1:3) - zold(1:3) - dt * vmid + fvec(4:6) = z(4:6) - zold(4:6) - dt * acc + end subroutine foimpl_residual + + ! Cartesian Boris, drift-kick-drift (leapfrog), CGS: + ! m dv/dt = (q/c) v x B, Omega = q B/(m c). + subroutine boris_step_cart(state, ierr) + type(FullOrbitState), intent(inout) :: state + integer, intent(out) :: ierr + real(dp) :: x(3), v(3), Bvec(3), Bmod, hcov(3) + real(dp) :: t(3), s(3), vprime(3), tmag2 + real(dp) :: Phi, dPhi(3), Efield(3) + real(dp) :: dt, qmc + + dt = state%dt + qmc = state%qm / c + x = state%z(1:3) + v = state%z(4:6) + + ! half position drift + x = x + 0.5_dp * dt * v + + call state%prov%eval_field(x, Bvec, Bmod, hcov, ierr) + if (ierr /= FO_OK) return + + ! half electric kick (E = -grad Phi); mocks return Phi=0 + call state%prov%eval_potential(x, Phi, dPhi) + Efield = -dPhi + v = v + 0.5_dp * dt * state%qm * Efield + + ! magnetic rotation (exact for constant B over the step) + t = qmc * Bvec * 0.5_dp * dt + tmag2 = dot_product(t, t) + s = 2.0_dp * t / (1.0_dp + tmag2) + vprime = v + cross(v, t) + v = v + cross(vprime, s) + + ! second half electric kick + v = v + 0.5_dp * dt * state%qm * Efield + + ! half position drift + x = x + 0.5_dp * dt * v + + state%z(1:3) = x + state%z(4:6) = v + ierr = FO_OK + end subroutine boris_step_cart + + ! Cylindrical Boris in coordinates (R,phi,Z), contravariant velocity in + ! z(4:6). Symmetric splitting: half geodesic kick / half drift / orthonormal + ! magnetic rotation / half drift / half geodesic kick. The geodesic kick + ! applies -Gamma^i_{mn} v^m v^n from the provider's Christoffel symbols. + subroutine boris_step_cyl(state, ierr) + type(FullOrbitState), intent(inout) :: state + integer, intent(out) :: ierr + real(dp) :: u(3), v(3), Bvec(3), Bmod, hcov(3) + real(dp) :: dt, qmc, R + real(dp) :: vorth(3), t(3), sv(3), vprime(3), tmag2 + + dt = state%dt + qmc = state%qm / c + u = state%z(1:3) + v = state%z(4:6) + + call geodesic_half_kick(state, u, v, 0.5_dp*dt) + + u = u + 0.5_dp * dt * v + R = u(1) + if (R <= 0.0_dp) then + ierr = FO_ERR_OUT_OF_DOMAIN + return + end if + + call state%prov%eval_field(u, Bvec, Bmod, hcov, ierr) + if (ierr /= FO_OK) return + + ! magnetic rotation in the local orthonormal frame + vorth = contravar_to_orth(v, R) + t = qmc * Bvec * 0.5_dp * dt + tmag2 = dot_product(t, t) + sv = 2.0_dp * t / (1.0_dp + tmag2) + vprime = vorth + cross(vorth, t) + vorth = vorth + cross(vprime, sv) + v = orth_to_contravar(vorth, R) + + u = u + 0.5_dp * dt * v + if (u(1) <= 0.0_dp) then + ierr = FO_ERR_OUT_OF_DOMAIN + return + end if + + call geodesic_half_kick(state, u, v, 0.5_dp*dt) + + state%z(1:3) = u + state%z(4:6) = v + ierr = FO_OK + end subroutine boris_step_cyl + + subroutine geodesic_half_kick(state, u, v, h) + type(FullOrbitState), intent(in) :: state + real(dp), intent(in) :: u(3), h + real(dp), intent(inout) :: v(3) + real(dp) :: Gamma(3,3,3), acc(3) + integer :: i, m, n + + call state%prov%christoffel(u, Gamma) + acc = 0.0_dp + do i = 1, 3 + do m = 1, 3 + do n = 1, 3 + acc(i) = acc(i) - Gamma(i,m,n) * v(m) * v(n) + end do + end do + end do + v = v + h * acc + end subroutine geodesic_half_kick + + ! Orthonormal cylindrical components (vR,vphi_phys,vZ) from contravariant + ! (v^R, v^phi, v^Z): vphi_phys = R v^phi. + pure function contravar_to_orth(v, R) result(vorth) + real(dp), intent(in) :: v(3), R + real(dp) :: vorth(3) + vorth = [v(1), R * v(2), v(3)] + end function contravar_to_orth + + pure function orth_to_contravar(vorth, R) result(v) + real(dp), intent(in) :: vorth(3), R + real(dp) :: v(3) + v = [vorth(1), vorth(2) / R, vorth(3)] + end function orth_to_contravar + + pure function vperp_sq_cyl(R, v, hcov) result(vperp2) + real(dp), intent(in) :: R, v(3), hcov(3) + real(dp) :: vperp2, vorth(3), hunit(3), vpar, vsq + vorth = contravar_to_orth(v, R) + ! hcov are covariant unit-field comps; orthonormal unit field is hcov/g_ii. + ! For the toroidal mock h_orth = (0,1,0). + hunit = [hcov(1), hcov(2) / R, hcov(3)] + vpar = dot_product(vorth, hunit) + vsq = dot_product(vorth, vorth) + vperp2 = max(vsq - vpar*vpar, 0.0_dp) + end function vperp_sq_cyl + + pure function cross(a, b) result(c) + real(dp), intent(in) :: a(3), b(3) + real(dp) :: c(3) + c(1) = a(2)*b(3) - a(3)*b(2) + c(2) = a(3)*b(1) - a(1)*b(3) + c(3) = a(1)*b(2) - a(2)*b(1) + end function cross + + ! Gyro-average to guiding-center scalars. For the analytic mocks this returns + ! position-coordinate scalars and the instantaneous pitch/speed; full + ! gyro-averaging over coordinate maps lands with the libneo PR. + subroutine convert_full_to_gc(state, s, theta, phi, lambda, v) + type(FullOrbitState), intent(in) :: state + real(dp), intent(out) :: s, theta, phi, lambda, v + real(dp) :: Bvec(3), Bmod, hcov(3), vpar, speed + real(dp) :: vorth(3), hunit(3) + integer :: ierr + + call state%prov%eval_field(state%z(1:3), Bvec, Bmod, hcov, ierr) + select case (state%ncoord) + case (COORD_CYL) + s = state%z(1); theta = state%z(3); phi = state%z(2) + vorth = contravar_to_orth(state%z(4:6), state%z(1)) + hunit = [hcov(1), hcov(2)/state%z(1), hcov(3)] + speed = sqrt(dot_product(vorth, vorth)) + vpar = dot_product(vorth, hunit) + case default + s = state%z(1); theta = state%z(2); phi = state%z(3) + speed = sqrt(dot_product(state%z(4:6), state%z(4:6))) + vpar = dot_product(state%z(4:6), hcov) + end select + v = speed + if (speed > 0.0_dp) then + lambda = vpar / speed + else + lambda = 0.0_dp + end if + end subroutine convert_full_to_gc + + ! 0.5*m*|v|^2 (+ q*Phi). Velocity magnitude uses the metric in curvilinear + ! coordinates so it is the physical speed. + function compute_energy(state) result(energy) + type(FullOrbitState), intent(in) :: state + real(dp) :: energy + real(dp) :: speed2, Phi, dPhi(3), vorth(3) + + select case (state%ncoord) + case (COORD_CYL) + vorth = contravar_to_orth(state%z(4:6), state%z(1)) + speed2 = dot_product(vorth, vorth) + case default + speed2 = dot_product(state%z(4:6), state%z(4:6)) + end select + call state%prov%eval_potential(state%z(1:3), Phi, dPhi) + energy = 0.5_dp * state%mass * speed2 + state%charge * Phi + end function compute_energy + +end module orbit_full diff --git a/src/orbit_full_mock_cart.f90 b/src/orbit_full_mock_cart.f90 new file mode 100644 index 00000000..04c70010 --- /dev/null +++ b/src/orbit_full_mock_cart.f90 @@ -0,0 +1,120 @@ +module orbit_full_mock_cart + ! Cartesian mock field/metric provider, no libneo dependency on the analytic + ! paths. Flat metric (g = I), zero Christoffel. B is either uniform, a linear + ! gradient, or (FIELD_COILS) a Biot-Savart evaluation of a coil set. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use orbit_full_provider, only: field_metric_provider_t, FO_ERR_FIELD + use neo_biotsavart, only: coils_t, compute_magnetic_field, compute_vector_potential + implicit none + private + public :: cartesian_provider_t + public :: FIELD_UNIFORM, FIELD_LINGRAD, FIELD_COILS + + integer, parameter :: FIELD_UNIFORM = 1 + integer, parameter :: FIELD_LINGRAD = 2 + integer, parameter :: FIELD_COILS = 3 + + type, extends(field_metric_provider_t), public :: cartesian_provider_t + integer :: field_kind = FIELD_UNIFORM + real(dp) :: B0(3) = [0.0_dp, 0.0_dp, 1.0_dp] + real(dp) :: gradB(3,3) = 0.0_dp ! B_i(x) = B0_i + sum_j gradB(i,j) x_j + type(coils_t), pointer :: coils => null() + contains + procedure :: eval_field => cart_eval_field + procedure :: metric => cart_metric + procedure :: christoffel => cart_christoffel + procedure :: eval_canfield => cart_eval_canfield + procedure :: eval_potential => cart_eval_potential + end type cartesian_provider_t + +contains + + subroutine cart_eval_field(self, x, Bvec, Bmod, hcov, ierr) + class(cartesian_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Bvec(3), Bmod, hcov(3) + integer, intent(out) :: ierr + integer :: i + + ierr = 0 + select case (self%field_kind) + case (FIELD_UNIFORM) + Bvec = self%B0 + case (FIELD_LINGRAD) + do i = 1, 3 + Bvec(i) = self%B0(i) + dot_product(self%gradB(i, :), x) + end do + case (FIELD_COILS) + if (.not. associated(self%coils)) then + ierr = FO_ERR_FIELD + Bvec = 0.0_dp; Bmod = 0.0_dp; hcov = 0.0_dp + return + end if + Bvec = compute_magnetic_field(self%coils, x) + case default + ierr = FO_ERR_FIELD + Bvec = 0.0_dp; Bmod = 0.0_dp; hcov = 0.0_dp + return + end select + + Bmod = sqrt(Bvec(1)**2 + Bvec(2)**2 + Bvec(3)**2) + if (Bmod <= 0.0_dp) then + ierr = FO_ERR_FIELD + hcov = 0.0_dp + return + end if + hcov = Bvec / Bmod + end subroutine cart_eval_field + + subroutine cart_metric(self, x, g, ginv, sqrtg) + class(cartesian_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg + integer :: i + + g = 0.0_dp + ginv = 0.0_dp + do i = 1, 3 + g(i,i) = 1.0_dp + ginv(i,i) = 1.0_dp + end do + sqrtg = 1.0_dp + end subroutine cart_metric + + subroutine cart_christoffel(self, x, Gamma) + class(cartesian_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Gamma(3,3,3) + + Gamma = 0.0_dp + end subroutine cart_christoffel + + ! Declared seam for the CPP/Pauli path; the Boris path never calls it. The + ! mock has no analytic canonical-field model, so it flags not-implemented. + subroutine cart_eval_canfield(self, x, Ath, Aph, hth, hph, Bmod, & + dAth, dAph, dhth, dhph, dBmod, & + d2Ath, d2Aph, d2hth, d2hph, d2Bmod, ierr) + class(cartesian_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Ath, Aph, hth, hph, Bmod + real(dp), intent(out) :: dAth(3), dAph(3), dhth(3), dhph(3), dBmod(3) + real(dp), intent(out) :: d2Ath(6), d2Aph(6), d2hth(6), d2hph(6), d2Bmod(6) + integer, intent(out) :: ierr + + Ath = 0.0_dp; Aph = 0.0_dp; hth = 0.0_dp; hph = 0.0_dp; Bmod = 0.0_dp + dAth = 0.0_dp; dAph = 0.0_dp; dhth = 0.0_dp; dhph = 0.0_dp; dBmod = 0.0_dp + d2Ath = 0.0_dp; d2Aph = 0.0_dp; d2hth = 0.0_dp; d2hph = 0.0_dp + d2Bmod = 0.0_dp + ierr = FO_ERR_FIELD + end subroutine cart_eval_canfield + + subroutine cart_eval_potential(self, x, Phi, dPhi) + class(cartesian_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Phi, dPhi(3) + + Phi = 0.0_dp + dPhi = 0.0_dp + end subroutine cart_eval_potential + +end module orbit_full_mock_cart diff --git a/src/orbit_full_mock_cyl.f90 b/src/orbit_full_mock_cyl.f90 new file mode 100644 index 00000000..b4d3f779 --- /dev/null +++ b/src/orbit_full_mock_cyl.f90 @@ -0,0 +1,112 @@ +module orbit_full_mock_cyl + ! Cylindrical mock provider in coordinates u = (R, phi, Z), no libneo + ! dependency. Analytic metric and closed-form Christoffel symbols, plus an + ! axisymmetric purely-toroidal 1/R field for the curvature + grad-B drift + ! oracle. + ! + ! g = diag(1, R^2, 1), ginv = diag(1, 1/R^2, 1), sqrtg = R + ! B = B0*R0/R in the orthonormal toroidal direction (|B| = B0*R0/R) + ! B_phi^cov = |B|*R = B0*R0 (covariant toroidal component) + ! + ! Nonzero Christoffel (second kind), Gamma(i,m,n) = Gamma^i_{mn}: + ! Gamma^R_{phi phi} = Gamma(1,2,2) = -R + ! Gamma^phi_{R phi} = Gamma(2,1,2) = 1/R + ! Gamma^phi_{phi R} = Gamma(2,2,1) = 1/R + use, intrinsic :: iso_fortran_env, only: dp => real64 + use orbit_full_provider, only: field_metric_provider_t, FO_ERR_FIELD + implicit none + private + public :: cylindrical_provider_t + + type, extends(field_metric_provider_t), public :: cylindrical_provider_t + real(dp) :: B0 = 1.0_dp + real(dp) :: R0 = 1.0_dp + contains + procedure :: eval_field => cyl_eval_field + procedure :: metric => cyl_metric + procedure :: christoffel => cyl_christoffel + procedure :: eval_canfield => cyl_eval_canfield + procedure :: eval_potential => cyl_eval_potential + end type cylindrical_provider_t + +contains + + ! Bvec returned in the orthonormal cylindrical frame (e_R, e_phi, e_Z); + ! hcov holds the covariant unit-field components h_i = B_i/|B|, so that + ! h_phi = R for this purely-toroidal field. + subroutine cyl_eval_field(self, x, Bvec, Bmod, hcov, ierr) + class(cylindrical_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Bvec(3), Bmod, hcov(3) + integer, intent(out) :: ierr + real(dp) :: R + + R = x(1) + if (R <= 0.0_dp) then + ierr = FO_ERR_FIELD + Bvec = 0.0_dp; Bmod = 0.0_dp; hcov = 0.0_dp + return + end if + ierr = 0 + Bmod = self%B0 * self%R0 / R + Bvec = [0.0_dp, Bmod, 0.0_dp] ! orthonormal toroidal + hcov = [0.0_dp, R, 0.0_dp] ! h_phi^cov = B_phi^cov/|B| = R + end subroutine cyl_eval_field + + subroutine cyl_metric(self, x, g, ginv, sqrtg) + class(cylindrical_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg + real(dp) :: R + + R = x(1) + g = 0.0_dp + ginv = 0.0_dp + g(1,1) = 1.0_dp; ginv(1,1) = 1.0_dp + g(2,2) = R*R; ginv(2,2) = 1.0_dp / (R*R) + g(3,3) = 1.0_dp; ginv(3,3) = 1.0_dp + sqrtg = R + end subroutine cyl_metric + + subroutine cyl_christoffel(self, x, Gamma) + class(cylindrical_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Gamma(3,3,3) + real(dp) :: R + + R = x(1) + Gamma = 0.0_dp + Gamma(1,2,2) = -R + Gamma(2,1,2) = 1.0_dp / R + Gamma(2,2,1) = 1.0_dp / R + end subroutine cyl_christoffel + + ! Declared CPP/Pauli seam; not exercised by the Boris path. Flags + ! not-implemented rather than returning a partial canonical model. + subroutine cyl_eval_canfield(self, x, Ath, Aph, hth, hph, Bmod, & + dAth, dAph, dhth, dhph, dBmod, & + d2Ath, d2Aph, d2hth, d2hph, d2Bmod, ierr) + class(cylindrical_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Ath, Aph, hth, hph, Bmod + real(dp), intent(out) :: dAth(3), dAph(3), dhth(3), dhph(3), dBmod(3) + real(dp), intent(out) :: d2Ath(6), d2Aph(6), d2hth(6), d2hph(6), d2Bmod(6) + integer, intent(out) :: ierr + + Ath = 0.0_dp; Aph = 0.0_dp; hth = 0.0_dp; hph = 0.0_dp; Bmod = 0.0_dp + dAth = 0.0_dp; dAph = 0.0_dp; dhth = 0.0_dp; dhph = 0.0_dp; dBmod = 0.0_dp + d2Ath = 0.0_dp; d2Aph = 0.0_dp; d2hth = 0.0_dp; d2hph = 0.0_dp + d2Bmod = 0.0_dp + ierr = FO_ERR_FIELD + end subroutine cyl_eval_canfield + + subroutine cyl_eval_potential(self, x, Phi, dPhi) + class(cylindrical_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Phi, dPhi(3) + + Phi = 0.0_dp + dPhi = 0.0_dp + end subroutine cyl_eval_potential + +end module orbit_full_mock_cyl diff --git a/src/orbit_full_provider.f90 b/src/orbit_full_provider.f90 new file mode 100644 index 00000000..8244f697 --- /dev/null +++ b/src/orbit_full_provider.f90 @@ -0,0 +1,77 @@ +module orbit_full_provider + ! Abstract field/metric provider for the full-orbit pushers. + ! + ! Kept in its own module so the Boris path carries no libneo symbol: the + ! mocks extend this type, orbit_full only depends on the abstract interface. + use, intrinsic :: iso_fortran_env, only: dp => real64 + implicit none + private + public :: field_metric_provider_t + public :: FO_OK, FO_ERR_FIELD, FO_ERR_NO_CONVERGE, FO_ERR_OUT_OF_DOMAIN + + ! Error codes shared by the provider seam and the integrator. Defined here so + ! the mocks can flag a not-implemented eval_canfield without depending on + ! orbit_full (which would create a use cycle); orbit_full re-exports them. + integer, parameter :: FO_OK = 0 + integer, parameter :: FO_ERR_FIELD = 1 ! provider field eval failed + integer, parameter :: FO_ERR_NO_CONVERGE = 2 ! CPP Newton did not converge + integer, parameter :: FO_ERR_OUT_OF_DOMAIN = 3 + + type, abstract :: field_metric_provider_t + contains + procedure(eval_field_if), deferred :: eval_field + procedure(metric_if), deferred :: metric + procedure(christoffel_if), deferred :: christoffel + procedure(eval_canfield_if), deferred :: eval_canfield + procedure(eval_potential_if), deferred :: eval_potential + end type field_metric_provider_t + + abstract interface + ! Physical B vector, |B|, and covariant unit-vector components h_i = B_i/|B|. + subroutine eval_field_if(self, x, Bvec, Bmod, hcov, ierr) + import :: field_metric_provider_t, dp + class(field_metric_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Bvec(3), Bmod, hcov(3) + integer, intent(out) :: ierr + end subroutine eval_field_if + + ! Metric for the curvilinear push: g_ij, g^ij, sqrt(det g). + subroutine metric_if(self, x, g, ginv, sqrtg) + import :: field_metric_provider_t, dp + class(field_metric_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: g(3,3), ginv(3,3), sqrtg + end subroutine metric_if + + ! Christoffel symbols of the second kind, Gamma(i,m,n) = Gamma^i_{mn}. + subroutine christoffel_if(self, x, Gamma) + import :: field_metric_provider_t, dp + class(field_metric_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Gamma(3,3,3) + end subroutine christoffel_if + + ! CPP / Pauli canonical field quantities, all in provider coordinates. + ! Second-derivative layout d2?(6): order (11,12,13,22,23,33). + subroutine eval_canfield_if(self, x, Ath, Aph, hth, hph, Bmod, & + dAth, dAph, dhth, dhph, dBmod, & + d2Ath, d2Aph, d2hth, d2hph, d2Bmod, ierr) + import :: field_metric_provider_t, dp + class(field_metric_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Ath, Aph, hth, hph, Bmod + real(dp), intent(out) :: dAth(3), dAph(3), dhth(3), dhph(3), dBmod(3) + real(dp), intent(out) :: d2Ath(6), d2Aph(6), d2hth(6), d2hph(6), d2Bmod(6) + integer, intent(out) :: ierr + end subroutine eval_canfield_if + + ! Electrostatic potential Phi and covariant gradient. Mocks return 0. + subroutine eval_potential_if(self, x, Phi, dPhi) + import :: field_metric_provider_t, dp + class(field_metric_provider_t), intent(in) :: self + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: Phi, dPhi(3) + end subroutine eval_potential_if + end interface +end module orbit_full_provider diff --git a/src/orbit_rk_core.f90 b/src/orbit_rk_core.f90 new file mode 100644 index 00000000..8c9f1d25 --- /dev/null +++ b/src/orbit_rk_core.f90 @@ -0,0 +1,299 @@ +module orbit_rk_core + ! Device-callable shared core for the implicit Gauss-collocation step used by + ! the guiding-center (GC) integrator and the Classical Pauli Particle (CPP) + ! pusher. The CPP guiding-center motion is the slow manifold of the Pauli + ! particle, so its degenerate-Lagrangian Gauss residual on field_can_t (mu + ! fixed) coincides with the GC Gauss residual; both map to one concrete + ! kernel here, selected by an integer model code. No procedure pointers and + ! no class() polymorphism in the hot path: dispatch is select case to + ! inlinable !$acc routine seq kernels, the way orbit_symplectic_euler1 + ! already shares the symplectic-Euler residual/Jacobian/Newton between the + ! CPU integrator and the OpenACC GPU kernel. + ! + ! The arithmetic is lifted byte-identically from the former f_rk_gauss / + ! jac_rk_gauss bodies in orbit_symplectic and the matching orbit_cpp bodies. + ! coeff_rk_gauss stays in orbit_symplectic_base (plain data) and is called + ! inside the kernels. The Newton shell here uses the device LU solve rk_solve; + ! the GC CPU path keeps dgesv in orbit_symplectic for byte-identical results. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: twopi + use field_can_mod, only: field_can_t, get_derivatives2, eval_field => evaluate + use orbit_symplectic_base, only: symplectic_integrator_t, coeff_rk_gauss + + implicit none + private + + ! Integer model codes for select-case dispatch (no proc pointers / class()). + integer, parameter, public :: MODEL_GC = 0, MODEL_CPP = 1, MODEL_FO_SYMPL = 2 + + public :: rk_gauss_residual, rk_gauss_jacobian, rk_gauss_newton, rk_solve + +contains + + ! Model-dispatched Gauss residual. Both MODEL_GC and MODEL_CPP evaluate the + ! field-canonical degenerate-Lagrangian residual (the CPP slow manifold equals + ! the GC equations at fixed mu), so they share one inlinable kernel. + subroutine rk_gauss_residual(model, si, fs, s, x, fvec) + !$acc routine seq + integer, intent(in) :: model + type(symplectic_integrator_t), intent(inout) :: si + integer, intent(in) :: s + type(field_can_t), intent(inout) :: fs(:) + real(dp), intent(in) :: x(4*s) + real(dp), intent(out) :: fvec(4*s) + + select case (model) + case (MODEL_GC, MODEL_CPP) + call gauss_canfield_residual(si, fs, s, x, fvec) + case default + fvec = 0d0 + end select + end subroutine rk_gauss_residual + + ! Model-dispatched Gauss Jacobian, analytic from the 2nd derivatives of + ! field_can_t. Shared kernel for MODEL_GC and MODEL_CPP. + subroutine rk_gauss_jacobian(model, si, fs, s, jac) + !$acc routine seq + integer, intent(in) :: model + type(symplectic_integrator_t), intent(in) :: si + integer, intent(in) :: s + type(field_can_t), intent(in) :: fs(:) + real(dp), intent(out) :: jac(4*s, 4*s) + + select case (model) + case (MODEL_GC, MODEL_CPP) + call gauss_canfield_jacobian(si, fs, s, jac) + case default + jac = 0d0 + end select + end subroutine rk_gauss_jacobian + + ! Concrete residual kernel: layout x(4*k-3:4*k) = (r,theta,phi,pphi) per stage. + ! Bodies lifted verbatim from f_rk_gauss / f_rk_gauss_cpp. + subroutine gauss_canfield_residual(si, fs, s, x, fvec) + !$acc routine seq + type(symplectic_integrator_t), intent(inout) :: si + integer, intent(in) :: s + type(field_can_t), intent(inout) :: fs(:) + real(dp), intent(in) :: x(4*s) + real(dp), intent(out) :: fvec(4*s) + + real(dp) :: a(s,s), b(s), c(s), Hprime(s) + integer :: k, l + + call coeff_rk_gauss(s, a, b, c) + + do k = 1, s + call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 2) + call get_derivatives2(fs(k), x(4*k)) + Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) + end do + + do k = 1, s + fvec(4*k-3) = fs(k)%pth - si%pthold + fvec(4*k-2) = x(4*k-2) - si%z(2) + fvec(4*k-1) = x(4*k-1) - si%z(3) + fvec(4*k) = x(4*k) - si%z(4) + do l = 1, s + fvec(4*k-3) = fvec(4*k-3) + si%dt*a(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot + fvec(4*k-2) = fvec(4*k-2) - si%dt*a(k,l)*Hprime(l) ! thdot + fvec(4*k-1) = fvec(4*k-1) - si%dt*a(k,l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph ! phdot + fvec(4*k) = fvec(4*k) + si%dt*a(k,l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot + end do + end do + end subroutine gauss_canfield_residual + + ! Concrete Jacobian kernel. Bodies lifted verbatim from jac_rk_gauss / + ! jac_rk_gauss_cpp. + subroutine gauss_canfield_jacobian(si, fs, s, jac) + !$acc routine seq + type(symplectic_integrator_t), intent(in) :: si + integer, intent(in) :: s + type(field_can_t), intent(in) :: fs(:) + real(dp), intent(out) :: jac(4*s, 4*s) + + real(dp) :: a(s,s), b(s), c(s), Hprime(s), dHprime(4*s) + integer :: k, l, m + + call coeff_rk_gauss(s, a, b, c) + + do k = 1, s + m = 4*k + Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) + dHprime(m-3) = (fs(k)%d2H(1) - Hprime(k)*fs(k)%d2pth(1))/fs(k)%dpth(1) ! d/dr + dHprime(m-2) = (fs(k)%d2H(2) - Hprime(k)*fs(k)%d2pth(2))/fs(k)%dpth(1) ! d/dth + dHprime(m-1) = (fs(k)%d2H(3) - Hprime(k)*fs(k)%d2pth(3))/fs(k)%dpth(1) ! d/dph + dHprime(m) = (fs(k)%d2H(7) - Hprime(k)*fs(k)%d2pth(7))/fs(k)%dpth(1) ! d/dpph + end do + + jac = 0d0 + + do k = 1, s + m = 4*k + jac(m-3, m-3) = fs(k)%dpth(1) + jac(m-3, m-2) = fs(k)%dpth(2) + jac(m-3, m-1) = fs(k)%dpth(3) + jac(m-3, m) = fs(k)%dpth(4) + jac(m-2, m-2) = 1d0 + jac(m-1, m-1) = 1d0 + jac(m, m) = 1d0 + end do + + do l = 1, s + do k = 1, s + m = 4*k + jac(m-3, 4*l-3) = jac(m-3, 4*l-3) & ! d/dr + + si%dt*a(k,l)*(fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-3)) + jac(m-3, 4*l-2) = jac(m-3, 4*l-2) & ! d/dth + + si%dt*a(k,l)*(fs(l)%d2H(4) - fs(l)%d2pth(4)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-2)) + jac(m-3, 4*l-1) = jac(m-3, 4*l-1) & ! d/dph + + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-1)) + jac(m-3, 4*l) = jac(m-3, 4*l) & ! d/dpph + + si%dt*a(k,l)*(fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l)) + + jac(m-2, 4*l-3) = jac(m-2, 4*l-3) - si%dt*a(k,l)*dHprime(4*l-3) ! d/dr + jac(m-2, 4*l-2) = jac(m-2, 4*l-2) - si%dt*a(k,l)*dHprime(4*l-2) ! d/dth + jac(m-2, 4*l-1) = jac(m-2, 4*l-1) - si%dt*a(k,l)*dHprime(4*l-1) ! d/dph + jac(m-2, 4*l) = jac(m-2, 4*l) - si%dt*a(k,l)*dHprime(4*l) ! d/dpph + + jac(m-1, 4*l-3) = jac(m-1, 4*l-3) & ! d/dr + - si%dt*a(k,l)*(-fs(l)%dhph(1)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & + + (fs(l)%dvpar(1) - dHprime(4*l-3)*fs(l)%hth - Hprime(l)*fs(l)%dhth(1))/fs(l)%hph) + jac(m-1, 4*l-2) = jac(m-1, 4*l-2) & ! d/dth + - si%dt*a(k,l)*(-fs(l)%dhph(2)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & + + (fs(l)%dvpar(2) - dHprime(4*l-2)*fs(l)%hth - Hprime(l)*fs(l)%dhth(2))/fs(l)%hph) + jac(m-1, 4*l-1) = jac(m-1, 4*l-1) & ! d/dph + - si%dt*a(k,l)*(-fs(l)%dhph(3)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & + + (fs(l)%dvpar(3) - dHprime(4*l-1)*fs(l)%hth - Hprime(l)*fs(l)%dhth(3))/fs(l)%hph) + jac(m-1, 4*l) = jac(m-1, 4*l) & ! d/dpph + - si%dt*a(k,l)*((fs(l)%dvpar(4) - dHprime(4*l)*fs(l)%hth)/fs(l)%hph) + + jac(m, 4*l-3) = jac(m, 4*l-3) & ! d/dr + + si%dt*a(k,l)*(fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-3)) + jac(m, 4*l-2) = jac(m, 4*l-2) & ! d/dth + + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-2)) + jac(m, 4*l-1) = jac(m, 4*l-1) & ! d/dph + + si%dt*a(k,l)*(fs(l)%d2H(6) - fs(l)%d2pth(6)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-1)) + jac(m, 4*l) = jac(m, 4*l) & ! d/dpph + + si%dt*a(k,l)*(fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l)) + end do + end do + end subroutine gauss_canfield_jacobian + + ! Device-portable dense LU solve A x = rhs with partial pivoting, in place on + ! rhs. n <= 4*S_MAX (n <= 16 for s <= 4). info = 0 on success, else the failing + ! pivot column. Replaces dgesv on the device; the GC CPU path keeps dgesv. + pure subroutine rk_solve(n, A, rhs, info) + !$acc routine seq + integer, intent(in) :: n + real(dp), intent(inout) :: A(n,n), rhs(n) + integer, intent(out) :: info + integer :: i, j, k, ipiv + real(dp) :: piv, amax, factor, tmp + + info = 0 + do k = 1, n + ipiv = k + amax = abs(A(k,k)) + do i = k+1, n + if (abs(A(i,k)) > amax) then + amax = abs(A(i,k)) + ipiv = i + end if + end do + if (amax == 0d0) then + info = k + return + end if + if (ipiv /= k) then + do j = 1, n + tmp = A(k,j); A(k,j) = A(ipiv,j); A(ipiv,j) = tmp + end do + tmp = rhs(k); rhs(k) = rhs(ipiv); rhs(ipiv) = tmp + end if + piv = A(k,k) + do i = k+1, n + factor = A(i,k)/piv + A(i,k) = factor + do j = k+1, n + A(i,j) = A(i,j) - factor*A(k,j) + end do + rhs(i) = rhs(i) - factor*rhs(k) + end do + end do + + do i = n, 1, -1 + tmp = rhs(i) + do j = i+1, n + tmp = tmp - A(i,j)*rhs(j) + end do + rhs(i) = tmp/A(i,i) + end do + end subroutine rk_solve + + ! Device-callable Newton iteration for the Gauss step. Mirrors the + ! newton_rk_gauss control flow (atol/rtol/tolref, boundary guards, maxit) with + ! the device LU solver rk_solve in place of dgesv. No event counters here: + ! hit_maxit is returned so the CPU caller can record EVT_RK_GAUSS_MAXIT, + ! exactly as the GPU newton1 path omits fort.6601. + subroutine rk_gauss_newton(model, si, fs, s, x, atol, rtol, maxit, xlast, & + hit_maxit) + !$acc routine seq + integer, intent(in) :: model + type(symplectic_integrator_t), intent(inout) :: si + integer, intent(in) :: s + type(field_can_t), intent(inout) :: fs(:) + real(dp), intent(inout) :: x(4*s) + real(dp), intent(in) :: atol, rtol + integer, intent(in) :: maxit + real(dp), intent(out) :: xlast(4*s) + logical, intent(out) :: hit_maxit + + real(dp) :: fvec(4*s), fjac(4*s, 4*s) + real(dp) :: xabs(4*s), tolref(4*s), fabs(4*s) + integer :: kit, ks, info + logical :: conv + + hit_maxit = .false. + do kit = 1, maxit + + ! Check if radius left the boundary + do ks = 1, s + if (x(4*ks-3) > 1d0) return + ! Transient guard for intermediate iterates; the converged-negative + ! case is handled by the caller via a chart switch (#370). + if (x(4*ks-3) < 0.0d0) x(4*ks-3) = 0.01d0 + end do + + call rk_gauss_residual(model, si, fs, s, x, fvec) + call rk_gauss_jacobian(model, si, fs, s, fjac) + fabs = abs(fvec) + xlast = x + call rk_solve(4*s, fjac, fvec, info) + if (info /= 0) return + ! after solution: fvec = (xold-xnew)_Newton + x = x - fvec + xabs = abs(x - xlast) + + do ks = 1, s + tolref(4*ks-3) = 1d0 + tolref(4*ks-2) = twopi + tolref(4*ks-1) = twopi + tolref(4*ks) = abs(xlast(4*ks)) + end do + + conv = .true. + do ks = 1, 4*s + if (fabs(ks) >= atol) conv = .false. + end do + if (conv) return + conv = .true. + do ks = 1, 4*s + if (xabs(ks) >= rtol*tolref(ks)) conv = .false. + end do + if (conv) return + end do + hit_maxit = .true. + end subroutine rk_gauss_newton + +end module orbit_rk_core diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index 99ddb904..50ea0e3e 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -14,6 +14,7 @@ module orbit_symplectic use orbit_symplectic_euler1, only: sympl_euler1_residual, sympl_euler1_jacobian, & sympl_euler1_newton_iter, sympl_euler1_extrapolate_field, & sympl_euler1_advance_angles +use orbit_rk_core, only: rk_gauss_residual, rk_gauss_jacobian, MODEL_GC use vector_potentail_mod, only: torflux use lapack_interfaces, only: dgesv use diag_counters, only: count_event, EVT_NEWTON1_MAXIT, EVT_NEWTON2_MAXIT, & @@ -505,38 +506,16 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast) ! Gauss-Legendre Runge-Kutta method with s internal stages (n=4*s variables) ! recursive subroutine f_rk_gauss(si, fs, s, x, fvec) - ! + ! GC Gauss residual: thin wrapper over the shared core (MODEL_GC). The + ! arithmetic lives in orbit_rk_core::gauss_canfield_residual, shared verbatim + ! with the CPP pusher. type(symplectic_integrator_t), intent(inout) :: si integer, intent(in) :: s type(field_can_t), intent(inout) :: fs(:) real(dp), intent(in) :: x(4*s) ! = (rend, thend, phend, pphend) real(dp), intent(out) :: fvec(4*s) - real(dp) :: a(s,s), b(s), c(s), Hprime(s) - integer :: k,l ! counters - - call coeff_rk_gauss(s, a, b, c) ! TODO: move this to preprocessing - - ! evaluate stages - do k = 1, s - call eval_field(fs(k), x(4*k-3), x(4*k-2), x(4*k-1), 2) - call get_derivatives2(fs(k), x(4*k)) - Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) - end do - - do k = 1, s - fvec(4*k-3) = fs(k)%pth - si%pthold - fvec(4*k-2) = x(4*k-2) - si%z(2) - fvec(4*k-1) = x(4*k-1) - si%z(3) - fvec(4*k) = x(4*k) - si%z(4) - do l = 1, s - fvec(4*k-3) = fvec(4*k-3) + si%dt*a(k,l)*(fs(l)%dH(2) - Hprime(l)*fs(l)%dpth(2)) ! pthdot - fvec(4*k-2) = fvec(4*k-2) - si%dt*a(k,l)*Hprime(l) ! thdot - fvec(4*k-1) = fvec(4*k-1) - si%dt*a(k,l)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph ! phdot - fvec(4*k) = fvec(4*k) + si%dt*a(k,l)*(fs(l)%dH(3) - Hprime(l)*fs(l)%dpth(3)) ! pphdot - end do - end do - + call rk_gauss_residual(MODEL_GC, si, fs, s, x, fvec) end subroutine f_rk_gauss @@ -544,79 +523,14 @@ end subroutine f_rk_gauss !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! recursive subroutine jac_rk_gauss(si, fs, s, jac) - ! + ! GC Gauss Jacobian: thin wrapper over the shared core (MODEL_GC). The + ! arithmetic lives in orbit_rk_core::gauss_canfield_jacobian. type(symplectic_integrator_t), intent(in) :: si integer, intent(in) :: s type(field_can_t), intent(in) :: fs(:) real(dp), intent(out) :: jac(4*s, 4*s) - real(dp) :: a(s,s), b(s), c(s), Hprime(s), dHprime(4*s) - integer :: k,l,m ! counters - - call coeff_rk_gauss(s, a, b, c) ! TODO: move this to preprocessing - - ! evaluate stages - do k = 1, s - m=4*k - Hprime(k) = fs(k)%dH(1)/fs(k)%dpth(1) - dHprime(m-3) = (fs(k)%d2H(1) - Hprime(k)*fs(k)%d2pth(1))/fs(k)%dpth(1) ! d/dr - dHprime(m-2) = (fs(k)%d2H(2) - Hprime(k)*fs(k)%d2pth(2))/fs(k)%dpth(1) ! d/dth - dHprime(m-1) = (fs(k)%d2H(3) - Hprime(k)*fs(k)%d2pth(3))/fs(k)%dpth(1) ! d/dph - dHprime(m) = (fs(k)%d2H(7) - Hprime(k)*fs(k)%d2pth(7))/fs(k)%dpth(1) ! d/dpph - end do - - jac = 0d0 - - do k = 1, s - m=4*k - jac(m-3, m-3) = fs(k)%dpth(1) - jac(m-3, m-2) = fs(k)%dpth(2) - jac(m-3, m-1) = fs(k)%dpth(3) - jac(m-3, m) = fs(k)%dpth(4) - jac(m-2, m-2) = 1d0 - jac(m-1, m-1) = 1d0 - jac(m, m) = 1d0 - end do - - do l = 1, s - do k = 1, s - m=4*k - jac(m-3, 4*l-3) = jac(m-3, 4*l-3) & ! d/dr - + si%dt*a(k,l)*(fs(l)%d2H(2) - fs(l)%d2pth(2)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-3)) - jac(m-3, 4*l-2) = jac(m-3, 4*l-2) & ! d/dth - + si%dt*a(k,l)*(fs(l)%d2H(4) - fs(l)%d2pth(4)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-2)) - jac(m-3, 4*l-1) = jac(m-3, 4*l-1) & ! d/dph - + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l-1)) - jac(m-3, 4*l) = jac(m-3, 4*l) & ! d/dpph - + si%dt*a(k,l)*(fs(l)%d2H(8) - fs(l)%d2pth(8)*Hprime(l) - fs(l)%dpth(2)*dHprime(4*l)) - - jac(m-2, 4*l-3) = jac(m-2, 4*l-3) - si%dt*a(k,l)*dHprime(4*l-3) ! d/dr - jac(m-2, 4*l-2) = jac(m-2, 4*l-2) - si%dt*a(k,l)*dHprime(4*l-2) ! d/dth - jac(m-2, 4*l-1) = jac(m-2, 4*l-1) - si%dt*a(k,l)*dHprime(4*l-1) ! d/dph - jac(m-2, 4*l) = jac(m-2, 4*l) - si%dt*a(k,l)*dHprime(4*l) ! d/dpph - - jac(m-1, 4*l-3) = jac(m-1, 4*l-3) & ! d/dr - - si%dt*a(k,l)*(-fs(l)%dhph(1)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & - + (fs(l)%dvpar(1) - dHprime(4*l-3)*fs(l)%hth - Hprime(l)*fs(l)%dhth(1))/fs(l)%hph) - jac(m-1, 4*l-2) = jac(m-1, 4*l-2) & ! d/dth - - si%dt*a(k,l)*(-fs(l)%dhph(2)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & - + (fs(l)%dvpar(2) - dHprime(4*l-2)*fs(l)%hth - Hprime(l)*fs(l)%dhth(2))/fs(l)%hph) - jac(m-1, 4*l-1) = jac(m-1, 4*l-1) & ! d/dph - - si%dt*a(k,l)*(-fs(l)%dhph(3)*(fs(l)%vpar - Hprime(l)*fs(l)%hth)/fs(l)%hph**2 & - + (fs(l)%dvpar(3) - dHprime(4*l-1)*fs(l)%hth - Hprime(l)*fs(l)%dhth(3))/fs(l)%hph) - jac(m-1, 4*l) = jac(m-1, 4*l) & ! d/dpph - - si%dt*a(k,l)*((fs(l)%dvpar(4) - dHprime(4*l)*fs(l)%hth)/fs(l)%hph) - - jac(m, 4*l-3) = jac(m, 4*l-3) & ! d/dr - + si%dt*a(k,l)*(fs(l)%d2H(3) - fs(l)%d2pth(3)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-3)) - jac(m, 4*l-2) = jac(m, 4*l-2) & ! d/dth - + si%dt*a(k,l)*(fs(l)%d2H(5) - fs(l)%d2pth(5)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-2)) - jac(m, 4*l-1) = jac(m, 4*l-1) & ! d/dph - + si%dt*a(k,l)*(fs(l)%d2H(6) - fs(l)%d2pth(6)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l-1)) - jac(m, 4*l) = jac(m, 4*l) & ! d/dpph - + si%dt*a(k,l)*(fs(l)%d2H(9) - fs(l)%d2pth(9)*Hprime(l) - fs(l)%dpth(3)*dHprime(4*l)) - end do - end do + call rk_gauss_jacobian(MODEL_GC, si, fs, s, jac) end subroutine jac_rk_gauss diff --git a/src/params.f90 b/src/params.f90 index 00b9afa0..c8d1da92 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -45,6 +45,10 @@ module params integer :: integmode = EXPL_IMPL_EULER + ! Orbit model selector: 0 guiding-center (default, symplectic GC path), + ! 1 Pauli/CPP full orbit, 2 Boris full orbit. See src/orbit_full.f90. + integer :: orbit_model = 0 + integer :: kpart = 0 ! progress counter for particles real(dp) :: relerr = 1d-13 @@ -109,7 +113,7 @@ module params trace_time, num_surf, sbeg, phibeg, thetabeg, contr_pp, & facE_al, npoiper2, n_e, n_d, netcdffile, ns_s, ns_tp, multharm, & isw_field_type, generate_start_only, startmode, grid_density, & - special_ants_file, integmode, relerr, tcut, nturns, debug, & + special_ants_file, integmode, orbit_model, relerr, tcut, nturns, debug, & class_plot, cut_in_per, fast_class, vmec_B_scale, & vmec_RZ_scale, swcoll, deterministic, old_axis_healing, & old_axis_healing_boundary, am1, am2, Z1, Z2, & diff --git a/src/simple_main.f90 b/src/simple_main.f90 index f4502f4b..3bf67b11 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -875,6 +875,9 @@ end subroutine trace_orbit subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local) use alpha_lifetime_sub, only: orbit_timestep_axis use orbit_symplectic, only: orbit_timestep_sympl + use orbit_cpp, only: orbit_timestep_cpp, cpp_stages_from_mode + use orbit_full, only: ORBIT_PAULI + use params, only: orbit_model type(tracer_t), intent(inout) :: anorb real(dp), intent(inout) :: z(5) @@ -889,7 +892,16 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local) call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr_orbit) else if (swcoll) call update_momentum(anorb, z) - call orbit_timestep_sympl(anorb%si, anorb%f, ierr_orbit) + ! Dispatch by integer orbit_model: GC (default) uses the + ! symplectic GC pusher; PAULI (CPP) integrates the same 4D + ! canonical state with mu held fixed on the slow manifold. + select case (orbit_model) + case (ORBIT_PAULI) + call orbit_timestep_cpp(anorb%si, anorb%f, & + cpp_stages_from_mode(integmode), ierr_orbit) + case default + call orbit_timestep_sympl(anorb%si, anorb%f, ierr_orbit) + end select call to_standard_z_coordinates(anorb, z) end if if (swcoll) call collide(z, dtaumin) ! Collisions diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 9796d0a3..4c39a52c 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -570,6 +570,34 @@ add_executable(test_orbit_symplectic_base.x test_orbit_symplectic_base.f90) target_link_libraries(test_orbit_symplectic_base.x simple) add_test(NAME test_orbit_symplectic_base COMMAND test_orbit_symplectic_base.x) +add_executable(test_full_orbit.x test_full_orbit.f90) +target_link_libraries(test_full_orbit.x simple) +add_test(NAME test_full_orbit COMMAND test_full_orbit.x) +set_tests_properties(test_full_orbit PROPERTIES LABELS "unit") + +# Implicit-midpoint full-orbit pusher on the analytic mocks (no field file). +add_executable(test_fo_symplectic.x test_fo_symplectic.f90) +target_link_libraries(test_fo_symplectic.x simple) +add_test(NAME test_fo_symplectic COMMAND test_fo_symplectic.x) +set_tests_properties(test_fo_symplectic PROPERTIES LABELS "unit" TIMEOUT 120) + +# orbit_model config parsing + dispatch keys (Wave-1 followup #1). +add_executable(test_orbit_model_dispatch.x test_orbit_model_dispatch.f90) +target_link_libraries(test_orbit_model_dispatch.x simple) +add_test(NAME test_orbit_model_dispatch COMMAND test_orbit_model_dispatch.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_orbit_model_dispatch PROPERTIES LABELS "unit") + +# CPP (flux-canonical) cross-checks against the BOOZER chart of the QA wout. +# Run from the binary dir where the wout.nc symlink resolves. +foreach(cpp_test test_cpp_equals_gc_largestep test_cpp_invariants) + add_executable(${cpp_test}.x ${cpp_test}.f90) + target_link_libraries(${cpp_test}.x simple) + add_test(NAME ${cpp_test} COMMAND ${cpp_test}.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties(${cpp_test} PROPERTIES LABELS "integration" TIMEOUT 120) +endforeach() + add_executable(test_field_base.x test_field_base.f90) target_link_libraries(test_field_base.x simple) add_test(NAME test_field_base COMMAND test_field_base.x) diff --git a/test/tests/test_cpp_equals_gc_largestep.f90 b/test/tests/test_cpp_equals_gc_largestep.f90 new file mode 100644 index 00000000..163b0b71 --- /dev/null +++ b/test/tests/test_cpp_equals_gc_largestep.f90 @@ -0,0 +1,110 @@ +program test_cpp_equals_gc_largestep + ! CPP (flux-canonical, mu fixed) must reproduce the GC symplectic trajectory + ! at GC-sized steps. The CPP Gauss residual is the GC degenerate-Lagrangian + ! Euler-Lagrange system specialized to fixed mu, so on the same canonical + ! chart, same stage, same dt the two integrators advance the 4D state + ! z=(r,theta,phi,pphi) identically up to the solver tolerance. + ! + ! Cheap real field: BOOZER chart on the QA wout. Cross-check (strongest + ! available): max over the orbit of |z_CPP - z_GC| < 1e-10. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: twopi + use simple_main, only: init_field + use simple, only: tracer_t, init_sympl, init_params + use params, only: isw_field_type, field_input, coord_input, integmode + use new_vmec_stuff_mod, only: rmajor + use magfie_sub, only: BOOZER + use field_can_mod, only: field_can_t, evaluate + use orbit_symplectic_base, only: symplectic_integrator_t, GAUSS1, GAUSS2 + use orbit_symplectic, only: orbit_timestep_sympl + use orbit_cpp, only: orbit_timestep_cpp, orbit_cpp_init, cpp_stages_from_mode + + implicit none + + integer, parameter :: ans_s = 5, ans_tp = 5, amultharm = 5 + integer :: nfail + type(tracer_t) :: norb + + nfail = 0 + + isw_field_type = BOOZER + field_input = 'wout.nc' + coord_input = 'wout.nc' + call init_field(norb, 'wout.nc', ans_s, ans_tp, amultharm, GAUSS1) + call init_params(norb, 2, 4, 3.5e6_dp, 256, 1, 1.0e-12_dp) + if (.not. associated(evaluate)) then + print *, 'evaluate pointer not associated for BOOZER field' + error stop 1 + end if + + call run_compare(norb, GAUSS1, 'GAUSS1', nfail) + call run_compare(norb, GAUSS2, 'GAUSS2', nfail) + + if (nfail == 0) then + print *, 'ALL CPP==GC LARGE-STEP TESTS PASSED' + else + print *, 'CPP==GC TESTS FAILED: ', nfail + error stop 1 + end if + +contains + + subroutine run_compare(norb, mode, tag, nfail) + type(tracer_t), intent(inout) :: norb + integer, intent(in) :: mode + character(*), intent(in) :: tag + integer, intent(inout) :: nfail + + type(symplectic_integrator_t) :: si_gc, si_cpp + type(field_can_t) :: f_gc, f_cpp + real(dp) :: z0(5), dtau, dtaumin, relerr, rbig, maxdiff, d + integer :: it, ierr_gc, ierr_cpp, s, nstep + + rbig = rmajor*1.0e2_dp + ! GC-sized step: 256 points per torus (production-class resolution). + dtau = twopi*rbig/256.0_dp + dtaumin = dtau + relerr = 1.0e-13_dp + nstep = 300 + s = cpp_stages_from_mode(mode) + + ! Trapped-particle-class IC. + z0 = [0.4_dp, 0.7_dp, 0.1_dp, 1.0_dp, 0.2_dp] + + integmode = mode + norb%relerr = relerr + + call init_sympl(si_gc, f_gc, z0, dtau, dtaumin, relerr, mode) + call init_sympl(si_cpp, f_cpp, z0, dtau, dtaumin, relerr, mode) + call orbit_cpp_init(si_cpp, f_cpp) + + maxdiff = 0.0_dp + do it = 1, nstep + call orbit_timestep_sympl(si_gc, f_gc, ierr_gc) + call orbit_timestep_cpp(si_cpp, f_cpp, s, ierr_cpp) + if (ierr_gc /= 0 .or. ierr_cpp /= 0) then + print '(A,A,A,I0,A,I0)', ' ', tag, ': early exit ierr_gc=', ierr_gc, & + ' ierr_cpp=', ierr_cpp + exit + end if + d = maxval(abs(si_cpp%z - si_gc%z)) + maxdiff = max(maxdiff, d) + end do + + print '(A,A,A,ES12.4)', ' ', tag, ': max |z_CPP - z_GC| = ', maxdiff + call check(tag//': CPP matches GC to Newton tol', maxdiff < 1.0e-10_dp, nfail) + end subroutine run_compare + + subroutine check(name, ok, nfail) + character(*), intent(in) :: name + logical, intent(in) :: ok + integer, intent(inout) :: nfail + if (ok) then + print '(A,A)', 'PASS ', name + else + print '(A,A)', 'FAIL ', name + nfail = nfail + 1 + end if + end subroutine check + +end program test_cpp_equals_gc_largestep diff --git a/test/tests/test_cpp_invariants.f90 b/test/tests/test_cpp_invariants.f90 new file mode 100644 index 00000000..0bab3d5e --- /dev/null +++ b/test/tests/test_cpp_invariants.f90 @@ -0,0 +1,156 @@ +program test_cpp_invariants + ! Invariant conservation for the CPP pusher on the BOOZER chart of the QA wout + ! (cheap real field). CPP is the GC degenerate-Lagrangian scheme with mu held + ! fixed, so its conserved-quantity behavior must match the validated GC + ! integrator. Asserts: + ! - mu is a fixed parameter -> identically conserved (byte ==). + ! - energy H = vpar^2/2 + mu*Bmod oscillation, energy secular drift, and the + ! canonical p_phi excursion each match GC to a tight relative margin (CPP + ! is no worse than GC at structure preservation), and are bounded with no + ! secular energy growth (modified-Hamiltonian property). + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: twopi + use simple_main, only: init_field + use simple, only: tracer_t, init_sympl, init_params + use params, only: isw_field_type, field_input, coord_input, integmode + use new_vmec_stuff_mod, only: rmajor + use magfie_sub, only: BOOZER + use field_can_mod, only: field_can_t, evaluate, get_val + use orbit_symplectic_base, only: symplectic_integrator_t, GAUSS1 + use orbit_symplectic, only: orbit_timestep_sympl + use orbit_cpp, only: orbit_timestep_cpp, orbit_cpp_init, cpp_stages_from_mode + + implicit none + + integer, parameter :: ans_s = 5, ans_tp = 5, amultharm = 5 + integer :: nfail + type(tracer_t) :: norb + + nfail = 0 + + isw_field_type = BOOZER + field_input = 'wout.nc' + coord_input = 'wout.nc' + call init_field(norb, 'wout.nc', ans_s, ans_tp, amultharm, GAUSS1) + call init_params(norb, 2, 4, 3.5e6_dp, 256, 1, 1.0e-12_dp) + if (.not. associated(evaluate)) then + print *, 'evaluate pointer not associated for BOOZER field' + error stop 1 + end if + + call run_invariants(norb, nfail) + + if (nfail == 0) then + print *, 'ALL CPP INVARIANT TESTS PASSED' + else + print *, 'CPP INVARIANT TESTS FAILED: ', nfail + error stop 1 + end if + +contains + + subroutine run_invariants(norb, nfail) + type(tracer_t), intent(inout) :: norb + integer, intent(inout) :: nfail + + real(dp) :: z0(5), dtau, relerr, rbig + real(dp) :: mu_cpp, mu_cpp_end + real(dp) :: osc_cpp, sec_cpp, pdev_cpp + real(dp) :: osc_gc, sec_gc, pdev_gc + integer :: s, nstep + + rbig = rmajor*1.0e2_dp + dtau = twopi*rbig/256.0_dp + relerr = 1.0e-13_dp + nstep = 3000 + s = cpp_stages_from_mode(GAUSS1) + + z0 = [0.5_dp, 0.0_dp, 0.0_dp, 1.0_dp, 0.1_dp] + integmode = GAUSS1 + norb%relerr = relerr + + call run_one(z0, dtau, relerr, s, .true., osc_cpp, sec_cpp, pdev_cpp, & + mu_cpp, mu_cpp_end) + call run_one(z0, dtau, relerr, s, .false., osc_gc, sec_gc, pdev_gc, & + mu_cpp, mu_cpp_end) + + print '(A,ES12.4)', ' mu (fixed) = ', mu_cpp + print '(A,ES12.4,A,ES12.4)', ' energy osc: CPP=', osc_cpp, ' GC=', osc_gc + print '(A,ES12.4,A,ES12.4)', ' energy secular: CPP=', sec_cpp, ' GC=', sec_gc + print '(A,ES12.4,A,ES12.4)', ' pphi excursion: CPP=', pdev_cpp, ' GC=', pdev_gc + + ! mu is a fixed parameter: it must not move at all. + call check('mu identically conserved', mu_cpp == mu_cpp_end, nfail) + ! no secular energy growth (structure preservation). + call check('energy no secular drift', sec_cpp < 0.5_dp*osc_cpp + 1.0e-12_dp, & + nfail) + ! CPP conserves no worse than the validated GC integrator (relative margin). + call check('energy osc <= GC', osc_cpp <= osc_gc*(1.0_dp + 1.0e-6_dp) + 1.0e-12_dp, & + nfail) + call check('energy secular <= GC', sec_cpp <= sec_gc*(1.0_dp + 1.0e-6_dp) + 1.0e-12_dp, & + nfail) + call check('pphi excursion ~ GC', & + abs(pdev_cpp - pdev_gc) <= 1.0e-6_dp*max(pdev_gc, 1.0e-12_dp) + 1.0e-10_dp, & + nfail) + end subroutine run_invariants + + subroutine run_one(z0, dtau, relerr, s, use_cpp, osc, sec, pdev, mu0, muend) + real(dp), intent(in) :: z0(5), dtau, relerr + integer, intent(in) :: s + logical, intent(in) :: use_cpp + real(dp), intent(out) :: osc, sec, pdev, mu0, muend + + type(symplectic_integrator_t) :: si + type(field_can_t) :: f + real(dp) :: H0, H, Hmin, Hmax, pphi0, H_first, H_last + integer :: it, ierr, nhalf, n1, n2 + + call init_sympl(si, f, z0, dtau, dtau, relerr, GAUSS1) + if (use_cpp) call orbit_cpp_init(si, f) + + mu0 = f%mu + pphi0 = si%z(4) + call get_val(f, si%z(4)) + H0 = f%H; Hmin = H0; Hmax = H0; pdev = 0.0_dp + nhalf = 3000/2 + H_first = 0.0_dp; H_last = 0.0_dp; n1 = 0; n2 = 0 + + do it = 1, 3000 + if (use_cpp) then + call orbit_timestep_cpp(si, f, s, ierr) + else + call orbit_timestep_sympl(si, f, ierr) + end if + if (ierr /= 0) then + print '(A,L1,A,I0)', ' early exit use_cpp=', use_cpp, ' at step ', it + exit + end if + call get_val(f, si%z(4)) + H = f%H + Hmin = min(Hmin, H); Hmax = max(Hmax, H) + pdev = max(pdev, abs(si%z(4) - pphi0)) + if (it <= nhalf) then + H_first = H_first + H; n1 = n1 + 1 + else + H_last = H_last + H; n2 = n2 + 1 + end if + end do + + muend = f%mu + osc = (Hmax - Hmin)/abs(H0) + sec = abs(H_last/max(n2,1) - H_first/max(n1,1))/abs(H0) + end subroutine run_one + + subroutine check(name, ok, nfail) + character(*), intent(in) :: name + logical, intent(in) :: ok + integer, intent(inout) :: nfail + if (ok) then + print '(A,A)', 'PASS ', name + else + print '(A,A)', 'FAIL ', name + nfail = nfail + 1 + end if + end subroutine check + +end program test_cpp_invariants diff --git a/test/tests/test_fo_symplectic.f90 b/test/tests/test_fo_symplectic.f90 new file mode 100644 index 00000000..ead5a290 --- /dev/null +++ b/test/tests/test_fo_symplectic.f90 @@ -0,0 +1,205 @@ +program test_fo_symplectic + ! Behavioral tests for the implicit-midpoint curvilinear full-orbit pusher + ! (ORBIT_FOSYMPL, foimpl_step) on the analytic mock providers. Oracles: + ! 1. uniform-B: |v| and energy conserved to solver tolerance; closed circle. + ! 2. cylindrical 1/R curvature drift vs analytic v_d (geodesic + Lorentz). + ! 3. long-run energy: no secular drift (structure-preserving midpoint), the + ! property an explicit RK does not have. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: c, twopi, p_mass, e_charge + use orbit_full, only: FullOrbitState, init_full_orbit_state, & + timestep_full_orbit, compute_energy, ORBIT_FOSYMPL, COORD_CART, COORD_CYL, & + FO_OK + use orbit_full_mock_cart, only: cartesian_provider_t, FIELD_UNIFORM + use orbit_full_mock_cyl, only: cylindrical_provider_t + implicit none + + integer :: nfail + nfail = 0 + + call test_uniform_gyration(nfail) + call test_cyl_curvature_drift(nfail) + call test_energy_no_secular_drift(nfail) + + if (nfail == 0) then + print *, 'ALL FO-SYMPLECTIC TESTS PASSED' + else + print *, 'FO-SYMPLECTIC TESTS FAILED: ', nfail + error stop 1 + end if + +contains + + subroutine check(name, ok, nfail) + character(*), intent(in) :: name + logical, intent(in) :: ok + integer, intent(inout) :: nfail + if (ok) then + print '(A,A)', 'PASS ', name + else + print '(A,A)', 'FAIL ', name + nfail = nfail + 1 + end if + end subroutine check + + subroutine test_uniform_gyration(nfail) + integer, intent(inout) :: nfail + type(cartesian_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, vperp, vpar, speed0 + real(dp) :: Omega, period, dt, rL + real(dp) :: x0(3), v0(3), xstart(3) + real(dp) :: e0, e1, errpos, errv + integer :: i, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + prov%field_kind = FIELD_UNIFORM + prov%B0 = [0.0_dp, 0.0_dp, B0] + + vperp = 1.0d7 + vpar = 3.0d6 + speed0 = sqrt(vperp**2 + vpar**2) + Omega = charge * B0 / (mass * c) + period = twopi / Omega + rL = mass * c * vperp / (charge * B0) + nstep = 400 + dt = period / nstep + + x0 = [0.0_dp, 0.0_dp, 0.0_dp] + v0 = [vperp, 0.0_dp, vpar] + call init_full_orbit_state(st, x0, v0, ORBIT_FOSYMPL, COORD_CART, & + mass, charge, dt, prov) + xstart = st%z(1:3) + e0 = compute_energy(st) + + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('uniform: timestep ierr', .false., nfail) + return + end if + end do + + e1 = compute_energy(st) + errv = abs(sqrt(dot_product(st%z(4:6), st%z(4:6))) - speed0) / speed0 + errpos = sqrt((st%z(1)-xstart(1))**2 + (st%z(2)-xstart(2))**2) + + print '(A,ES12.4,A,ES12.4)', ' uniform: |v| relerr=', errv, & + ' return-pos err=', errpos + print '(A,ES12.4)', ' uniform: dE/E=', abs(e1-e0)/e0 + + call check('uniform: |v| constant', errv < 1d-9, nfail) + call check('uniform: energy constant', abs(e1-e0)/e0 < 1d-9, nfail) + call check('uniform: return to start', errpos < 1d-2*rL, nfail) + end subroutine test_uniform_gyration + + subroutine test_cyl_curvature_drift(nfail) + integer, intent(inout) :: nfail + type(cylindrical_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, R0, R, Bloc + real(dp) :: Omega, period, dt, vd_exact, vd_meas + real(dp) :: u0(3), w0(3), z0, t_total, vpar_in, vperp_in + integer :: i, nstep_per, nper_run, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + R0 = 200.0_dp + R = 200.0_dp + prov%B0 = B0 + prov%R0 = R0 + + vpar_in = 1.0d7 + vperp_in = 1.0d5 + Bloc = B0 * R0 / R + Omega = charge * Bloc / (mass * c) + period = twopi / Omega + nstep_per = 300 + nper_run = 40 + nstep = nstep_per * nper_run + dt = period / nstep_per + + vd_exact = (mass * c) / (charge * Bloc * R) * & + (vpar_in**2 + 0.5_dp * vperp_in**2) + + u0 = [R, 0.0_dp, 0.0_dp] + w0 = [vperp_in, vpar_in / R, 0.0_dp] + call init_full_orbit_state(st, u0, w0, ORBIT_FOSYMPL, COORD_CYL, & + mass, charge, dt, prov) + z0 = st%z(3) + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('cyl curvature: timestep ierr', .false., nfail) + return + end if + end do + t_total = nstep * dt + vd_meas = (st%z(3) - z0) / t_total + + print '(A,ES12.4,A,ES12.4)', ' cyl curvature: vd_exact=', vd_exact, & + ' vd_meas=', vd_meas + call check('cyl curvature: vertical drift', & + abs(vd_meas - vd_exact) < 0.10_dp*abs(vd_exact), nfail) + end subroutine test_cyl_curvature_drift + + subroutine test_energy_no_secular_drift(nfail) + ! Long run in uniform B: the midpoint energy stays bounded with no secular + ! growth. Compare mean energy over first vs last quarter of the run. + integer, intent(inout) :: nfail + type(cartesian_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, vperp, vpar + real(dp) :: Omega, period, dt, e0, e, emin, emax + real(dp) :: efirst, elast, x0(3), v0(3) + integer :: i, nstep, nq, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + prov%field_kind = FIELD_UNIFORM + prov%B0 = [0.0_dp, 0.0_dp, B0] + + vperp = 1.0d7 + vpar = 2.0d6 + Omega = charge * B0 / (mass * c) + period = twopi / Omega + nstep = 200 * 200 + dt = period / 200 + + x0 = [0.0_dp, 0.0_dp, 0.0_dp] + v0 = [vperp, 0.0_dp, vpar] + call init_full_orbit_state(st, x0, v0, ORBIT_FOSYMPL, COORD_CART, & + mass, charge, dt, prov) + e0 = compute_energy(st) + emin = e0; emax = e0 + nq = nstep/4 + efirst = 0.0_dp; elast = 0.0_dp + + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('no-drift: timestep ierr', .false., nfail) + return + end if + e = compute_energy(st) + emin = min(emin, e); emax = max(emax, e) + if (i <= nq) efirst = efirst + e + if (i > nstep - nq) elast = elast + e + end do + efirst = efirst / nq + elast = elast / nq + + print '(A,ES12.4)', ' no-drift: energy bound (emax-emin)/e0=', (emax-emin)/e0 + print '(A,ES12.4)', ' no-drift: secular |_last-_first|/e0=', & + abs(elast - efirst)/e0 + + call check('fosympl: energy bounded long run', (emax-emin)/e0 < 1d-8, nfail) + call check('fosympl: no secular energy drift', & + abs(elast - efirst)/e0 < 1d-10, nfail) + end subroutine test_energy_no_secular_drift + +end program test_fo_symplectic diff --git a/test/tests/test_full_orbit.f90 b/test/tests/test_full_orbit.f90 new file mode 100644 index 00000000..61892e60 --- /dev/null +++ b/test/tests/test_full_orbit.f90 @@ -0,0 +1,392 @@ +program test_full_orbit + ! Behavioral tests for the full-orbit Boris pusher using only the analytic + ! Cartesian and cylindrical mock providers (no libneo field). Oracles: + ! 1. uniform-B exact gyration: |v|, vpar, energy, mu conserved; closed + ! circle of Larmor radius r_L = m c vperp/(q B), period T = 2pi m c/(qB). + ! 2. Cartesian linear grad-B drift vs analytic v_gradB. + ! 3. cylindrical 1/R curvature + grad-B drift vs analytic v_d, separated + ! into pure curvature (vperp=0) and pure grad-B (vpar=0). + ! 4. mu adiabatic invariance over many gyroperiods. + ! 5. cylindrical Christoffel mock vs closed form. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use util, only: c, twopi, pi, p_mass, e_charge + use orbit_full, only: FullOrbitState, init_full_orbit_state, & + timestep_full_orbit, convert_full_to_gc, compute_energy, & + ORBIT_BORIS, COORD_CART, COORD_CYL, FO_OK + use orbit_full_mock_cart, only: cartesian_provider_t, FIELD_UNIFORM, FIELD_LINGRAD + use orbit_full_mock_cyl, only: cylindrical_provider_t + implicit none + + integer :: nfail + nfail = 0 + + call test_uniform_gyration(nfail) + call test_cart_gradb_drift(nfail) + call test_cyl_curvature_drift(nfail) + call test_cyl_gradb_drift(nfail) + call test_mu_invariance(nfail) + call test_cyl_christoffel(nfail) + + if (nfail == 0) then + print *, 'ALL FULL-ORBIT TESTS PASSED' + else + print *, 'FULL-ORBIT TESTS FAILED: ', nfail + error stop 1 + end if + +contains + + subroutine check(name, ok, nfail) + character(*), intent(in) :: name + logical, intent(in) :: ok + integer, intent(inout) :: nfail + if (ok) then + print '(A,A)', 'PASS ', name + else + print '(A,A)', 'FAIL ', name + nfail = nfail + 1 + end if + end subroutine check + + subroutine test_uniform_gyration(nfail) + integer, intent(inout) :: nfail + type(cartesian_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, vperp, vpar, speed0 + real(dp) :: Omega, period, dt, rL + real(dp) :: x0(3), v0(3), xstart(3) + real(dp) :: e0, e1, mu0, vz0, errpos, errv + integer :: i, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + prov%field_kind = FIELD_UNIFORM + prov%B0 = [0.0_dp, 0.0_dp, B0] + + vperp = 1.0d7 + vpar = 3.0d6 + speed0 = sqrt(vperp**2 + vpar**2) + Omega = charge * B0 / (mass * c) + period = twopi / Omega + rL = mass * c * vperp / (charge * B0) + nstep = 400 + dt = period / nstep + + x0 = [0.0_dp, 0.0_dp, 0.0_dp] + v0 = [vperp, 0.0_dp, vpar] + call init_full_orbit_state(st, x0, v0, ORBIT_BORIS, COORD_CART, & + mass, charge, dt, prov) + xstart = st%z(1:3) + e0 = compute_energy(st) + mu0 = st%mu + vz0 = st%z(6) + + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('uniform: timestep ierr', .false., nfail) + return + end if + end do + + e1 = compute_energy(st) + errv = abs(sqrt(dot_product(st%z(4:6), st%z(4:6))) - speed0) / speed0 + ! transverse return: z advances by vpar*period, so only (x,y) close. + errpos = sqrt((st%z(1)-xstart(1))**2 + (st%z(2)-xstart(2))**2) + + print '(A,ES12.4,A,ES12.4)', ' uniform: r_L=', rL, ' period=', period + print '(A,ES12.4,A,ES12.4)', ' uniform: |v| relerr=', errv, & + ' return-pos err=', errpos + print '(A,ES12.4)', ' uniform: dE/E=', abs(e1-e0)/e0 + print '(A,ES12.4)', ' uniform: dmu/mu=', abs(st%mu-mu0)/mu0 + + call check('uniform: |v| constant', errv < 1d-10, nfail) + call check('uniform: energy constant', abs(e1-e0)/e0 < 1d-9, nfail) + call check('uniform: vpar=vz constant', abs(st%z(6)-vz0) < 1d-6*speed0, nfail) + ! closed circle: error scales with dt; Boris is 2nd order -> bound on rL. + call check('uniform: return to start', errpos < 1d-3*rL, nfail) + call check('uniform: mu constant', abs(st%mu-mu0)/mu0 < 1d-9, nfail) + end subroutine test_uniform_gyration + + subroutine test_cart_gradb_drift(nfail) + integer, intent(inout) :: nfail + type(cartesian_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, g, vperp, vpar + real(dp) :: Omega, period, dt, vd_exact, vd_meas + real(dp) :: x0(3), v0(3), ygc0, ygc1, t_total + integer :: i, nper_run, nstep_per, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + g = 1.0d2 ! G/cm gradient of B_z along x + prov%field_kind = FIELD_LINGRAD + prov%B0 = [0.0_dp, 0.0_dp, B0] + prov%gradB = 0.0_dp + prov%gradB(3,1) = g ! B_z = B0 + g*x + + vperp = 1.0d7 + vpar = 0.0_dp + Omega = charge * B0 / (mass * c) + period = twopi / Omega + nstep_per = 200 + nper_run = 2000 + nstep = nstep_per * nper_run + dt = period / nstep_per + + ! analytic grad-B drift: v = (m c vperp^2)/(2 q B0) * (g/B0), along +e_y for + ! q>0, B along +z, grad|B| along +x. + vd_exact = mass * c * vperp**2 / (2.0_dp * charge * B0) * (g / B0) + + x0 = [0.0_dp, 0.0_dp, 0.0_dp] + v0 = [vperp, 0.0_dp, vpar] + call init_full_orbit_state(st, x0, v0, ORBIT_BORIS, COORD_CART, & + mass, charge, dt, prov) + ygc0 = guiding_center_y_cart(st) + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('cart gradB: timestep ierr', .false., nfail) + return + end if + end do + ygc1 = guiding_center_y_cart(st) + t_total = nstep * dt + ! guiding-center y-drift removes the gyration, leaving the secular drift. + vd_meas = (ygc1 - ygc0) / t_total + + print '(A,ES12.4,A,ES12.4)', ' cart gradB: vd_exact=', vd_exact, & + ' vd_meas=', vd_meas + call check('cart gradB: drift sign/magnitude', & + abs(vd_meas - vd_exact) < 0.05_dp*abs(vd_exact), nfail) + end subroutine test_cart_gradb_drift + + ! Guiding-center y from the instantaneous Cartesian state: + ! x_gc = x - (1/Omega) (b x v), Omega = qB/(mc), b = B/|B|. + function guiding_center_y_cart(st) result(ygc) + type(FullOrbitState), intent(in) :: st + real(dp) :: ygc + real(dp) :: Bvec(3), Bmod, hcov(3), Omega, rho(3) + integer :: ierr + call st%prov%eval_field(st%z(1:3), Bvec, Bmod, hcov, ierr) + Omega = st%qm * Bmod / c + rho = cross_local(hcov, st%z(4:6)) / Omega + ygc = st%z(2) - rho(2) + end function guiding_center_y_cart + + pure function cross_local(a, b) result(cc) + real(dp), intent(in) :: a(3), b(3) + real(dp) :: cc(3) + cc(1) = a(2)*b(3) - a(3)*b(2) + cc(2) = a(3)*b(1) - a(1)*b(3) + cc(3) = a(1)*b(2) - a(2)*b(1) + end function cross_local + + subroutine test_cyl_curvature_drift(nfail) + ! Pure curvature: vperp -> 0 (small), vpar finite. v_d = (mc/qBR) vpar^2. + integer, intent(inout) :: nfail + call run_cyl_drift(nfail, 'cyl curvature', vpar_in=1.0d7, vperp_in=1.0d5) + end subroutine test_cyl_curvature_drift + + subroutine test_cyl_gradb_drift(nfail) + ! Pure grad-B: vpar -> 0 (small), vperp finite. v_d = (mc/qBR)(vperp^2/2). + integer, intent(inout) :: nfail + call run_cyl_drift(nfail, 'cyl gradB', vpar_in=1.0d5, vperp_in=1.0d7) + end subroutine test_cyl_gradb_drift + + subroutine run_cyl_drift(nfail, tag, vpar_in, vperp_in) + integer, intent(inout) :: nfail + character(*), intent(in) :: tag + real(dp), intent(in) :: vpar_in, vperp_in + type(cylindrical_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, R0, R, Bloc + real(dp) :: Omega, period, dt, vd_exact, vd_meas + real(dp) :: u0(3), w0(3), z0, t_total + integer :: i, nstep_per, nper_run, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + R0 = 200.0_dp + R = 200.0_dp + prov%B0 = B0 + prov%R0 = R0 + + Bloc = B0 * R0 / R ! = B0 here since R=R0 + Omega = charge * Bloc / (mass * c) + period = twopi / Omega + nstep_per = 300 + nper_run = 40 + nstep = nstep_per * nper_run + dt = period / nstep_per + + ! v_d = (m c)/(q B R) * (vpar^2 + vperp^2/2), along +Z for q>0, B toroidal. + vd_exact = (mass * c) / (charge * Bloc * R) * & + (vpar_in**2 + 0.5_dp * vperp_in**2) + + ! contravariant velocity in (R,phi,Z): orthonormal vphi_phys=vpar (toroidal + ! is field direction); vperp put into v^R (radial). v^phi = vpar/R. + u0 = [R, 0.0_dp, 0.0_dp] + w0 = [vperp_in, vpar_in / R, 0.0_dp] + call init_full_orbit_state(st, u0, w0, ORBIT_BORIS, COORD_CYL, & + mass, charge, dt, prov) + z0 = st%z(3) + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check(tag//': timestep ierr', .false., nfail) + return + end if + end do + t_total = nstep * dt + vd_meas = (st%z(3) - z0) / t_total + + print '(A,A,A,ES12.4,A,ES12.4)', ' ', tag, ': vd_exact=', vd_exact, & + ' vd_meas=', vd_meas + call check(tag//': vertical drift', & + abs(vd_meas - vd_exact) < 0.10_dp*abs(vd_exact), nfail) + end subroutine run_cyl_drift + + subroutine test_mu_invariance(nfail) + integer, intent(inout) :: nfail + type(cartesian_provider_t), target :: prov + type(FullOrbitState) :: st + real(dp) :: mass, charge, B0, g, vperp, vpar + real(dp) :: Omega, period, dt, mu0, mumax_dev, mu_now + real(dp) :: Bvec(3), Bmod, hcov(3), vperp2, vpar_now + real(dp) :: x0(3), v0(3) + integer :: i, nstep, ierr + + mass = 4.0_dp * p_mass + charge = 2.0_dp * e_charge + B0 = 1.0d4 + g = 1.0d0 + prov%field_kind = FIELD_LINGRAD + prov%B0 = [0.0_dp, 0.0_dp, B0] + prov%gradB = 0.0_dp + prov%gradB(3,1) = g + + vperp = 1.0d7 + vpar = 2.0d6 + Omega = charge * B0 / (mass * c) + period = twopi / Omega + nstep = 200 * 60 + dt = period / 200 + + x0 = [0.0_dp, 0.0_dp, 0.0_dp] + v0 = [vperp, 0.0_dp, vpar] + call init_full_orbit_state(st, x0, v0, ORBIT_BORIS, COORD_CART, & + mass, charge, dt, prov) + mu0 = st%mu + mumax_dev = 0.0_dp + do i = 1, nstep + call timestep_full_orbit(st, ierr) + if (ierr /= FO_OK) then + call check('mu inv: timestep ierr', .false., nfail) + return + end if + call prov%eval_field(st%z(1:3), Bvec, Bmod, hcov, ierr) + vpar_now = dot_product(st%z(4:6), hcov) + vperp2 = dot_product(st%z(4:6), st%z(4:6)) - vpar_now**2 + mu_now = mass * vperp2 / (2.0_dp * Bmod) + mumax_dev = max(mumax_dev, abs(mu_now - mu0) / mu0) + end do + + print '(A,ES12.4)', ' mu inv: max rel deviation=', mumax_dev + call check('mu adiabatic invariance', mumax_dev < 1d-3, nfail) + end subroutine test_mu_invariance + + subroutine test_cyl_christoffel(nfail) + integer, intent(inout) :: nfail + type(cylindrical_provider_t) :: prov + real(dp) :: Gamma(3,3,3), x(3), R + real(dp) :: gfd(3,3,3), err + logical :: ok + + prov%B0 = 1.0_dp + prov%R0 = 1.0_dp + R = 1.7_dp + x = [R, 0.3_dp, -0.5_dp] + call prov%christoffel(x, Gamma) + + ! closed-form check + ok = abs(Gamma(1,2,2) - (-R)) < 1d-12 .and. & + abs(Gamma(2,1,2) - 1.0_dp/R) < 1d-12 .and. & + abs(Gamma(2,2,1) - 1.0_dp/R) < 1d-12 + call check('cyl christoffel: closed form entries', ok, nfail) + + ! symmetry Gamma^i_{mn} = Gamma^i_{nm} + ok = christoffel_symmetric(Gamma) + call check('cyl christoffel: symmetry', ok, nfail) + + ! all other entries zero + call check('cyl christoffel: only known nonzeros', & + only_known_nonzero(Gamma, R), nfail) + + ! finite-difference metric -> Gamma, compare to closed form + call christoffel_fd(prov, x, gfd) + err = maxval(abs(gfd - Gamma)) + print '(A,ES12.4)', ' cyl christoffel: max FD-vs-analytic err=', err + call check('cyl christoffel: FD agrees', err < 1d-5, nfail) + end subroutine test_cyl_christoffel + + logical function christoffel_symmetric(Gamma) result(ok) + real(dp), intent(in) :: Gamma(3,3,3) + integer :: i, m, n + ok = .true. + do i = 1, 3 + do m = 1, 3 + do n = 1, 3 + if (abs(Gamma(i,m,n) - Gamma(i,n,m)) > 1d-14) ok = .false. + end do + end do + end do + end function christoffel_symmetric + + logical function only_known_nonzero(Gamma, R) result(ok) + real(dp), intent(in) :: Gamma(3,3,3), R + real(dp) :: g(3,3,3) + g = Gamma + g(1,2,2) = 0.0_dp + g(2,1,2) = 0.0_dp + g(2,2,1) = 0.0_dp + ok = maxval(abs(g)) < 1d-14 + end function only_known_nonzero + + subroutine christoffel_fd(prov, x, gfd) + type(cylindrical_provider_t), intent(in) :: prov + real(dp), intent(in) :: x(3) + real(dp), intent(out) :: gfd(3,3,3) + real(dp) :: h, gp(3,3), gm(3,3), ginv(3,3), sqrtg + real(dp) :: dg(3,3,3), gnow(3,3), xp(3) + integer :: i, j, k, l + h = 1d-4 + + call prov%metric(x, gnow, ginv, sqrtg) + ! dg(k,i,j) = d g_ij / d x_k + do k = 1, 3 + xp = x; xp(k) = x(k) + h + call prov%metric(xp, gp, ginv, sqrtg) + xp = x; xp(k) = x(k) - h + call prov%metric(xp, gm, ginv, sqrtg) + dg(k,:,:) = (gp - gm) / (2.0_dp*h) + end do + call prov%metric(x, gnow, ginv, sqrtg) + + gfd = 0.0_dp + do i = 1, 3 + do j = 1, 3 ! j = m + do k = 1, 3 ! k = n + do l = 1, 3 + gfd(i,j,k) = gfd(i,j,k) + 0.5_dp*ginv(i,l) * & + (dg(j,l,k) + dg(k,l,j) - dg(l,j,k)) + end do + end do + end do + end do + end subroutine christoffel_fd + +end program test_full_orbit diff --git a/test/tests/test_orbit_model_dispatch.f90 b/test/tests/test_orbit_model_dispatch.f90 new file mode 100644 index 00000000..00ba1c7a --- /dev/null +++ b/test/tests/test_orbit_model_dispatch.f90 @@ -0,0 +1,66 @@ +program test_orbit_model_dispatch + ! Wave-1 followup #1: orbit_model is read from the config namelist and the + ! macrostep dispatch maps it to the right pusher. This test writes a minimal + ! namelist, parses it via params%read_config, and asserts: + ! - orbit_model is parsed (default 0 = GC; here set to ORBIT_PAULI). + ! - cpp_stages_from_mode maps GAUSS1..4 to stage counts 1..4 (the CPP branch + ! dispatch key), proving the integer-coded select-case path is wired. + use, intrinsic :: iso_fortran_env, only: dp => real64 + use params, only: orbit_model, integmode, read_config + use orbit_full, only: ORBIT_GC, ORBIT_PAULI, ORBIT_BORIS, ORBIT_FOSYMPL + use orbit_symplectic_base, only: GAUSS1, GAUSS2, GAUSS3, GAUSS4 + use orbit_cpp, only: cpp_stages_from_mode + + implicit none + + integer :: nfail, u + character(256) :: cfgfile + + nfail = 0 + cfgfile = 'test_orbit_model_dispatch.in' + + open (newunit=u, file=cfgfile, status='replace', action='write') + write (u, '(A)') '&config' + write (u, '(A)') ' orbit_model = 1' + write (u, '(A)') ' integmode = 4' + write (u, '(A)') '/' + close (u) + + call read_config(cfgfile) + + call check('orbit_model parsed as ORBIT_PAULI', orbit_model == ORBIT_PAULI, nfail) + call check('integmode parsed as GAUSS1', integmode == GAUSS1, nfail) + + ! The dispatch keys are distinct integers (no overlap). + call check('orbit model codes distinct', & + ORBIT_GC == 0 .and. ORBIT_PAULI == 1 .and. ORBIT_BORIS == 2 .and. & + ORBIT_FOSYMPL == 3, nfail) + + ! Stage mapping that the CPP select-case dispatch uses. + call check('GAUSS1 -> 1 stage', cpp_stages_from_mode(GAUSS1) == 1, nfail) + call check('GAUSS2 -> 2 stages', cpp_stages_from_mode(GAUSS2) == 2, nfail) + call check('GAUSS3 -> 3 stages', cpp_stages_from_mode(GAUSS3) == 3, nfail) + call check('GAUSS4 -> 4 stages', cpp_stages_from_mode(GAUSS4) == 4, nfail) + + if (nfail == 0) then + print *, 'ALL ORBIT-MODEL DISPATCH TESTS PASSED' + else + print *, 'ORBIT-MODEL DISPATCH TESTS FAILED: ', nfail + error stop 1 + end if + +contains + + subroutine check(name, ok, nfail) + character(*), intent(in) :: name + logical, intent(in) :: ok + integer, intent(inout) :: nfail + if (ok) then + print '(A,A)', 'PASS ', name + else + print '(A,A)', 'FAIL ', name + nfail = nfail + 1 + end if + end subroutine check + +end program test_orbit_model_dispatch