experiment: claude ported everything - #82
Draft
JorgeG94 wants to merge 46 commits into
Draft
Conversation
- Block all local arrays in find_N2 with an nj j-dimension; move the find_N2 call before the inner OMP j-loop in set_diffusivity and add an internal j-loop inside find_N2 - Make h_bot/k_bot non-optional in find_N2 - Add find_rho_bottom generic interface in MOM_interface_heights, renaming the existing routine to find_rho_bottom_1d and adding find_rho_bottom_2d that accepts blocked arrays (nj j-dimension); move find_rho_bottom call outside find_N2's j-loop - Replace each outer j-loop in find_N2 and find_rho_bottom_2d with individual do j / do i nested loops placed next to each i-loop Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…_to_Kd Add njblock j-dimension to local and dummy arrays in calculate_bkgnd_mixing and find_TKE_to_Kd, move j-loops into both routines, and decouple inner loops to individual j-loops in the caller. Also block TKE_to_Kd, maxTKE, and kb arrays in set_diffusivity. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Compute dz once unconditionally in set_diffusivity and pass it as a required argument to find_N2, find_TKE_to_Kd, and add_LOTW_BBL_diffusivity, removing the redundant per-call thickness_to_dz invocations in each routine. set_BBL_TKE retains its own thickness_to_dz since it is called externally. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…P loop Give each branch its own j-loop outside the OMP parallel region, with the int_tides diagnostics kept inside that loop alongside the call. Remove stray print statements. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Add nj j-block dimension to local arrays (htot, rho_htot, gh_sum_top, Rho_top, TKE, I2decay, do_i), promote all dummy array args to blocked, replace scalar j with jstart/jend/nj, and fission the outer j-loop into three separate ji/jki loop nests. Move the call before the OMP region and pass full blocked arrays. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Moves add_LOTW_BBL_diffusivity, VBF%Kd_BBL copy, limit_dissipation, Kd_add, Kd_int/Kd_lay copy, Kd_Work, and Kd_Work_added blocks out of the OMP parallel do into individual do j=jstart,jend loops. Removes the now-empty OMP parallel do directive and its enclosing j-loop. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Record each clock's name as cpu_clock_id() registers it, then open an NVTX range in cpu_clock_begin() and close it in cpu_clock_end(). Every existing cpu_clock_id() name becomes a named range in an nsys timeline with no call-site changes. The range is opened before, and closed after, the mpp clock so that it encloses it. Begin and end test the same condition, so starts and ends stay balanced for handles that were never named or that fall outside the table. The code is guarded behind -DMOM_USE_NVTX, which is undefined by default: without it the preprocessed source is byte-identical to before and no extra library is linked. Profiling builds need -DMOM_USE_NVTX and -cudalib=nvtx. Applied to both the FMS1 and FMS2 infra shims, which are byte-identical to each other. This also covers the halo and message-passing clocks, which MOM_domain_infra begins and ends through this module directly rather than through the MOM_cpu_clock framework wrapper. Based on Edward Yang's ae67665 ("add nvtx markers to clocks", branch benchmark_ALE_nvtx_clocks), which introduced the name-table approach and the range nesting; this adds the build guard so non-NVHPC builds are unaffected. Co-authored-by: Edward Yang <edward.yang@anu.edu.au> Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
First increment of the MOM_kappa_shear GPU port. Converts the two
thickness-to-vertex interpolation loops in Calc_kappa_shear_vertex
(both branches of vertex_shear_OBC_bug) from !$OMP parallel do to
do concurrent, offloading them to the device.
Residency:
- h is host-authoritative in the (host-only) diabatic stack and is
mutated on host during diabatic (applyBoundaryFluxesInOut), so the
globally-resident device copy (CS%h) is refreshed with update to(h)
before the kernels; a plain map(to:) does not refresh a present
object.
- G%mask2dCu/Cv/T are already device-resident (initialize_MOM).
- h_at_u/h_at_v are device workspace (map alloc), copied back with
update from for the still-host per-column solver below. enter/exit
data are mirrored and refcount-neutral (release, not delete).
The per-column solver (kappa_shear_column/find_kappa_tke) stays on
host for now; later increments move it to the device, at which point
the update from(h_at_u,h_at_v) round-trip is removed.
Verification:
- benchmark_ALE ocean.stats bit-identical to the dev/gpu reference (np=1).
- nsys --stats=true confirms both do concurrent loops run on the GPU
(calc_kappa_shear_vertex_{550,553}_gpu, 12 launches each).
- Reviewed; no code defect found. The >=2-rank bitwise run remains the
mandatory pre-merge gate and is still pending (single-GPU node).
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Groundwork for whole-column GPU kernels (kappa_shear next) that must
compute density derivatives from INSIDE a device region, where the
polymorphic EOS%type%... dispatch is not usable.
- Expose the existing this-free _loc elementals
calculate_density_derivs_elem_{Roquet_rho,buggy_Wright}_loc as public
and !$omp declare target.
- Add calculate_density_derivs_elem_loc(form_of_EOS, T, S, pressure,
drho_dT, drho_dS): a declare-target subroutine that select-case's on
the integer form and calls the matching _loc kernel, returning mks
derivatives. Unit rescaling and any `scale` factor remain the
caller's responsibility (as in calculate_density_derivs_1d). Forms
without a device _loc kernel hit the default branch; a device-using
module must FATAL at init on a GPU build before reaching it.
This is the reusable, form-general enabler chosen over a Roquet-scoped
inline; it is aligned with doc 06 6.4's select-case direction but adds
only the device-callable point path (not the full polymorphism removal).
No functional change: nothing calls the dispatcher yet.
Verification: benchmark_ALE ocean.stats bit-identical to the dev/gpu
reference (np=1). Dispatcher is exercised (and Fable-reviewed) in the
follow-up that wires kappa_shear to it.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…r on GPU builds Wires kappa_shear_column's Boussinesq density-derivative computation to the device-callable path added in 47ac931, so the column solver can later run inside a target region (the polymorphic calculate_density_derivs interface is not device-callable). - MOM_EOS: add get_EOS_form_and_scaling(EOS, form, kg_m3_to_R, C_to_degC, S_to_ppt, RL2_T2_to_Pa) so callers can reproduce, host-side, the unit rescaling that calculate_density_derivs_1d applies (EOS_type components are private). - kappa_shear_column: under #ifdef __NVCOMPILER_OPENMP_GPU, replace the calculate_density_derivs array call with an inline K=2..nzc loop calling calculate_density_derivs_elem_loc + the exact scale math of calculate_density_derivs_1d(dom=(/2,nzc/), scale=-g_R0). Non-GPU builds keep the original polymorphic call verbatim for all EOS forms. On GPU builds only ROQUET_RHO and WRIGHT have device _loc kernels; others FATAL. Verification: benchmark_ALE (ROQUET_RHO, Boussinesq) ocean.stats bit-identical to the dev/gpu reference (np=1) with the GPU path active on host. Reviewed; bitwise equivalence confirmed operation-by-operation, verdict ship. Follow-ups required by the offload step (increment 3), per review: hoist get_EOS_form_and_scaling + the form FATAL out of kappa_shear_column to the driver (both are host-only and cannot run in a device region), passing the form/scaling in; and FATAL on non-Boussinesq on GPU builds (that branch still uses polymorphic EOS calls). Untested-by-benchmark but code-identical: the non-unity unit-rescale branch and the WRIGHT dispatcher case. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…d prep) Prerequisite for offloading the column solver (increment 3): get_EOS_form_and_scaling and MOM_error are host-only and cannot run inside a device region, so move them out of kappa_shear_column up into both drivers (Calculate_kappa_shear and Calc_kappa_shear_vertex). - Resolve eos_form + the four unit-scaling factors once per driver call (guarded by use_temperature), and thread them into kappa_shear_column as new intent(in) args. - Move the "unsupported EOS form on GPU" FATAL to the drivers, and add a FATAL on non-Boussinesq on GPU builds (that branch still uses polymorphic EOS calls). - Add the five hoisted scalars to the shared() clause of both drivers' !$OMP parallel do (they are loop-invariant and read-only in the region; default(private) would leave them uninitialized inside the loop -- caught by the bitwise check). No functional change: same values, computed once per call instead of once per column. Verification: benchmark_ALE (vertex driver, ROQUET_RHO, Boussinesq) ocean.stats bit-identical to the dev/gpu reference (np=1). Reviewed; verdict ship, bitwise safety confirmed by inspection for both drivers. Not runtime-exercised by the benchmark (inspection-only): the h-point driver, the non-unity unit-scaling branch. Notes for the offload step: the driver FATAL guards execution not compilation, so declare target on kappa_shear_column will require #ifdef-excluding the non-Boussinesq branch; CS (pointer dummy) and tv (pointer-laden) must be demoted/pre-extracted before declare target. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Prepares the per-column solver for offload: the three workers (kappa_shear_column, find_kappa_tke, calculate_projected_state) are now !$omp declare target, so the driver column loop can run inside a target region in the next step. Device code is generated for all three (confirmed by -Minfo). Execution is still on the host (the column loop is still !$OMP parallel do), so this is bitwise-neutral. To make them device-compilable: - Demote CS from `pointer` to `type(Kappa_shear_CS), intent(in)` in kappa_shear_column and find_kappa_tke (the type is all scalars bar a diag pointer that the column math never touches). - Hoist use_temperature to a kappa_shear_column argument (resolved in the drivers) so the device path no longer reads the tv%T pointer via associated(); tv is then untouched on GPU builds. - #ifndef __NVCOMPILER_OPENMP_GPU around the non-Boussinesq branch's polymorphic calculate_specific_vol_derivs/calculate_density calls so they are excluded from device compilation (unreachable on GPU anyway -- the driver FATALs on non-Boussinesq). Notably the 26 early exit/return sites in the two solvers compile for the device without error on NVHPC 26.5 (the 25.11 exit-under-parallel miscompile appears fixed); their runtime correctness on device is gated by the bitwise check in the offload step, not pre-emptive rewrites. Verification: benchmark_ALE ocean.stats bit-identical to the dev/gpu reference (np=1). Next step flips the driver !$OMP parallel do to a target teams loop and adds residency; the solver then runs on the GPU. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Promote the per-J 2-D vertex slabs (u_2d..rho_2d) in Calc_kappa_shear_vertex to 3-D arrays (u_slab..rho_slab) computed in device do-concurrent passes ahead of the column loop. The loop bodies are verbatim except that I_hwt is inlined as a reciprocal multiply and the use_temperature branch is hoisted so each pass is purely elementwise; answers are bitwise-identical (verified np=1 against repro_ocean.stats on benchmark_ALE, V100/NVHPC 26.5). The key residency fix relative to the earlier reverted attempt: u_in/v_in and T_in/S_in (tv%T, tv%S) are already device-resident (mapped in MOM.F90), so a map(to:) on them does not copy - and the diabatic stack modifies T/S on the host before set_diffusivity, leaving the device copies stale. They now get an explicit "target update to" refresh, which resolves the NaN-diffusivity crash the first attempt hit. dz_3d is a fresh local each call, so its map(to:) does copy. h_at_u/h_at_v are now consumed on the device by the slab passes, so increment 1's copy-back is dropped; the slabs are copied back for the (still host) per-column solver loop, which reads them as X_slab(I,J,k). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Convert Calc_kappa_shear_vertex's per-column loop from a host OMP parallel do over J into an "omp target teams loop collapse(2)" over (J,I), calling the declare-target kappa_shear_column on the device. Bitwise-identical (np=1, benchmark_ALE vs repro_ocean.stats, V100/NVHPC 26.5), which also demonstrates that the solver's 26 early exit/return sites are runtime-correct on 26.5. Restructuring, all bitwise-neutral: - The per-J staging slices kappa_2d/tke_2d become device-resident 3-D staging arrays kappa_3d/tke_3d (map(alloc:)), written per column and consumed by new elementwise do-concurrent write-back passes for kappa_vertex/tke_io/kv_io. - The surface pressure is interpolated to vertices on the host ahead of the region (surface_pres_2d, verbatim expressions), keeping the possibly-null p_surf pointer and the CS%psurf_bug branch out of device code. - The h_vert fill stays on the host (it is only consumed by the host tracer-point averaging when VS_ThicknessMean=true), fed by a guarded update from(h_slab). - All per-column scratch is private at the loop; loop-invariant scalars are firstprivate; CS gets a per-call shallow map(to:) in the driver (all scalars plus a diag pointer that device code never dereferences). kappa_vertex, tke_io and kv_io are mapped to: (not alloc:) so host-set halo/initial values survive the full-array update from. diag_N2/S2_init/mean are mapped after their conditional host zeroing and copied back under one guard matching their consumers (post_data/Bchksum). Known limitation, to be fixed next: the solver's automatic arrays (sized nzc or SZK_(GV)) are device-heap allocated, so runs currently need NVCOMPILER_ACC_CUDA_HEAPSIZE=1073741824 (the default heap is too small: "DEVICE FORTRAN AUTO ALLOCATION FAILED"). Fixed-size locals will remove this. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… allocs)
The vertex driver's per-column private scratch and the column routines'
automatic locals are sized by GV%ke/nzc, which is not a compile-time constant,
so nvfortran allocated them from the device heap on every column call
(NVFORTRAN-W-0155). That exhausted the default device heap ("DEVICE FORTRAN
AUTO ALLOCATION FAILED" unless NVCOMPILER_ACC_CUDA_HEAPSIZE was raised) and
serialized columns on the device allocator (~2x total model slowdown).
In GPU builds (#ifdef __NVCOMPILER_OPENMP_GPU) these arrays now use a
compile-time ceiling, the new module parameter GPU_nk_max=128, so each device
thread gets stack ("local memory") arrays; kappa_shear_init fatals if
GV%ke > GPU_nk_max when the scheme is active. CPU builds keep the exact-size
declarations. No code queries these arrays' sizes and all accesses are
indexed to nzc/nz, so the enlargement is bitwise-neutral: verified np=1
against repro_ocean.stats at the default device heap, wall time back to the
pre-offload baseline (98s -> 45s), and the W-0155 warnings are gone.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…residency Review of the column-offload series found one latent halo bug: kv_io's actual argument (visc%Kv_shear_Bu) is already persistently device-resident, mapped in set_visc_init BEFORE the restart-reproducibility pass_var halo update, so this routine's map(to: kv_io) copies nothing and the full-array "update from" wrote stale pre-pass_var device halos over the corrected host halos - silently resurrecting the restart non-reproducibility in non-symmetric mode that the pass_var was added to fix. Invisible to a cold-start np=1 run. Fixed with an explicit "update to(tke_io, kv_io)" after the maps (tke_io defensively, so the routine does not depend on it staying unmapped elsewhere). Also corrects the increment-2 residency comment and drops the redundant "update to(T_in, S_in)": the actuals are the caller's convection-filtered T_f/S_f, fresh host locals each call, so their map(to:) already copies (they are not tv%T/tv%S as the comment claimed). This also establishes that the earlier reverted slab attempt's NaN was not T/S staleness; the operative fix was restructuring the interpolation (hoisting the use_temperature branch out of the do concurrent). And the GPU_nk_max FATAL message now notes it covers the shared tracer-point column routines too. Bit-identical at np=1 (benchmark_ALE vs repro_ocean.stats, V100/NVHPC 26.5). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Co-authored-by: Utheri Wagura <Utheri.Wagura@noaa.gov>
Groundwork for offloading the finite-volume pressure-gradient density integrals (int_density_dz_generic_plm next), which must evaluate the equation of state from INSIDE a device region, where the polymorphic EOS%type%density_elem / density_anomaly_elem dispatch is not usable. - Expose the existing this-free _loc elementals as public and !$omp declare target: density_elem_Roquet_rho_loc and density_anomaly_elem_Roquet_rho_loc (Roquet_rho), and density_elem_buggy_Wright_loc (Wright). buggy_Wright has no anomaly _loc kernel, so its anomaly branch stays unsupported here. - Add calculate_density_elem_loc(form_of_EOS, T, S, pressure, use_rho_ref, rho_ref): a declare-target function that select-case's on the integer form and returns, in mks units, either in-situ density or (when use_rho_ref) the anomaly relative to rho_ref, by forwarding to the matching _loc kernel. It is the density analogue of calculate_density_derivs_elem_loc (47ac931) and, like it, leaves unit rescaling (EOS%*_to_* factors) and any `scale` factor to the caller, exactly as calculate_density_elem does. Combinations without a device-callable kernel -- an unsupported form, or the anomaly branch of buggy_Wright -- return 0; a device-using module must FATAL at init on a GPU build before reaching such a combination. Bitwise-safety: this is purely additive. Nothing calls the dispatcher yet, and the exposed elementals are the identical functions the polymorphic host path already wraps (density_anomaly_elem_Roquet_rho => density_anomaly_elem_Roquet_rho_loc, etc.), so no existing arithmetic moves or changes. Scope note (reviewer-directed): no new Wright kernels are added here. benchmark_ALE (ROQUET_RHO) cannot exercise a new Wright anomaly kernel, and untested numerics is scope creep; parity can come later with a config that exercises it. Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical to the dev/gpu reference (np=1); DEBUG=True redundancy gate clean (0 'redundant'). "(Ocean pressure force)" clock unchanged at ~20.8 s/24 calls, as expected for an additive change. Co-Authored-By: Claude Opus <noreply@anthropic.com>
…fload prep) Host-side bitwise-safe refactor of the finite-volume PGF density integrals, preparing int_density_dz_generic_plm for offload (increment C next): replace the polymorphic array calculate_density calls -- which are not callable from a device region -- with an element-wise path built on the device-callable dispatcher calculate_density_elem_loc (327f88470). Still host-executed here; directives come in C. - MOM_EOS: get_EOS_form_and_scaling gains an optional R_to_kg_m3 out-arg (backward-compatible; kappa_shear's two callers are unaffected). Needed to reproduce calculate_density_1d's unity fast-path test and non-unity rescale exactly. - int_density_dz_generic_plm: resolve eos_form + the five unit-scaling factors once on the host (the accessor and MOM_error are not device-callable), and precompute eos_unity mirroring calculate_density_1d's no-rescale test. Under #ifdef __NVCOMPILER_OPENMP_GPU, replace the three non-Stanley calculate_density array calls (phase 1 vertical r5; phase 2 x-integral r15; phase 3 y-integral r15) with element loops calling calculate_density_elem_loc, reproducing calculate_density_1d(..., EOSdom_*[, rho_ref]) bit-for-bit: same _loc kernel per point, same rescale math (Ta=C_to_degC*T, ..., rho*=kg_m3_to_R, rho_ref->R_to_kg_m3*rho_ref). The eos_unity branch is hoisted outside the point loop so the loop bodies are branch-free (offload-ready; trap 3). Non-GPU builds keep the original array call verbatim in #else. The Stanley calls, and the int_density_dz_generic_ppm / _pcm siblings, are untouched. - On GPU builds, FATAL host-side (before any device region) when the non-Stanley path hits an unsupported EOS form / anomaly mode: the dispatcher covers only ROQUET_RHO (density and anomaly) and buggy_Wright (density only). Bitwise-safety: (i,n)/(i,n15) loops cover exactly the EOSdom_h5 / EOSdom_q15 / EOSdom_h15 index ranges the array calls work on; per-point density values are order-independent, so the loop order does not matter; the _loc kernels are the identical functions the polymorphic path already wraps. Multiplying by an unscaled factor of exactly 1.0 is bit-exact, so the non-unity branch matches calculate_density_1d's guarded rho_scale multiply. Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical to the dev/gpu reference (np=1); the #ifdef GPU element path is the executed arithmetic on the GPU build. DEBUG=True redundancy gate clean (0 'redundant'). "(Ocean pressure force)" clock ~20.4 s/24 calls, unchanged (host execution; directives land in C). Not runtime-exercised by the benchmark (inspection-only, code-identical to the array path): the Stanley branch, the use_rho_ref=.false. (inaccurate-form) branch, and the non-unity unit-scaling branch. Co-Authored-By: Claude Opus <noreply@anthropic.com>
First quadrature-phase offload of the finite-volume PGF density integrals. The phase-1 vertical-integral loop (-> dpa, intz_dpa) now runs in a !$omp target teams loop collapse(2) over (j,i) when the configuration is device-supported; the host loop is unchanged for every other case. - New host-resolved gate offload_phase1 = (.not.use_stanley_eos) .and. eos_unity .and. use_rho_ref, set only on GPU builds (#ifdef). When false (CPU build, Stanley, non-unity units, or inaccurate-form) the original host loop runs verbatim -- no code duplicated (single else block). - Device kernel: per-(i,j) size-5 quadrature locals T5l/S5l/p5l/r5l (compile-time sized, so no device auto-allocation) reproduce the host point-for-point, call the declare-target dispatcher calculate_density_elem_loc, and write dpa(i,j)/intz_dpa(i,j) via Boole's rule. dpa/intz_dpa are independent across (i,j). - Residency: e, dpa, intz_dpa are already device-resident (caller maps); the region carries no map() clauses and relies on presence. z0pres (a host-computed local) is mapped to/released around the region. The layer edge reconstructions T_t/T_b/S_t/S_b are mapped to() once per call by the caller (PressureForce_FV_Bouss) before the k-loop and released after (guarded to the plm/ppm path); they are read-only and identical for all k. - After the region, target update from(dpa) (and intz_dpa when present) refreshes the host: phases 2 & 3 (still host) read dpa, and the caller's kept "update to(dpa,...)" would otherwise overwrite the device values with a stale host copy. This intermediate copy-back is removed once phases 2/3 also run on device (increment D) and the caller round-trip is deleted (E). Bitwise-safety: the size-5 locals hold the identical per-point values as the old i*5+n slices; rho_anom/dpa/intz_dpa are computed independently per (i,j) from those points, so collapsing the (j,i) loops reorders no FP op; the _loc kernel and Boole coefficients are unchanged. The device path is entered only for the unscaled, accurate, non-Stanley configuration; all others keep the verbatim host arithmetic. Sole caller: PressureForce_FV_Bouss (audited: only call site of int_density_dz_generic_plm), where all no-clause arrays are device-present, so the compiler's implicit map(tofrom:) are runtime no-ops. Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/incC_verify.log. Execution proof: NVCOMPILER_ACC_NOTIFY=1 shows the F1L639 kernel launching (tmp_local_artifacts/notify_incC.log). Compile proof: -Minfo reports "Generating NVIDIA GPU code / Loop parallelized across teams collapse(2)" at the teams loop. DEBUG=True redundancy gate clean, 0 'redundant' (tmp_local_artifacts/debug_incC.log). "(Ocean pressure force)" clock 19.49 s/24 (was 20.40), consistent with phase 1 on device; the ~2-3 s target awaits D+E (round-trip removal). Not benchmark-exercised (inspection-only): Stanley, inaccurate-form, and non-unity-scaling paths (host, unchanged). Co-Authored-By: Claude Opus <noreply@anthropic.com>
…grals) to GPU Second and third quadrature-phase offload of the finite-volume PGF density integrals. The x-integral (-> intx_dpa) and y-integral (-> inty_dpa) loops now each run in a !$omp target teams loop collapse(2) when the configuration is device-supported; the host loops are unchanged for every other case. - New host-resolved gates offload_phase2/offload_phase3 = offload_phase1's condition, kept as separate flags so phase 3 can be forced off for x-only failure bisection. - Device kernels: per-cell size-15 q-point locals T15l/S15l/p15l/r15l and dz_x_l/dz_y_l (compile-time sized) reproduce the host loops point-for-point, call the declare-target dispatcher calculate_density_elem_loc for all 15 points, and write intx_dpa(I,j)/inty_dpa(i,J) via Boole's rule. Both the hWght>0 mass-weighting branch and its else branch are carried verbatim from the host path (they compute the same set of private corner scalars, so there is no divergent conditional array write). intx_dpa/inty_dpa are independent across cells. - Residency: dpa (written on the device by phase 1) is now read by phases 2/3 on the device with no host round-trip. z0pres is mapped once to bracket all three phases (its enter-data moved above phase 1, release moved after phase 3). e/dpa/intx_dpa/inty_dpa are already device-resident; T_t/T_b/S_t/S_b are mapped by the caller (unchanged from the phase-1 increment). - Phase 1's interim "update from(dpa,intz_dpa)" is replaced by a single consolidated "update from(dpa, intz_dpa, intx_dpa, inty_dpa)" after phase 3, so the caller's kept "update to(dpa,...)" still sees current host values. That consolidated copy-back and the caller round-trip are removed in increment E. Bitwise-safety: the size-15 locals hold the identical per-point values as the old i*15+(m-2)*5 slices; the pressure recurrence over n and the Boole quadrature are copied verbatim and run serially within each cell's thread; intx_dpa/inty_dpa are computed independently per cell, so collapsing the (j,I)/(J,i) loops reorders no FP op. The device path is entered only for the unscaled, accurate, non-Stanley configuration; all others keep the verbatim host arithmetic (host else-paths are byte-identical to the prior increment, verified by diff). Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/incD_verify.log. Compile proof: -Minfo generates GPU kernels F1L750 (x) and F1L946 (y), each "Loop parallelized across teams collapse(2)". DEBUG=True redundancy gate clean, 0 'redundant' (tmp_local_artifacts/debug_incD.log). "(Ocean pressure force)" clock 9.66 s/24 (was 19.49), halved by moving the horizontal integrals on-device; the remaining host round-trip is removed in increment E. Not benchmark-exercised (inspection-only): Stanley, inaccurate-form, and non-unity-scaling host paths. Co-Authored-By: Claude Opus <noreply@anthropic.com>
…increment E, caller side)
Removes the caller half of the int_density_dz_generic_plm device round-trip:
the four integral outputs (dpa, intz_dpa, intx_dpa, inty_dpa) are now consumed
entirely in the caller's own device regions (interface-pressure accumulation and
the PFu/PFv loops), so the host no longer pushes them to the device before use.
- New host-resolved gate idz_on_device, set only on GPU builds (#ifdef
__NVCOMPILER_OPENMP_GPU), that mirrors int_density_dz_generic_plm's internal
offload_phase1 condition EXACTLY on the generic_plm call path: use_EOS .and.
use_ALE .and. Recon in {1,3} .and. .not.use_stanley_pgf .and.
.not.use_inaccurate_pgf_rho_anom .and. the four EOS unit-rescale factors == 1
(queried via get_EOS_form_and_scaling). When true, the routine leaves its
outputs device-resident; when false (CPU build, Recon scheme 2, Stanley,
inaccurate-form, non-unity units) the host round-trip is kept verbatim.
- Drop the per-call "defensive" update to(dpa,intx_dpa,inty_dpa,intz_dpa) when
idz_on_device: the generic_plm kernels already wrote them on the device, so the
copy would clobber them with a stale host copy. Kept (guarded) for the host path.
- Narrow the full update from(e): e is built on the device, and under
idz_on_device the blanket ALE pull is skipped. The pull is retained for its
genuine host consumers -- host-run int_density, the correction/reset-intxpa
reference-interface searches, the MassWt diagnostics, and (this is the
stale-host fix) the SAL/tides "new answers" blocks that read e(:,:,1) and
read-modify e(:,:,K) on the host.
- Add a guarded update from(intx_dpa,inty_dpa) at routine exit for the SAL/tides
horizontal-gradient diagnostics (the only host readers of those two arrays),
guarded to idz_on_device .and. .not.bq_sal_tides .and. (calculate_SAL .or. tides).
Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical
to the dev/gpu reference (np=1); log tmp_local_artifacts/incE_fix_verify.log. The
SAL/tides and correction/reset guard terms are false in benchmark_ALE, so this is
inspection-only for those configs; "(Ocean pressure force)" clock ~2.3 s/24 (was
~20). Note the SAL/tides-new-answers GPU write-back (host e/pa mods pushed back to
the device) is a separate, pre-existing gap -- idz_on_device is the fully-GPU,
non-SAL/tides PGF path. Remaining gates: np>=2 multi-GPU verify (1 V100 here) and a
plain-CPU compile of the #else branches.
Co-Authored-By: Claude Opus <noreply@anthropic.com>
…ncrement E) Removes the routine half of the device round-trip. Phases 1-3 already compute dpa/intz_dpa/intx_dpa/inty_dpa on the device; this drops the consolidated device->host copy-back at the end of the routine. The sole caller (PressureForce_FV_Bouss) now consumes all four outputs in its own device regions (prior commit), so no host refresh is needed. Only the per-call z0pres scratch release is kept. Also rewrites the offload_phase2/offload_phase3 comment: they must track offload_phase1 now that the round-trip is gone. Decoupling a phase would leave its host-computed integral stranded on the host while the caller (which mirrors this gate as idz_on_device and no longer copies back) reads a stale device copy. The "settable independently for x-only failure bisection" affordance no longer holds. Bitwise-safety: no arithmetic changes; only OpenMP data-movement directives are removed. On CPU builds (offload_phase1/2/3 all false) the removed blocks were already skipped, so host builds are unaffected. Verification: benchmark_ALE (ROQUET_RHO, Boussinesq, PLM) ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/incE_fix_verify.log. Co-Authored-By: Claude Opus <noreply@anthropic.com>
… the GPU Mirrors the Calc_kappa_shear_vertex offload for the tracer-point driver used when VERTEX_SHEAR=False (the plain benchmark case). All the shared machinery was already in place -- the three column workers are declare target, the density-derivs element dispatcher covers WRIGHT and ROQUET_RHO, eos_form/scaling are hoisted host-side in this driver too, and GPU_nk_max sizes the per-column scratch -- so this is only the final loop-flip. - Hoist thickness_to_dz out of the j-loop to a whole-domain dz_3d (halo_size=0), and precompute surface_pres_2d from the (possibly-null) p_surf pointer on the host, ahead of the device region. - Flip the column loop from !$OMP parallel do (host) to !$omp target teams loop collapse(2) over (j,i). The per-J 2-D slabs (h_2d/u_2d/...) are gone: the device loop reads the mapped 3-D inputs (h, u_in, v_in, tv%T, tv%S, dz_3d) directly -- the h-point driver has no vertex interpolation, so no slab staging is needed. Results stage in device-only kappa_3d/tke_3d (was the per-J kappa_2d/tke_2d); a masked do concurrent then writes kappa_io/tke_io/kv_io. - Per-column private scratch (Idz/h_lay/.../kappa/tke/.../kc/kf) is GPU_nk_max-sized under #ifdef __NVCOMPILER_OPENMP_GPU (compile-time constant -> device stack, not the serialized device-heap path); SZK_(GV) on CPU builds. Bitwise-safe (no size()/bounds intrinsics on these in the routine). - Residency: h host-authoritative -> map(to:)+update to; u_in/v_in (set_diffusivity's u_h/v_h) and dz_3d are fresh host locals -> map(to:); CS shallow map(to:); outputs map(to:)+update to (preserve host halos) then update from; diag_* mapped after their conditional host zeroing, copy-back guarded by the same predicate. tv%T/tv%S are read directly on the device here (this path passes raw tv, no full_convection filtering, unlike the vertex driver's fresh T_f/S_f): the T/S data is persistently mapped via MOM's CS%tv but on a different descriptor, so this routine's tv dummy is explicitly mapped map(to: tv, tv%T, tv%S) + update to (the MOM_set_viscosity idiom) so the device can resolve tv%T/tv%S -- update to alone does not attach the descriptor. Verification: benchmark (VERTEX_SHEAR=False, EQN_OF_STATE=WRIGHT, Boussinesq) ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/kshpoint_verify2.log. Execution proof: NVCOMPILER_ACC_NOTIFY shows the device kernel calculate_kappa_shear line 328 (target teams loop) launching grid=64800 (360x180 columns) once per call, plus the write-back do concurrent at line 469. "(Ocean kappa_shear)" ~0.78 s/6 hits (the benchmark kappa_shear workload is small; the kernel is teams-only block=(1,1,1), the same occupancy headroom the vertex driver carries). Untested-but-code-identical: the non-unity unit-rescale branch and the non-temperature (rho) branch. Deferred pre-merge gate: >=2-GPU / np>=2 verify (this node has 1 V100). Co-Authored-By: Claude Opus <noreply@anthropic.com> Co-authored-by: Utheri Wagura <Utheri.Wagura@noaa.gov>
…ix W-0155) The two !$omp target teams loop regions in tracer_epipycnal_ML_diff had per-column private arrays declared dimension(SZK_(GV)) / dimension(SZK_(GV)*2), i.e. sized by the runtime GV%ke. NVFORTRAN cannot place a runtime-sized private array in local (stack) memory (NVFORTRAN-W-0155) and spills it to slow device global memory, with a "results may be incorrect" warning. Fix mirrors the kappa_shear GPU_nk_max pattern (KNOWLEDGE row 21): under #ifdef __NVCOMPILER_OPENMP_GPU the per-column scratch (h_demand_L/R, h_used_L/R, h_supply_frac_L/R real, and kbs_Lp/Rp integer + left_set/right_set logical, the last four dimension GPU_nk_max*2) is given the compile-time-constant bound GPU_nk_max=128; CPU builds keep the exact SZK_(GV) sizes. A GPU-only FATAL guard in tracer_epipycnal_ML_diff catches GV%ke > GPU_nk_max instead of silently corrupting. Bitwise-safe: no size()/ubound()/lbound() on any of these arrays (loops are bound by num_srt(i,j) <= nz), so enlarging the declared extent changes no arithmetic. Pre-existing issue (not from the pressure-force/kappa stack); surfaced by a clean rebuild. Verification: benchmark_ALE (which exercises "Ocean epipycnal diffuse tracer") ocean.stats bit-identical to the dev/gpu reference (np=1); the ten W-0155 lines are gone from the build log. Deferred: >=2-GPU verify (1 V100 here). Co-Authored-By: Claude Opus <noreply@anthropic.com>
…the ZSTAR regrid
Verbatim refactor (no OpenMP directives, no behaviour change) that prepares the ALE
ZSTAR regrid for GPU offload. The per-column regrid kernels dereference the regridding
control structure on entry, and regridding_CS holds a nested zlike_CS (a pointer) whose
coordinateResolution is allocatable -- exactly the nested-derived-type-with-allocatable
shape that is fragile to map onto the device. Rather than map it, the scalar/array
parameters are hoisted so the kernels take them explicitly.
- coord_zlike: new build_zstar_column_loc(nk, min_thick, coordinateResolution, ...) holds
the former build_zstar_column body verbatim; build_zstar_column becomes a thin wrapper
passing CS%nk/CS%min_thickness/CS%coordinateResolution. _loc is public.
- MOM_regridding: new filtered_grid_motion_loc(nk_tgt, filt_shallow, filt_deep,
old_grid_wt, ...) and adjust_interface_motion_loc(nk_tgt, min_thick, ...) hold the
former bodies verbatim (CS%nk -> nk_tgt, CS%depth_of_time_filter_* -> filt_*,
CS%old_grid_weight -> old_grid_wt, CS%min_thickness -> min_thick); the CS-taking
filtered_grid_motion / adjust_interface_motion become thin wrappers.
No call site changes anywhere: every existing caller (build_{zstar,sigma,rho,...}_grid,
MOM_diag_remap) keeps using the CS-taking wrappers unchanged. This is the first of two
commits; the device directives + residency come next (GPU-gated).
Bitwise-safety: the _loc bodies are the originals with CS%x textually replaced by the
matching argument; the wrappers pass exactly those CS fields, so the arithmetic and the
values are identical. Verification: benchmark_ALE ocean.stats bit-identical to the
dev/gpu reference (np=1) -- CPU-path gate (the _loc kernels still run on the host; no
directives yet). log tmp_local_artifacts/ale_refactor_verify3.log.
Co-Authored-By: Claude Opus <noreply@anthropic.com>
Second of two commits (the first extracted the device-callable _loc kernels). This offloads build_zstar_grid's per-column loop to the device for the ZSTAR coordinate (the mode benchmark_ALE uses), building on the _loc refactor. - build_zstar_column_loc / filtered_grid_motion_loc / adjust_interface_motion_loc are now !$omp declare target. The "should never happen" grid sign/tangling/negative-thickness MOM_error + write + stop asserts, and the debug-only blocks (debug=.false.) and their z_act scratch, are #ifndef __NVCOMPILER_OPENMP_GPU (kept on the host, compiled out of device builds). filtered_grid_motion_loc's massless-column early return is expressed as an if-guard (if (.not. massless)) since device loops cannot return early. - build_zstar_column_loc's zScale is now a required positional arg (z_scale) ahead of the two optional shelf args, so the device call sites need no keyword arguments (nvfortran rejects keyword-form args inside a target region). - build_zstar_grid: the !$OMP parallel do over (j,i) becomes a !$omp target teams loop collapse(2); the land-mask `cycle` is an if-guard; the regridding_CS scalars + coordinateResolution are hoisted into plain locals on the host (the private component is fetched via the new get_zlike_coord_res accessor) so the device loop never dereferences the nested (pointer + allocatable) control structure; per-column zNew/zOld are fixed-size (GPU_nk_max=128, module parameter) under #ifdef with a FATAL guard when CS%nk/GV%ke exceed it. - Residency: h is host-authoritative in the (host) ALE stack -> update to(h); nom_depth_H and the hoisted coordRes -> map(to:); dzInterface is device-written and copied back (update from) for the subsequent host remap that reads it (ALE_remap_set_h_vel_via_dz, remap_OBC). calc_h_new_by_dz and the rest of ALE stay on the host for now. Verification: benchmark_ALE (ZSTAR, PLM, ROQUET_RHO, Boussinesq) ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/ale_port_verify3.log. Execution proof: NVCOMPILER_ACC_NOTIFY shows the device kernel build_zstar_grid line 1728 launching grid=2432 per call (teams-only block=(1,1,1); the _loc kernels are inlined). Untested-but-code-identical: the ice-shelf (frac_shelf_h) branch and the RHO/SIGMA/HYCOM coordinate builders (host, unchanged). Deferred pre-merge gate: >=2-GPU / np>=2 verify (this node has 1 V100). Co-Authored-By: Claude Opus <noreply@anthropic.com>
…t/port-set_diffusivity
First device increment of the set_diffusivity port, on Ed's j-blocked groundwork. find_N2 is the dominant, self-contained, EOS-bearing piece of the active benchmark_ALE path (find_TKE_to_Kd and set_density_ratios are never called there and stay host). - EOS: the per-interface calculate_density_derivs (ROQUET_RHO) is evaluated on the device through the declare-target calculate_density_derivs_elem_loc dispatcher. eos_form + unit scaling are resolved once on the host (get_EOS_form_and_scaling, with a GPU-only FATAL if the form is not ROQUET_RHO/WRIGHT); the scaling is always applied, so the unscaled (unity) case reproduces the former array call bit-for-bit. This is the kappa_shear increment-2b pattern. - Loop forms: the pressure recurrence (serial down/up-K) + the two bottom-boundary-layer do k=nz,2,-1 loops become per-column !$omp target teams loop collapse(2) over (j,i) with a serial k inside; the former cross-column do_any/exit early-out is dropped -- each device thread runs its own k loop and terminates via its do_i flag, bit-for-bit equivalent. The elementwise N2_lay/N2_int, dRho_int init/copy and dz_BBL_avg loops become do concurrent. - Residency: the full-domain EOS inputs (T_f/S_f and the tv descriptor) are mapped ONCE by the caller (set_diffusivity) around the whole jstart loop, and tv%T/tv%S are refreshed there (update to) since they are the diabatic-mutated persistently-mapped fields; h is already persistently device-resident. find_N2 itself maps only the per-block scratch (dz/pres/h_amp + the 3-D block arrays); doing the full-domain maps per block would re-copy them on every j-row (njblock=1) and dominate the runtime. dRho_int/N2_int/N2_lay/N2_bot (find_N2 outputs) and pres/dz_BBL_avg (read by the host find_rho_bottom that follows) are update from before the per-block region closes. find_rho_bottom stays on the host. Verification: benchmark_ALE ocean.stats bit-identical to the dev/gpu reference (np=1); logs tmp_local_artifacts/setdiff_findN2_verify.log (+ ..._perf_verify.log for the residency fix). Execution proof: NVCOMPILER_ACC_NOTIFY shows the device kernels find_n2 lines ~1284 (pres+EOS), ~1333 and ~1368 (bottom-BL) launching per block. Timing: "(Ocean set_diffusivity)" 3.4 s/24 (host baseline 3.2 s; find_N2-on-GPU is ~break-even at this size because the bkgnd/assembly/add_LOTW_BBL tail is still host, so find_N2 pays a host<->device boundary each call and the njblock=1 kernels are teams-only/low-occupancy; the win materialises once the surrounding path is on-device). Remaining set_diffusivity tail left on host by design. Deferred: >=2-GPU verify (1 V100 here). Co-Authored-By: Claude Opus <noreply@anthropic.com>
…le (increment 1) First increment of the CPU-preserving ePBL GPU port (the ePBL half of the mislabeled "(Ocean vert remap)" region). Bitwise-neutral enablers: !$omp declare target the pure per-column kernels so ePBL_column can call them once it (and the driver loop) are offloaded in later increments. The same source serves the host and device paths. - MOM_energetic_PBL.F90: module-scope declare target on the pure scalar helpers exp_decay_TKE_adjust, find_PE_chg, find_PE_chg_orig, find_Kd_from_PE_chg, find_mstar, mstar_Langmuir (no host-sinks, no runtime-sized array locals -> nothing to gate). - MOM_wave_interface.F90: declare target on get_Langmuir_Number, get_StokesSL_LiFoxKemper, ust_2_u10_coare3p5 (the LF17 statistical-wave Langmuir path that benchmark_ALE uses: USE_LA_LI2016=True, USE_WAVES=False). get_Langmuir_Number is #ifdef-split: on GPU builds only the LF17 branch (call get_StokesSL_LiFoxKemper) is compiled -- the wave-model branches (TESTPROF/ SURFBANDS/DHH85, their US_H/VS_H(SZK) + allocatable StkBand scratch, the Get_SL_Average_* calls, the Null-method + misalignment MOM_errors) are host-only. This reproduces the host LF17 result bit-for-bit (Dpt_LASL is unused on the LF17 path; LA_Misalignment defaults false). ust_2_u10_coare3p5's vonKar<0 config MOM_error is #ifndef-gated off device. TODO before merge: a GPU build currently compiles the LF17-only get_Langmuir_Number for all callers; add a host-side FATAL guard for GPU builds with a non-LF17 WaveMethod or LA_Misalignment (benign for benchmark_ALE, which is LF17/no-misalignment). Verification: benchmark_ALE ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/epbl_inc1_verify.log. This device-compiles the full ePBL+wave kernel set (the key de-risking) while being bitwise-neutral (device versions unused until ePBL_column + the driver loop are offloaded). Remaining: increment 2 = declare target ePBL_column + fixed-size its SZK locals; increment 3 = flip the energetic_PBL driver loop to a target region + residency. Deferred: >=2-GPU verify. Co-Authored-By: Claude Opus <noreply@anthropic.com>
…ment 2) Second increment of the CPU-preserving ePBL GPU port. ePBL_column (the ~1061-line iterative single-column solver: outer MLD iteration + down-K TKE budget with an inner Newton solve) is now !$omp declare target, so it device-compiles with its declare-target helper kernels (find_mstar, find_PE_chg*, get_Langmuir_Number, ...) from increment 1. Still bitwise-neutral: ePBL_column is host-called until the driver loop is offloaded (increment 3); on GPU builds its host version uses the #ifdef paths below, which are bit-identical for benchmark_ALE. - GPU_nk_max=128 module parameter; ePBL_column's per-column private automatic arrays (pres_Z/hb_hs, the dT_to_dColHt.../hp_a.../Sh_b block, MixLen_shape/h_dz_int/Kddt_h, and mech_TKE_k/conv_PErel_k/ nstar_k/dT_expect/dS_expect/num_itts) are given compile-time-constant sizes under #ifdef __NVCOMPILER_OPENMP_GPU (KNOWLEDGE row 21), keeping SZK_(GV) on CPU builds. The dummy args and the dimension(20) iteration arrays are already fine. (ePBL_BBL_column's analogous locals are fixed-sized too via the shared decls; it is dead in benchmark_ALE, not declare target.) - The equation-discovery (eqdisc) mixing-length/velocity calls (kappa_eqdisc/get_eqdisc_v0/v0h) are #ifndef-gated off device (dead in benchmark_ALE: EPBL_EQD_* absent -> .false.); the transLay mixing-length + v0_scale=0 paths run, bit-identical. - FATAL guards in energetic_PBL_init on GPU builds: GV%ke > GPU_nk_max, and any eqdisc* enabled. TODO before merge (carried from increment 1): host-side guard for GPU builds with a non-LF17 WaveMethod / LA_Misalignment (get_Langmuir_Number compiles the LF17-only body on GPU). Verification: benchmark_ALE ocean.stats bit-identical to the dev/gpu reference (np=1); log tmp_local_artifacts/epbl_inc2_verify.log. ePBL_column device-compiles cleanly (the key de-risking). Remaining: increment 3 = flip the energetic_PBL driver loop to a target teams loop + one driver-level data region (residency) so the solver runs on the device. Deferred: >=2-GPU verify. Co-Authored-By: Claude Opus <noreply@anthropic.com>
Flip the energetic_PBL driver (do j / do i over ocean columns calling the now-declare-target ePBL_column) to a single `!$omp target teams loop collapse(2)` with one driver-level data region. Columns are read directly from the 3-d inputs (2-d staging slabs and the per-j Kd_2d staging removed, bit-identical); thickness_to_dz is hoisted whole-domain (dz_3d); the ~14 per-column scratch arrays are teams-loop private + fixed-size under GPU_nk_max. Residency: map(to:) the inputs + CS + tv/tv%T/tv%S (with update-to for the map-to-present trap) + Waves; map(from:) Kd_int / ML_depth / BBL_depth. Forcing (fluxes%ustar) and the MLD first guess are hoisted to plain host arrays so the device never touches a fluxes% pointer member or a CS allocatable. Dead-in-benchmark_ALE branches (BBL mixing, options_diff, stochastics, per-column diagnostics, report_avg_its) are #ifndef'd off the device path and backed by runtime FATAL guards so an unsupported config fails loudly rather than silently wrong. CPU path is preserved verbatim under #else. Device-callable Langmuir chain (MOM_wave_interface): the LF17 path dereferenced the Waves *pointer* on-device, which nvfortran will not attach under -gpu=mem:separate. Replace it with a pointer-free `wave_LF17_params` bundle (allocatable-free, maps cleanly) populated host-side by set_wave_LF17_params and threaded through get_Langmuir_Number_LF17 -> get_StokesSL_LiFoxKemper -> ust_2_u10_coare3p5. get_Langmuir_Number_LF17 uses the in-routine `!$omp declare target` block form (the module name-list form silently failed to emit its device symbol). get_Langmuir_Number is restored to its full CPU form (no longer declare-target). STATUS: ePBL executes correctly on the device (no illegal-address; T/S/mass bit-identical to displayed precision). NOT YET bitwise: the energy diagnostic differs at ~1e-13, the documented host-vs-device x**(1/n) transcendental floor (KNOWLEDGE row 5) from the 4 `(...)**C1_3` cube-roots in the ePBL_column vstar budget and the Langmuir `LA**(-1.33)` power. Closing it needs cuberoot() for the **C1_3 sites + an explicit-form rewrite of the arbitrary power, applied to the one CPU+GPU source, plus a regenerated golden (the fix changes the CPU answer). Deferred by decision. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ispatcher (increment 1) Offload prep for the thickness_diffuse_full GPU port. Swap the three live calculate_density_derivs generic array calls (the u-point slopes :933, the v-point slopes :1246, and the surface u block :1545) for the element-wise device-callable dispatcher calculate_density_derivs_elem_loc, folded into the existing fill loops so the evaluated range matches the array call exactly. get_EOS_form_and_scaling resolves the form + unit factors once (host); the scaling is always applied (unity-exact for the unscaled case), reproducing the former array calls bit-for-bit. Adds the same GPU-build FATAL guard as find_N2 for EOS forms without a device kernel. This mirrors the find_N2/kappa EOS-dispatcher pattern; the calculate_density _second_derivs (Stanley) path is left untouched (dead in benchmark_ALE). No device directives execute yet (elem_loc runs on the host here) -- this is prep so the u/v slope loops can be offloaded next without a polymorphic EOS call on the device. The GPU port of this file follows the structure of Ed Yang's edoyango/port/thickness_diffuse branch (the (i,j,k)-blocking and directive placement) as a reference; his branch is not merged (it is stale and overlaps our density/EOS work) but informs the approach. BITWISE: neutral -- benchmark_ALE np=1 ocean.stats is byte-identical to the pre-increment (ePBL-device) baseline (tmp_local_artifacts/td_baseline_ocean .stats). The only residual vs the canonical golden repro_ocean.stats is the separately-deferred ePBL x**(1/n) transcendental floor (~1e-13 in energy), which this commit does not touch. Co-authored-by: Edward Yang <edward.yang@anu.edu.au>
…d prep) Promote the ten per-column slope-scratch arrays of thickness_diffuse_full from 2-D (i,K) to 3-D (i,j,K) so the u/v slope+flux loops can later be fully collapsed onto the device with the scratch device-resident across the slope->flux passes: drdi_u/drdj_v, drdkDe_u/_v, dzN2_u/_v, Sfn_unlim_u/_v, slope2_Ratio_u/_v (v-arrays carry the SZJB_ symmetric-memory j-extent). Declarations + every index site updated (~50 sites across the u-block and v-block). c2_dz_u/_v are left 2-D (FGNV path, dead in benchmark_ALE, host-only). No directives yet: this is a pure storage transform, computation unchanged. The promoted arrays remain in the CPU !$OMP parallel do private() clauses transitionally (per-thread copies, still bit-for-bit; the do-j loop writes disjoint j-slices) -- the target-teams-loop offload (increment 2b) replaces those clauses. BITWISE: neutral -- benchmark_ALE np=1 ocean.stats byte-identical to the pre-increment ePBL-device baseline (tmp_local_artifacts/td_baseline_ocean.stats, diff empty). Only residual vs the canonical golden is the separately-deferred ePBL transcendental floor. The GPU offload of this file follows the structure of Ed Yang's edoyango/port/thickness_diffuse (3-D scratch promotion + full collapse) as a reference. Co-authored-by: Edward Yang <edward.yang@anu.edu.au>
Offload the whole thickness_diffuse_full compute (setup + u/v slope + streamfunction + flux + MEKE work) to the device, one target-data region, building on the 3-D scratch of increment 2a and the elem_loc EOS of increment 1. - CPU-neutral refactor first (verified bitwise-neutral before any directive): fuse the per-K EOS-fill do-I loop into the slope do-I loop (scalarizing the per-I EOS temporaries) so the slope nest is a tight do j / do K / do I; split the outer do-j into slope and flux j-loops; swap the flux loops to do I / serial do K (per-column uhtot/vhtot recurrence order preserved); move the now-3-D scratch out of the OMP private() clauses to shared. - Loop forms (all `#ifdef __NVCOMPILER_OPENMP_GPU`, original !$OMP parallel do kept under #else so CPU threading is retained): slope loops -> `target teams loop collapse(3)` over (j,K,I) with per-iteration scalars private; the uhtot/vhtot flux recurrence + the h_avail_rsum setup running sum -> collapse(2) + serial do K (KNOWLEDGE 4.1 branch 5); the remaining elementwise/init/assembly loops -> collapse(2). 16 device kernels generate, no W-0155. - Dead-in-benchmark_ALE paths #ifndef'd off the device with host FATAL guards: use_stanley (second derivs), FGNV streamfn + streamfn_solver + c2_dz, skeb_use_gm/STOCH, non-Boussinesq tv%SpV_avg, and the slope/sfn diagnostics. The layer-1 v-point EOS array call is converted to the elem_loc dispatcher. - Residency: one target enter/exit data region; map(to:) h (+update to), e, dz, T, S, Kh_u/v, int_slope_u/v, CS, slope_x/y (present-guarded) and the host-zeroed Slope_/hN2 PE arrays; map(alloc:) pres, dz-scratch, h_avail*, uhtot/vhtot, Work_u/v, uhD/vhD and the 3-D slope scratch; map(from:) uhD/vhD. All exit maps use `release` (never `delete`). - Pointer/allocatable hoist (the map-to-present / Waves-pointer-trap class): MEKE%GM_src, CS%GMwork and tv%p_surf hoisted to plain local arrays (copied in, written/read on device, copied back) so no allocatable component of a pointer/derived type is dereferenced on the device. Two fixes found via the crash -> compute-sanitizer -> fix loop: (1) tv%p_surf IS associated in benchmark_ALE, so hoist it rather than FATAL; (2) an illegal address that surfaced downstream in set_viscous_BBL (pbv attach) was traced (sanitizer + a revert-to-2a control build) to this region -- fixed by making the whole region pure OpenMP target (no `do concurrent` inside the target-data region) and by using map(release:) instead of map(delete:) on exit (KNOWLEDGE 8: a per-call delete can zero an outer/persistent mapping refcount). Note tv/tv%T/tv%S are deliberately NOT mapped here: on device the density inputs are the host-filled local arrays T/S from vert_fill_TS (mapped to: fresh), so tv%T/tv%S are never read on device. BITWISE: neutral -- benchmark_ALE np=1 ocean.stats byte-identical to the pre-increment ePBL-device baseline (tmp_local_artifacts/td_baseline_ocean.stats, diff empty), verified on a clean from-scratch rebuild. Device execution confirmed (16 kernels; full 24-step run, no CUDA fault; compute-sanitizer memcheck clean). Only residual vs the canonical golden is the separately- deferred ePBL x**(1/n) transcendental floor. Structure follows Ed Yang's edoyango/port/thickness_diffuse (3-D scratch + full collapse); our EOS goes through the elem_loc dispatcher rather than his EOS-3d stack. Co-authored-by: Edward Yang <edward.yang@anu.edu.au>
…(TODO for Marshall) Comment-only. ePBL is the one offloaded routine on this branch that is not yet GPU/CPU bit-for-bit: the device run differs from host at ~1e-13 in the energy diagnostic (T/S/mass bit-identical) solely from device-vs-host libm/libdevice divergence in transcendental ops in its on-device kernels, which have no deterministic drop-in (unlike the x**(1/3) cube roots, already handled via cuberoot() on the answer_date>=20240101 path). Add a consolidated, greppable `TODO(gpu-bitwise-repro)` audit at the module level enumerating the live (benchmark_ALE) transcendental sites for the real matched-transcendental fix (being done by Marshall -- not duplicated here): find_PE_chg exp() (core, likely dominant), exp_decay_TKE_adjust exp(), find_mstar mstar_N log() (OM4 Ekman), mstar_Langmuir **(-1.33) (arbitrary power), and the MixLen pow(x,1.0) (trivial identity stopgap noted but not applied). MKE_src exp() is inert (MKE_TO_TKE_EFFIC=0). No numerical change, so the ePBL-device baseline and the no-new-diff offload gate are unaffected; an nvfortran-CPU reference build is staged at ocean_only/cpu_build for the GPU==CPU one-source verification once Marshall's fix lands. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… (TODO for Marshall) Comment-only. The LF17 Langmuir chain (get_StokesSL_LiFoxKemper / one_minus_exp_x / ust_2_u10_coare3p5) runs on-device when called from ePBL_column and contributes device-vs-host transcendental divergence to the ~1e-13 ePBL energy floor via LA -> mstar_LT. Add a greppable TODO(gpu-bitwise-repro) enumerating the live sites (exp in one_minus_exp_x; erfc in get_StokesSL r5 large-root branch; atan(1.0) for PI, usually folded; sqrt is exact) for Marshall's matched-transcendental pass, cross-referencing the ePBL audit. No numerical change. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…llable exp() The intrinsic exp() lowers to host libm on the CPU and CUDA libdevice on the GPU, whose last-bit rounding differs, so any do concurrent / omp target kernel that calls exp() is not bit-for-bit CPU==GPU -- the dominant source of the ~1e-13 ePBL energy-budget divergence on feat/port-ePBL (see the TODO(gpu-bitwise-repro) audit in MOM_energetic_PBL.F90). exp_reprod avoids the library call: Cody-Waite range reduction (x = k*ln2 + r with a two-part ln2 so k*ln2_hi is near-exact), a degree-12 Horner polynomial for exp(r) whose reciprocal-factorial coefficients are compile-time constant-folded (identical host/device), and scale(.,k) for 2**k. Every op is + - * / (plus nint/scale, exact), all bit-identical host vs device under -Mnofma -- so the result is reproducible by construction, the same way cuberoot/nth_root replace x**(1/n). declare target; compiles to device code (symbol confirmed); scale/nint verified device-callable. Accuracy: max relative error 3.1e-16 (~1.4 ULP) vs the intrinsic over x in [-70, 20]; unit test Test_exp_reprod added (tol 1e-14). Like cuberoot it is NOT bit-identical to intrinsic exp, so adopting it in a kernel changes answers and needs a golden regeneration. Nothing calls it yet (bitwise-neutral to commit); offered as the exp building block for the transcendental-repro pass. Log/erfc/arbitrary-pow still need the same treatment. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…llable log() Companion to exp_reprod: a device-callable natural log built from only +,-,*,/ (plus the exact exponent/fraction bit intrinsics), so it is bit-identical host vs device under -Mnofma. x = m*2**k via exponent/fraction; m reduced to [sqrt(1/2), sqrt(2)); log(m) = 2*(s + s^3/3 + ... ), s=(m-1)/(m+1), |s|<=0.172, as a degree-21 Horner poly in s^2 with compile-folded reciprocal-odd coeffs; then + k*ln2. Max relative error 3.7e-16 (~1.7 ULP) vs intrinsic log over x in [1e-30, 1e30]; unit test Test_log_reprod added (tol 1e-13). declare target; device code + symbol confirmed. Together exp_reprod + log_reprod make arbitrary real powers reproducible via x**y = exp_reprod(y*log_reprod(x)) -- covering the ePBL Langmuir Convect_Langmuir_Number**(-1.33) and the MixLen **1.0. Bitwise-neutral commit (nothing calls it yet). Building block for the transcendental-repro pass; erfc and site wiring + golden regen still to come. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…allable erfc() Third companion to exp_reprod/log_reprod, for x >= 0. Built from +,-,*,/ and exp_reprod only, so bit-identical host vs device under -Mnofma. Two ranges: x < 1.5 uses the erf Maclaurin series (48 fixed terms), erfc = 1 - erf; x >= 1.5 uses the incomplete-gamma continued fraction erfc = (x/sqrt(pi))*exp(-x^2)* Q-CF(1/2,x^2) via a FIXED-iteration (60) modified-Lentz recurrence (no early exit -> device-safe). Max relative error 5.8e-14 vs the intrinsic over x in (0.01,26] (worst near x~24 where erfc~1e-263 is physically zero; ~1e-15 across the physically relevant range). declare target; device code + symbol confirmed. Unit test Test_erfc_reprod added (tol 1e-12). This is the last transcendental in the ePBL/wave on-device path: get_StokesSL r5 calls erfc() (root_2kz>1e-3 branch), which was the sole remaining GPU-vs-CPU floor after exp/log/pow were wired to reprod kernels (~8.5e-20 residual). Bitwise-neutral commit (nothing calls it yet). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ls (GPU==CPU bitwise) Replace every device-executed intrinsic transcendental in the ePBL + LF17 Langmuir on-device path with the bit-reproducible kernels exp_reprod / log_reprod / erfc_reprod (committed just prior), so the whole path is built from +,-,*,/ (plus exact bit intrinsics) and is therefore bit-identical host vs device under -Mnofma. Only the transcendental symbol is swapped at each site; all other arithmetic is preserved character-for-character. MOM_energetic_PBL.F90: exp_decay_TKE_adjust + ePBL_column TKE-decay exp() -> exp_reprod; find_mstar mstar_N log() -> log_reprod; mstar_Langmuir Convect_Langmuir_Number**LT_enhance_exp -> exp_reprod(y*log_reprod(x)); MixLen **1.0 -> an ==1.0 identity fast-path. MOM_wave_interface.F90: one_minus_exp_x exp() -> exp_reprod; get_StokesSL r5 erfc() -> erfc_reprod; ust_2_u10_coare3p5 COARE-Cd log() -> log_reprod; PI = 4*atan(1) -> a literal (kills device atan). Dead branches (answer_date<2023/2024, RH18, BBL, Stanley) and the inert MKE_src exp (x MKE_TO_TKE_EFFIC=0) left as-is; sqrt is IEEE-exact host==device. RESULT: benchmark_ALE np=1 GPU build == a same-source nvfortran-CPU build (ocean_only/cpu_build, GPU-offload flags stripped) BIT-FOR-BIT (diff ocean.stats empty) -- ePBL is now fully GPU/CPU reproducible, closing the ~1e-13 energy-diagnostic floor that had made the branch non-bitwise since the ePBL offload. The reproducible kernels change answers slightly vs the intrinsics (like cuberoot vs x**(1/3)), so the canonical golden needs regenerating to this reference. These prototype kernels feed / may be superseded by Marshall's transcendental-repro pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…d kernels, GPU==CPU) Comment-only: flip the two TODO(gpu-bitwise-repro) audit blocks to RESOLVED now that the transcendentals are wired to exp_reprod/log_reprod/erfc_reprod and the GPU build reproduces a same-source nvfortran-CPU build bit-for-bit. Kept the site lists as a record; noted the prototype kernels may be superseded by Marshall's pass. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.