bodner but good by JorgeG94 · Pull Request #86 · MOM6-GPU/MOM6

JorgeG94 · 2026-07-28T23:57:50Z

This is a fresh implementation of the Bodner mixed layer restrat since I completely destroyed my previous branch when implementing blocking.

github-actions · 2026-07-29T00:15:39Z

GPU Port Coverage

Overall: 2261 / 5592 portable executed lines ported (40.4%)
Since base branch: +85 ported lines (+1.4 pp)

Files touched by this PR: 83 / 100 portable executed lines ported (83.0%)

_{Full per-file / per-routine breakdown: see the "gpu-port-report" job summary and artifact.}

Full line-by-line coverage report

JorgeG94 · 2026-08-10T01:14:21Z

With:

! Blank file in which we can put "overrides" for parameters
#override NIGLOBAL = 300
#override NJGLOBAL = 300
#override DT=800.0
#override DT_THERM=1600.0
#override DAYMAX=0.5

compiled with :

CC = mpicc
MPICC = mpicc
FC = mpifort
MPIFC = mpifort
PYTHON = python3

FCFLAGS = -O3 -g -qno-opt-dynamic-align -fp-model source -march=core-avx-i

even gives a bit of a speedup over CPU:

                             tmin                tmax               tavg
GPU ver: 60      0.234275      0.495331      0.301064 
CPU ver: 60      0.302732      0.475717      0.366187                  < dev/gfdl

tmax is a bit higher but the rest look very nice. Fine combing this PR a bit more...

edoyango

thanks for doing the cpu tests. the i-blocking still needs to be done. I've also added some other comments that should help get the cpu tmax down.

edoyango · 2026-08-13T23:26:23Z

-    keep_going = .true.
-    do k=1,nz
+
+    do jstart=js-1,je+1,njblock


can you block the i-dimension please?

edoyango · 2026-08-13T23:26:55Z

+    if (nkblock < nz) then
+      ! The early-exit test below compares htot against big_H on the host, and big_H has been
+      ! computed on the device, so fetch it once.  With a single k block that test never runs.
+      !$omp target update from(big_H)
+    endif


Suggested change

if (nkblock < nz) then

! The early-exit test below compares htot against big_H on the host, and big_H has been

! computed on the device, so fetch it once. With a single k block that test never runs.

!$omp target update from(big_H)

endif

! The early-exit test below compares htot against big_H on the host, and big_H has been

! computed on the device, so fetch it once. With a single k block that test never runs.

!$omp target update from(big_H) if (nkblock<nz)

edoyango · 2026-08-13T23:29:38Z

-      htot(i,j) = 0.0 ; Rml_int(i) = 0.0 ; SpV_int(i) = 0.0
+  do concurrent (k=1:nz, j=js-1:je+1, i=is-1:ie+1)
+    vol_dt_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-GV%Angstrom_H), 0.0)
+  enddo


This look should probably be left inside the blocked loop. Moving it outside probably reduces cache locality and does extra work (previously the early exist would've meant not all k values get populated). Might be why tmax is higher?

edoyango · 2026-08-13T23:30:29Z

+        kend = min(kstart+nkblock-1, nz)
+        EOSdom3(3,:) = [1, kend-kstart+1]
+        call calculate_density(tv%T(:,jstart:jend,kstart:kend), tv%S(:,jstart:jend,kstart:kend), &
+                               p_blk, rho_blk, tv%eqn_of_state, EOSdom3)


can you pass the whole T/S instead of slices, and instead update EOSdom(2,:) appropriately?

edoyango · 2026-08-13T23:34:05Z

-      dmu(k) = muza - muzb                  ! Change in mu(z) across layer [nondim]
-      ! dmu(k)*psi_mag is the transport in this layer [L2 H T-1 ~> m3 s-1]
+      muzb = mu(sigint, CS%MLE_tail_dh)       ! mu(z/MLD) for lower interface [nondim]
+      dmu = muza - muzb                     ! Change in mu(z) across layer [nondim]


can you keep dmu a vector and remove the 3d uhml?

edoyango · 2026-08-13T23:36:40Z

+    !$omp target update from(vhml)
+  endif
+
+  !$omp target exit data map(from: little_h, big_H, wpup, htot, buoy_av, uDml_diag, vDml_diag)


do all these have to be map(from)? can you guard them with if to reduce transfers?

edoyango · 2026-08-13T23:37:15Z

+    !$omp target update to(CS)
+    !$omp target enter data map(to: CS%Cr_space)
+    !$omp target enter data map(to: CS%MLD_filtered, CS%MLD_filtered_slow, CS%wpup_filtered)
+  endif


use inline if please

edoyango · 2026-08-13T23:41:18Z

                 G%HI, haloshift=1, unscale=GV%H_to_mks)
  endif

+  !$omp target enter data map(to: U_star_2d, h_MLD)


I think this interface for find_ustar is ported, so you can map(alloc: U_star_2d)

edoyango · 2026-08-13T23:42:42Z

+  enddo

  ! Estimate w'u' at h-points, with a floor to avoid division by zero later.
+  ! needed here, otehrwe death


Suggested change

! needed here, otehrwe death

edoyango · 2026-08-13T23:43:37Z

-      wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 + CS%nstar * w_star3))**2, CS%min_wstar2 ) &
+    do concurrent (j=js-1:je+1, i=is-1:ie+1)
+      wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 &
+                                 + CS%nstar * (max(0., -bflux(i,j)) * BLD(i,j))))**2, CS%min_wstar2 ) &


was there any reason for the code changes here? can you leave it as-is?

JorgeG94 · 2026-08-14T05:32:33Z

it seems some of the review fixes have affected CPU performance, it seems to be now stuck at

                                                        hits          tmin          tmax          tavg     
(Ocean mixed layer restrat)             60      0.257507      0.552579      0.323248

versus previously it was

                                                       hits          tmin          tmax          tavg        
(Ocean mixed layer restrat)             60      0.244403      0.477632      0.309130

with some per run noise, versus GDFL:

                                                       hits          tmin          tmax          tavg        
(Ocean mixed layer restrat)             60      0.302732      0.475717      0.366187

will investigate exactly what is causing it. I-Blocking is not helping

JorgeG94 · 2026-08-14T06:02:24Z

I think there's too much noise on Gadi if I run sufficiently things start to vary quite a bit. Will investigate...

edoyango · 2026-08-14T06:26:18Z

shrinking the tmp arrays to be block-sized where possible should help

Rewrite of mixedlayer_restrat_Bodner as do concurrent loops inside a single persistent device data region, replacing the earlier naive port. * The mixed-layer density integral evaluates the equation of state a block of MLE_NKBLOCK layers at a time rather than hoisting it to a whole-column 3D call. The block scratch is (SZI,SZJ,nkblock), so the CPU default of 1 keeps it to a single slab; the GPU default of 0 puts the column in one block. The row-at-a-time early exit becomes a slab-wide one, which only changes work that the htot < big_H guard already discarded. * No shadow copies of CS/GV/US scalars or of the input fields. The kernels read CS% and tv% directly, with those parents mapped alongside their components as nvfortran requires once any component is mapped. * One target data region covers the whole calculation, so nothing round-trips between phases. Arrays that only one branch or diagnostic touches (the block scratch, uhml/vhml, the host-side big_H and wpup variants, the chksums) are mapped or synced inside that branch. * dmu is stashed in uhml/vhml instead of a per-column array, which do concurrent cannot give each iteration without a per-thread allocation. mu is now pure so it can be called from a concurrent loop. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

h_MLD, BLD and bflux were left to implicit copyin, which the mapped-component switch disables along with derived-type copyin, so the kernels read uninitialized device memory and little_h/wpup came back NaN. BLD and bflux are mapped around the w'u' branches that dereference them; h_MLD joins the outer region. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Mapping a derived-type component makes nvfortran create a parent entry that supersedes an explicit map(to: CS) in the same directive. The component descriptors still attach -- CS%Cr_space read correctly -- but the struct body is never copied, so CS%mstar, CS%nstar and CS%min_wstar2 silently read as zero. wpup then collapsed to exactly zero, r_wpup went infinite, and psi_mag came out 0*Inf = NaN. Mapping the parents in their own enclosing region copies the struct body. benchmark_ALE now reproduces repro_ocean.stats bitwise. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

With MLE_NKBLOCK=1 the equation-of-state scratch was still a full (SZI,SZJ) slab, where upstream used a single row rho_ml(SZI). That loses the locality the original had, where the EOS fills a row and the integral consumes it straight back out of cache. Add MLE_NJBLOCK alongside MLE_NKBLOCK, following PGF_PLM_NJBLOCK: 1 and 1 on the CPU reduce the scratch to the one row upstream used, 0 and 0 on the GPU make it the whole halo-1 j domain by the whole column, which is a single kernel. Rml_int goes back to being the 1D host array upstream declares, with the tiled path carrying its running totals in Rml_blk. The early-exit test compares htot against big_H on the host, so big_H is fetched once when it can actually run. Without that it read a stale big_H, concluded every column was full and stopped calling the equation of state early, which truncated the integral -- benchmark_ALE energy moved in the 7th digit at NJBLOCK=NKBLOCK=1 while the untiled default stayed bitwise. benchmark_ALE reproduces repro_ocean.stats bitwise at both 0/0 and 1/1. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Dropping it reproduces the zeroed-scalar NaN even with the parent in its own enclosing region. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

The previous form leaned on `map(always, to: CS, tv)` in an enclosing structured region, a construct used nowhere else in MOM6. Replace it with what the rest of the GPU port does: unstructured `target enter data`, the control structure's array components mapped once in the init routine as hor_visc_init and friends do, and no structured `target data` regions anywhere. CS's scalars cannot be mapped alongside its arrays. Mapping any component of a derived type costs that type its implicit copyin, after which its scalars read as zero on the device while its array components still read correctly -- silent, and it surfaces only as a NaN much further downstream. So the eight scalars the kernels need are hoisted to plain locals, which is the practice the rest of the port already follows. Leaving the arrays unmapped instead is not an option: that faults with an illegal address on the first CS%MLD_filtered_slow read. benchmark_ALE reproduces repro_ocean.stats bitwise at MLE_NJBLOCK/NKBLOCK of both 0/0 and 1/1. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

The solo driver puts the whole MOM_control_struct on the device with `enter data map(alloc: MOM_CSp)`: the address range is claimed but its contents are never copied. Any later map(to:) of anything inside that range -- including the implicit copyin every kernel gets -- is then a present-table no-op, so kernels read zeros for the scalars of any control structure held as an inline member. Mapped array components were unaffected because their maps copy the heap data and attach the descriptor into the parent's device block; scalars have no attach equivalent. This was the mechanism behind the l_*/tau_* hoists. `target update` always copies regardless of presence, so issue `target update to(CS)` in mixedlayer_restrat_init once the parameters are set, BEFORE the component maps whose attach then fixes this struct's array descriptors in the device copy. With that, delete the eight hoisted locals and read CS% directly in the kernels, like every other module in the tree. Verified bitwise against repro_ocean.stats at both MLE_NJBLOCK/NKBLOCK 0/0 and 1/1 (overrides confirmed applied in MOM_parameter_doc.all). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Moves the BLD/bflux device mapping from mid-subroutine (with its stray placement note) into the enter-data cluster at the top of mixedlayer_restrat_Bodner, and releases them in the tail cluster with the other persistent arrays. These arrays are mapped nowhere else, so the mapping itself must stay; only its placement changes. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

find_ustar's mechanical-forcing path runs as do concurrent, so when U_star_2d is device-resident before the call the result is computed in place and the host-to-device copy after the call is unnecessary. The host copy is now only refreshed where it is actually read: the debug checksum and the id_ustar diagnostic. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Replaces the block-if around the initialization maps with if clauses on the directives themselves, and documents why the update to(CS) must precede the component maps: the driver maps the bare control structure without copying its contents, so the scalar components have to be pushed before array components are attached. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

The ported loop had inlined the w_star3 temporary with explicit parentheses reproducing its evaluation order. Restoring the named temporary as a do-concurrent local keeps the source identical to upstream while evaluating the exact same expression tree, so the answers are unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Replaces the unconditional exit-data map(from:) of seven scratch arrays with per-array target updates guarded by the diagnostic ids that actually consume them on the host, converts the uhml/vhml guards to inline if clauses, and releases all scratch without copying. The update from(h, uhtr, vhtr) moves out of mixedlayer_restrat_Bodner into the mixedlayer_restrat dispatch wrapper: the wrapper knows which scheme ran, and the host copies only need refreshing on the Bodner (device-side) path before the caller's halo update. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

The uhml/vhml stash existed only because do concurrent local() of an automatic array (bound from a dummy) crashes nvfortran with an illegal address during kernel execution (still present in 26.3). OpenMP private() handles the same array correctly, so the two velocity loops become omp target teams loop with dmu(SZK) private and their bodies restored verbatim from the original single-loop form: mu chain into dmu(:), CFL limiting, then uhml = dmu(k)*psi_mag. Bitwise-neutral: the stash stored dmu exactly and applied one multiply; the restored form applies the same single multiply. These loops can return to do concurrent when the compiler bug is fixed. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Three review items on the same loop nest, applied together so the walk reads coherently: - vol_dt_avail is folded back into the integral walk in both branches, restoring the original loop fusion (evaluated unconditionally ahead of the early-exit guard, exactly as the original computed it for every layer; the host fallback pushes it to the device with its other results). - The equation-of-state calls take the whole tv%T/tv%S with the tile expressed as absolute index ranges in EOSdom3, so no argument slices are made; the tile scratch carries full ij extent so the same ranges index it. - The i dimension is blocked: MLE_NIBLOCK joins MLE_NJBLOCK and MLE_NKBLOCK, with the i tile selected through EOSdom3(1,:) exactly as the j tile is. Defaults are 0 (whole i domain) on both CPU and GPU builds: benchmarks at these row lengths showed per-tile call overhead outweighing any cache benefit, and the runtime parameter remains for configurations whose rows outgrow cache. Per-column arithmetic and accumulation order are unchanged by tiling, so the answers are identical for every block-size choice. Also replaces a stale comment on the big_H branches with one that says which variants run on the host. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

The dev/gpu merge added a WAVE_ENHANCED_USTAR scaling of U_star_2d in mixedlayer_restrat_Bodner. It lands after find_ustar, which this port runs on the device, so the new host loop would read stale host memory and its result would never reach the kernels that consume U_star_2d. Lam2 is an unmapped host pointer, so keep the scaling on the host and round-trip U_star_2d around it, guarded on the flag since it defaults to false.

JorgeG94 marked this pull request as ready for review August 10, 2026 01:12

edoyango suggested changes Aug 13, 2026

View reviewed changes

JorgeG94 force-pushed the feat/bodner_mle branch 3 times, most recently from cc26279 to 15d0c8e Compare August 18, 2026 01:15

JorgeG94 and others added 19 commits August 31, 2026 21:59

Trim the map-clause commentary to the load-bearing line

addd99e

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Record that always is load-bearing on the CS/tv parent map

4dba21f

Dropping it reproduces the zeroed-scalar NaN even with the parent in its own enclosing region. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Drop leftover scaffolding from MOM.F90

9d6d039

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

fix fix

9beb8eb

fix doxygne

2a8ccd7

JorgeG94 force-pushed the feat/bodner_mle branch from 15d0c8e to 76f3d0e Compare September 1, 2026 05:54

Uh oh!

Conversation

JorgeG94 commented Jul 28, 2026

Uh oh!

github-actions Bot commented Jul 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

GPU Port Coverage

Uh oh!

JorgeG94 commented Aug 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

edoyango left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JorgeG94 commented Aug 14, 2026

Uh oh!

JorgeG94 commented Aug 14, 2026

Uh oh!

edoyango commented Aug 14, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

github-actions Bot commented Jul 29, 2026 •

edited

Loading

JorgeG94 commented Aug 10, 2026 •

edited

Loading