From 05164e7e973021e8ec9ee411757d39ec2f242b7e Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 4 Aug 2026 10:57:27 -0700 Subject: [PATCH 1/2] remove gpu code --- orbitize/__init__.py | 10 ----- orbitize/kepler.py | 84 +++++-------------------------------- orbitize/kernels/mikkola.cu | 45 -------------------- orbitize/kernels/newton.cu | 50 ---------------------- tests/test_kepler_solver.py | 45 +++++--------------- 5 files changed, 20 insertions(+), 214 deletions(-) delete mode 100644 orbitize/kernels/mikkola.cu delete mode 100644 orbitize/kernels/newton.cu diff --git a/orbitize/__init__.py b/orbitize/__init__.py index 3c3475fa..ec57693e 100644 --- a/orbitize/__init__.py +++ b/orbitize/__init__.py @@ -6,16 +6,6 @@ orbitize_dir = os.path.dirname(__file__) DATADIR = os.path.join(orbitize_dir, "example_data/") -# Detect a valid CUDA environment -try: - import pycuda.driver as cuda - import pycuda.autoinit - from pycuda.compiler import SourceModule - - cuda_ext = True -except: - cuda_ext = False - try: from . import _kepler diff --git a/orbitize/kepler.py b/orbitize/kepler.py index c7506686..bd5e7cd4 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -10,10 +10,6 @@ if cext: from . import _kepler -if cuda_ext: - # Configure GPU context for CUDA accelerated compute - from orbitize import gpu_context - kep_gpu_ctx = gpu_context.gpu_context() def tau_to_manom(date, sma, mtot, tau, tau_ref_epoch): """ @@ -54,7 +50,6 @@ def times2trueanom_and_eccanom( tolerance=1e-9, max_iter=100, use_c=True, - use_gpu=False, ): """ Convert times to true anomaly and eccentric anomaly by solving Kepler's Equation. @@ -69,7 +64,6 @@ def times2trueanom_and_eccanom( tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. use_c (bool, optional): Use the C solver if configured. Defaults to True - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False Returns: 2-tuple: @@ -91,7 +85,7 @@ def times2trueanom_and_eccanom( # # compute mean anomaly (size: n_orbs x n_dates) manom = tau_to_manom(epochs[:, None], sma, mtot, tau, tau_ref_epoch) # compute eccentric anomalies (size: n_orbs x n_dates) - eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) + eanom = _calc_ecc_anom(manom, ecc_arr, tolerance=tolerance, max_iter=max_iter, use_c=use_c) # compute the true anomalies (size: n_orbs x n_dates) # Note: matrix multiplication makes the shapes work out here and below @@ -103,7 +97,7 @@ def times2trueanom_and_eccanom( def calc_orbit( epochs, sma, ecc, inc, aop, pan, tau, plx, mtot, mass_for_Kamp=None, tau_ref_epoch=58849, tolerance=1e-9, - max_iter=100, use_c=True, use_gpu=False + max_iter=100, use_c=True ): """ @@ -130,7 +124,6 @@ def calc_orbit( tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. use_c (bool, optional): Use the C solver if configured. Defaults to True - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False Return: 3-tuple: @@ -151,7 +144,7 @@ def calc_orbit( mass_for_Kamp = mtot ecc - tanom, eanom = times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=tau_ref_epoch, tolerance=tolerance, max_iter=max_iter, use_c=use_c, use_gpu=use_gpu) + tanom, eanom = times2trueanom_and_eccanom(sma, epochs, mtot, ecc, tau, tau_ref_epoch=tau_ref_epoch, tolerance=tolerance, max_iter=max_iter, use_c=use_c) # compute 3-D orbital radius of second body (size: n_orbs x n_dates) radius = sma * (1.0 - ecc * np.cos(eanom)) @@ -184,7 +177,7 @@ def calc_orbit( vz = np.squeeze(vz)[()] return raoff, deoff, vz -def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gpu=False): +def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False): """ Computes the eccentric anomaly from the mean anomlay. Code from Rob De Rosa's orbit solver (e < 0.95 use Newton, e >= 0.95 use Mikkola) @@ -195,7 +188,6 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp tolerance (float, optional): absolute tolerance of iterative computation. Defaults to 1e-9. max_iter (int, optional): maximum number of iterations before switching. Defaults to 100. use_c (bool, optional): Use the C solver if configured. Defaults to False - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False Return: eanom (float/np.array): eccentric anomalies, same shape as manom @@ -231,16 +223,16 @@ def _calc_ecc_anom(manom, ecc, tolerance=1e-9, max_iter=100, use_c=False, use_gp # Now low eccentricities ind_low = np.where(~ecc_zero & ecc_low) if len(ind_low[0]) > 0: - eanom[ind_low] = _newton_solver_wrapper(manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c, use_gpu) + eanom[ind_low] = _newton_solver_wrapper(manom[ind_low], ecc[ind_low], tolerance, max_iter, use_c) # Now high eccentricities ind_high = np.where(~ecc_zero & ~ecc_low | (eanom == -1)) # The C and CUDA solvers return the unphysical value -1 if they fail to converge if len(ind_high[0]) > 0: - eanom[ind_high] = _mikkola_solver_wrapper(manom[ind_high], ecc[ind_high], use_c, use_gpu) + eanom[ind_high] = _mikkola_solver_wrapper(manom[ind_high], ecc[ind_high], use_c) return np.squeeze(eanom)[()] -def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu=False): +def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False): """ Wrapper for the various (Python, C, CUDA) implementations of the Newton-Raphson solver for eccentric anomaly. @@ -250,7 +242,6 @@ def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu ecc (np.array): array of eccentricities eanom0 (np.array, optional): array of first guess for eccentric anomaly, same shape as manom (optional) use_c (bool, optional): Use the C solver if configured. Defaults to False - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False Return: eanom (np.array): array of eccentric anomalies @@ -258,10 +249,7 @@ def _newton_solver_wrapper(manom, ecc, tolerance, max_iter, use_c=False, use_gpu """ eanom = np.empty_like(manom) - if cuda_ext and use_gpu: - # the CUDA solver returns eanom = -1 if it doesnt converge after max_iter iterations - eanom = _CUDA_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) - elif cext and use_c: + if cext and use_c: # the C solver returns eanom = -1 if it doesnt converge after max_iter iterations eanom = _kepler._c_newton_solver(manom, ecc, tolerance=tolerance, max_iter=max_iter) else: @@ -323,33 +311,8 @@ def _newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): return eanom -def _CUDA_newton_solver(manom, ecc, tolerance=1e-9, max_iter=100, eanom0=None): - """ - Helper function for calling the CUDA implementation of the Newton-Raphson solver for eccentric anomaly. - - Args: - manom (np.array): array of mean anomalies - ecc (np.array): array of eccentricities - eanom0 (np.array, optional): array of first guess for eccentric anomaly, same shape as manom (optional) - Return: - eanom (np.array): array of eccentric anomalies - - Written: Devin Cody, 2021 - """ - global kep_gpu_ctx - - # Ensure manom and ecc are np.array (might get passed as astropy.Table Columns instead) - manom = np.asarray(manom) - ecc = np.asarray(ecc) - eanom = np.empty_like(manom) - tolerance = np.asarray(tolerance, dtype = np.float64) - max_iter = np.asarray(max_iter) - - kep_gpu_ctx.newton(manom, ecc, eanom, eanom0, tolerance, max_iter) - return eanom - -def _mikkola_solver_wrapper(manom, ecc, use_c=False, use_gpu=False): +def _mikkola_solver_wrapper(manom, ecc, use_c=False): """ Wrapper for the various (Python, C, CUDA) implementations of Analtyical Mikkola solver @@ -357,8 +320,6 @@ def _mikkola_solver_wrapper(manom, ecc, use_c=False, use_gpu=False): manom (np.array): array of mean anomalies between 0 and 2pi ecc (np.array): eccentricity use_c (bool, optional): Use the C solver if configured. Defaults to False - use_gpu (bool, optional): Use the GPU solver if configured. Defaults to False - Return: eanom (np.array): array of eccentric anomalies @@ -368,9 +329,7 @@ def _mikkola_solver_wrapper(manom, ecc, use_c=False, use_gpu=False): ind_change = np.where(manom > np.pi) manom[ind_change] = (2.0 * np.pi) - manom[ind_change] - if cuda_ext and use_gpu: - eanom = _CUDA_mikkola_solver(manom, ecc) - elif cext and use_c: + if cext and use_c: eanom = _kepler._c_mikkola_solver(manom, ecc) else: eanom = _mikkola_solver(manom, ecc) @@ -417,26 +376,3 @@ def _mikkola_solver(manom, ecc): u4 = -f/(f1+0.5*f2*u3+(1.0/6.0)*f3*u3*u3+(1.0/24.0)*f4*(u3**3.0)) return (e0 + u4) - -def _CUDA_mikkola_solver(manom, ecc): - """ - Helper function for calling the CUDA implementation of the Analtyical Mikkola solver for the eccentric anomaly. - - Args: - manom (float or np.array): mean anomaly, must be between 0 and pi. - ecc (float or np.array): eccentricity - Return: - eanom (np.array): array of eccentric anomalies - - Written: Devin Cody, 2021 - """ - global kep_gpu_ctx - - # Ensure manom and ecc are np.array (might get passed as astropy.Table Columns instead) - manom = np.asarray(manom) - ecc = np.asarray(ecc) - eanom = np.empty_like(manom) - - kep_gpu_ctx.mikkola(manom, ecc, eanom) - - return eanom diff --git a/orbitize/kernels/mikkola.cu b/orbitize/kernels/mikkola.cu deleted file mode 100644 index 3eb940c0..00000000 --- a/orbitize/kernels/mikkola.cu +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef M_PI -#define M_PI 3.14159265358979323846 /* pi */ -#endif - -__global__ void mikkola_gpu(const double *manom, const double *ecc, double *eanom){ - /* - Vectorized C Analtyical Mikkola solver for the eccentric anomaly. - See: S. Mikkola. 1987. Celestial Mechanics, 40, 329-334. - Adapted from IDL routine keplereq.pro by Rob De Rosa http://www.lpl.arizona.edu/~bjackson/idl_code/keplereq.pro - Args: - manom (double[]): mean anomaly, must be between 0 and pi. - ecc (double[]): eccentricity - eanom0 (double[]): array for eccentric anomaly - Return: - None: eanom (double[]): is changed by reference - Written: Devin Cody, 2019 - */ - - int i = threadIdx.x + blockIdx.x*blockDim.x; - double alpha, beta, aux, z, s0, s1, se0, ce0; - double f, f1, f2, f3, f4, u1, u2, u3; - - alpha = (1.0 - ecc[i]) / ((4.0 * ecc[i]) + 0.5); - beta = (0.5 * manom[i]) / ((4.0 * ecc[i]) + 0.5); - - aux = sqrt(beta*beta + alpha*alpha*alpha); - z = pow(fabs(beta + aux), (1.0/3.0)); - - s0 = z - (alpha/z); - s1 = s0 - (0.078*(pow(s0, 5))) / (1.0 + ecc[i]); - eanom[i] = manom[i] + (ecc[i] * (3.0*s1 - 4.0*(s1*s1*s1))); - - se0=sin(eanom[i]); - ce0=cos(eanom[i]); - - f = eanom[i]-ecc[i]*se0-manom[i]; - f1 = 1.0-ecc[i]*ce0; - f2 = ecc[i]*se0; - f3 = ecc[i]*ce0; - f4 = -f2; - u1 = -f/f1; - u2 = -f/(f1+0.5*f2*u1); - u3 = -f/(f1+0.5*f2*u2+(1.0/6.0)*f3*u2*u2); - eanom[i] += -f/(f1+0.5*f2*u3+(1.0/6.0)*f3*u3*u3+(1.0/24.0)*f4*(u3*u3*u3)); -} \ No newline at end of file diff --git a/orbitize/kernels/newton.cu b/orbitize/kernels/newton.cu deleted file mode 100644 index 3f8667bb..00000000 --- a/orbitize/kernels/newton.cu +++ /dev/null @@ -1,50 +0,0 @@ - -#ifndef M_PI -#define M_PI 3.14159265358979323846 /* pi */ -#endif - - -__global__ void newton_gpu(const double *manom, - const double *ecc, - double *eanom, - const int *max_iter, - const double *tol){ - /* - Vectorized C++ Newton-Raphson solver for eccentric anomaly. - Args: - manom (double[]): array of mean anomalies - ecc (double[]): array of eccentricities - eanom0 (double[]): array of first guess for eccentric anomaly, same shape as manom (optional) - Return: - None: eanom is changed by reference - Written: Devin Cody, 2018 - */ - int i = threadIdx.x + blockIdx.x*blockDim.x; - double diff; - int niter = 0; - int half_max = *max_iter/2.0; // divide convergence->max_iter by 2 using bit shift - - // Let's do one iteration to start with - eanom[i] -= (eanom[i] - (ecc[i] * sin(eanom[i])) - manom[i]) / (1.0 - (ecc[i] * cos(eanom[i]))); - diff = (eanom[i] - (ecc[i] * sin(eanom[i])) - manom[i]) / (1.0 - (ecc[i] * cos(eanom[i]))); - - while ((fabs(diff) > *tol) && (niter <= *max_iter)){ - eanom[i] -= diff; - - // If it hasn't converged after half the iterations are done, try starting from pi - if (niter == half_max) { - eanom[i] = M_PI; - } - - diff = (eanom[i] - (ecc[i] * sin(eanom[i])) - manom[i]) / (1.0 - (ecc[i] * cos(eanom[i]))); - niter += 1; - } - - // If it has not converged, set eccentricity to -1 to signal that it needs to be - // solved using the analytical version. Note this behavior is a bit different from the - // numpy implementation - if (niter >= *max_iter){ - printf("%f %f %f %f >= %d iter\n", manom[i], eanom[i], diff, ecc[i], *max_iter); - eanom[i] = -1; - } -} diff --git a/tests/test_kepler_solver.py b/tests/test_kepler_solver.py index fa234ced..4410345e 100644 --- a/tests/test_kepler_solver.py +++ b/tests/test_kepler_solver.py @@ -8,7 +8,6 @@ import os import numpy as np import orbitize.kepler as kepler -from orbitize import cuda_ext from orbitize import cext threshold = 1e-5 @@ -17,7 +16,7 @@ def angle_diff(ang1, ang2): # Return the difference between two angles return np.arctan2(np.sin(ang1 - ang2), np.cos(ang1 - ang2)) -def test_analytical_ecc_anom_solver(use_c = False, use_gpu = False): +def test_analytical_ecc_anom_solver(use_c = False): """ Test orbitize.kepler._calc_ecc_anom() in the analytical solver regime (e > 0.95) by comparing the mean anomaly computed from _calc_ecc_anom() output vs the input mean anomaly @@ -25,12 +24,12 @@ def test_analytical_ecc_anom_solver(use_c = False, use_gpu = False): mean_anoms = np.linspace(0,2.0*np.pi,1000) eccs = np.linspace(0.95,0.999999,100) for ee in eccs: - ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c, use_gpu = use_gpu) + ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c) calc_mm = (ecc_anoms - ee*np.sin(ecc_anoms)) % (2*np.pi) # plug solutions into Kepler's equation for meas, truth in zip(calc_mm, mean_anoms): assert angle_diff(meas, truth) == pytest.approx(0.0, abs=threshold) -def test_iterative_ecc_anom_solver(use_c = False, use_gpu = False): +def test_iterative_ecc_anom_solver(use_c = False): """ Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from _calc_ecc_anom() output vs the input mean anomaly @@ -38,7 +37,7 @@ def test_iterative_ecc_anom_solver(use_c = False, use_gpu = False): mean_anoms = np.linspace(0,2.0*np.pi,100) eccs = np.linspace(0,0.9499999,100) for ee in eccs: - ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c, use_gpu = use_gpu) + ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c=use_c) calc_ma = (ecc_anoms - ee*np.sin(ecc_anoms)) % (2*np.pi) # plug solutions into Kepler's equation for meas, truth in zip(calc_ma, mean_anoms): assert angle_diff(meas, truth) == pytest.approx(0.0, abs=threshold) @@ -52,11 +51,6 @@ def test_c_ecc_anom_solver(): test_iterative_ecc_anom_solver(use_c = True) test_analytical_ecc_anom_solver(use_c = True) -def test_pycuda_ecc_anom_solver(): - if cuda_ext: - test_iterative_ecc_anom_solver(use_gpu = True) - test_analytical_ecc_anom_solver(use_gpu = True) - def test_orbit_e03(): @@ -233,7 +227,7 @@ def test_orbit_scalar(): assert true_deoff == pytest.approx(deoffs, abs=threshold) assert true_vz == pytest.approx(vzs, abs=1e-8) -def profile_iterative_ecc_anom_solver(n_orbits = 1000, use_c = True, use_gpu = False): +def profile_iterative_ecc_anom_solver(n_orbits = 1000, use_c = True): """ Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from _calc_ecc_anom() output vs the input mean anomaly @@ -242,9 +236,9 @@ def profile_iterative_ecc_anom_solver(n_orbits = 1000, use_c = True, use_gpu = F mean_anoms=np.linspace(0, 2.0*np.pi,n_orbits) eccs=np.linspace(0,0.9499999, n_orbits) for ee in eccs: - ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c = use_c, use_gpu = use_gpu) + ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, tolerance=1e-9, use_c = use_c) -def profile_mikkola_ecc_anom_solver(n_orbits = 1000, use_c = True, use_gpu = False): +def profile_mikkola_ecc_anom_solver(n_orbits = 1000, use_c = True): """ Test orbitize.kepler._calc_ecc_anom() in the iterative solver regime (e < 0.95) by comparing the mean anomaly computed from _calc_ecc_anom() output vs the input mean anomaly @@ -252,22 +246,12 @@ def profile_mikkola_ecc_anom_solver(n_orbits = 1000, use_c = True, use_gpu = Fal mean_anoms=np.linspace(0, 2.0*np.pi,n_orbits) eccs=np.linspace(.95,0.999999, n_orbits) for ee in eccs: - ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, use_c = use_c, use_gpu = use_gpu) + ecc_anoms = kepler._calc_ecc_anom(mean_anoms, ee, use_c = use_c) def profile_all(n_orbits, print_profiles = False): profile_name = "Profile.prof" n_print_lines = 15 d = dict() - - if cuda_ext: - cProfile.runctx("profile_iterative_ecc_anom_solver(n_orbits = n_orbits, use_c = False, use_gpu = True)", globals(), locals(), profile_name) - s = pstats.Stats(profile_name) - if print_profiles: - print("Profiling Newton: CUDA with {} orbits".format(n_orbits**2)) - s.strip_dirs().sort_stats("time").print_stats(n_print_lines) - d["Newton GPU Solver"] = s.__dict__["total_tt"] - else: - print("System not configured for CUDA") if cext: cProfile.runctx("profile_iterative_ecc_anom_solver(n_orbits = n_orbits, use_c = True)", globals(), locals(), profile_name) @@ -286,23 +270,15 @@ def profile_all(n_orbits, print_profiles = False): s.strip_dirs().sort_stats("time").print_stats(n_print_lines) d["Newton Python Solver"] = s.__dict__["total_tt"] - if cuda_ext: - cProfile.runctx("profile_mikkola_ecc_anom_solver(n_orbits = n_orbits, use_c = False, use_gpu = True)", globals(), locals(), profile_name) - s = pstats.Stats(profile_name) - if print_profiles: - print("Profiling Mikkola: CUDA with {} orbits".format(n_orbits**2)) - s.strip_dirs().sort_stats("time").print_stats(n_print_lines) - d["Mikkola GPU Solver"] = s.__dict__["total_tt"] - if cext: - cProfile.runctx("profile_mikkola_ecc_anom_solver(n_orbits = n_orbits, use_c = True, use_gpu = False)", globals(), locals(), profile_name) + cProfile.runctx("profile_mikkola_ecc_anom_solver(n_orbits = n_orbits, use_c = True)", globals(), locals(), profile_name) s = pstats.Stats(profile_name) if print_profiles: print("Profiling Mikkola: C with {} orbits".format(n_orbits**2)) s.strip_dirs().sort_stats("time").print_stats(n_print_lines) d["Mikkola C Solver"] = s.__dict__["total_tt"] - cProfile.runctx("profile_mikkola_ecc_anom_solver(n_orbits = n_orbits, use_c = False, use_gpu = False)", globals(), locals(), profile_name) + cProfile.runctx("profile_mikkola_ecc_anom_solver(n_orbits = n_orbits, use_c = False)", globals(), locals(), profile_name) s = pstats.Stats(profile_name) if print_profiles: print("Profiling Mikkola: Python with {} orbits".format(n_orbits**2)) @@ -327,7 +303,6 @@ def profile_all(n_orbits, print_profiles = False): test_analytical_ecc_anom_solver() test_iterative_ecc_anom_solver() test_c_ecc_anom_solver() - test_pycuda_ecc_anom_solver() test_orbit_e03() test_orbit_e03_array() test_orbit_e99() From b4f14e37c9ae8b52d33744134f7d1d7d901a52c6 Mon Sep 17 00:00:00 2001 From: Sarah Blunt Date: Tue, 4 Aug 2026 10:59:30 -0700 Subject: [PATCH 2/2] remove cuda import --- orbitize/kepler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/kepler.py b/orbitize/kepler.py index bd5e7cd4..3365bc3a 100644 --- a/orbitize/kepler.py +++ b/orbitize/kepler.py @@ -5,7 +5,7 @@ import astropy.units as u import astropy.constants as consts -from orbitize import cuda_ext, cext +from orbitize import cext if cext: from . import _kepler