diff --git a/src/megatron/bridge/models/qwen_omni/qwen3_omni_step.py b/src/megatron/bridge/models/qwen_omni/qwen3_omni_step.py
index ac6efda805..fa0688711d 100644
--- a/src/megatron/bridge/models/qwen_omni/qwen3_omni_step.py
+++ b/src/megatron/bridge/models/qwen_omni/qwen3_omni_step.py
@@ -27,6 +27,7 @@
 from megatron.bridge.training.losses import (
     create_masked_next_token_loss_function as _create_loss_function,
 )
+from megatron.bridge.training.utils.flop_utils import accumulate_flops_metadata
 from megatron.bridge.training.utils.padding_utils import (
     pad_or_truncate_2d_to_len,
     pad_or_truncate_attn_to_len,
@@ -181,6 +182,16 @@ def forward_step(
     if pg_collection.cp.size() > 1:
         raise NotImplementedError("Qwen3-Omni training supports SP/EP, but CP is not supported yet.")
 
+    # Accumulate FLOPS metadata across micro-batches. Qwen3-Omni does not pack
+    # within a batch, so cu_seqlens is absent and the helper falls back to
+    # BSHD math for the attention term. Vision-patch tracking still applies.
+    accumulate_flops_metadata(
+        state,
+        tokens,
+        image_grid_thw=multimodal_inputs.get("image_grid_thw") if isinstance(multimodal_inputs, dict) else None,
+        video_grid_thw=multimodal_inputs.get("video_grid_thw") if isinstance(multimodal_inputs, dict) else None,
+    )
+
     forward_args = {
         "input_ids": tokens,
         "position_ids": position_ids,
diff --git a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/model.py b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/model.py
index 3ee0309f25..2ac3bb772b 100644
--- a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/model.py
+++ b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/model.py
@@ -51,6 +51,58 @@
 )
 from megatron.bridge.models.qwen_vl.modelling_qwen3_vl.vision_model import Qwen3VLVisionModel
 
+try:
+    import transformer_engine_torch as tex
+except ImportError:
+    tex = None
+
+
+def _compact_thd_cp_index(
+    packed_seq_params: PackedSeqParams,
+    total_tokens: int,
+    cp_size: int,
+    cp_rank: int,
+) -> torch.Tensor | None:
+    """Return the CP-local token indices for an already-packed THD stream."""
+    if cp_size <= 1:
+        return None
+    if tex is None:
+        raise RuntimeError("QWEN3VL_THD_COMPACT_PACKING with CP>1 requires transformer_engine_torch")
+    return tex.thd_get_partitioned_indices(
+        packed_seq_params.cu_seqlens_q_padded,
+        total_tokens,
+        cp_size,
+        cp_rank,
+    )
+
+
+def _pack_position_ids_to_compact_thd(
+    position_ids: torch.Tensor,
+    packed_seq_params: PackedSeqParams,
+) -> torch.Tensor:
+    """Pack BSHD MRoPE position IDs into the same compact THD layout as input_ids.
+
+    Qwen3-VL's ``get_rope_index`` understands per-sample vision placeholder
+    order best in BSHD form. Compact THD computes those BSHD position IDs first,
+    then copies each real sample span into the physical padded THD segment.
+    Alignment padding keeps the default zero position IDs and is masked later.
+    """
+    real_lengths = packed_seq_params.cu_seqlens_q[1:] - packed_seq_params.cu_seqlens_q[:-1]
+    cu_seqlens_padded = packed_seq_params.cu_seqlens_q_padded
+    total_padded_len = int(cu_seqlens_padded[-1].item())
+    packed_position_ids = torch.zeros(
+        position_ids.size(0),
+        1,
+        total_padded_len,
+        dtype=position_ids.dtype,
+        device=position_ids.device,
+    )
+    for batch_idx in range(real_lengths.numel()):
+        real_len = int(real_lengths[batch_idx].item())
+        start = int(cu_seqlens_padded[batch_idx].item())
+        packed_position_ids[:, 0, start : start + real_len] = position_ids[:, batch_idx, :real_len]
+    return packed_position_ids
+
 
 class Qwen3VLModel(MegatronModule):
     """Qwen3VL multi-modal model.
@@ -352,6 +404,10 @@ def forward(
         inference_context: object | None = None,
         runtime_gather_output: bool | None = None,
         mm_token_type_ids: torch.Tensor = None,
+        moe_padding_mask: torch.Tensor = None,
+        qwen3vl_thd_compact_packing: bool = False,
+        qwen3vl_compact_input_ids_bshd: torch.Tensor = None,
+        qwen3vl_compact_attention_mask_bshd: torch.Tensor = None,
         **kwargs,
     ) -> torch.Tensor:
         """Forward function of the Qwen3VL model.
@@ -399,6 +455,19 @@ def forward(
         # so it must be a real tensor. For packed sequences we use the THD-format
         # input_ids_thd (updated below); for regular sequences we use input_ids as-is.
         lm_input_ids = input_ids
+        # In compact THD mode qwen3_vl_step has already changed input_ids from
+        # BSHD [B, S] to THD [1, T]. Cache the full physical length before any
+        # CP-local index_select so we can avoid double-splitting labels/masks.
+        compact_total_tokens = input_ids.size(1) if qwen3vl_thd_compact_packing else None
+        compact_cp_index = None
+        if qwen3vl_thd_compact_packing:
+            assert packed_seq_params is not None, "QWEN3VL_THD_COMPACT_PACKING requires packed_seq_params"
+            compact_cp_index = _compact_thd_cp_index(
+                packed_seq_params,
+                compact_total_tokens,
+                cp_size,
+                cp_rank,
+            )
 
         if self.pre_process:
             # can reorganize_inputs at dataset
@@ -503,7 +572,10 @@ def forward(
 
             if combined_embeddings is not None and cp_size > 1 and packed_seq_params is None:
                 combined_embeddings = split_data_cp_rank(combined_embeddings, cp_size, 0, cp_rank)
-            if packed_seq_params is not None:
+            if packed_seq_params is not None and not qwen3vl_thd_compact_packing:
+                # Legacy Qwen3-VL THD path: model.forward receives padded BSHD
+                # input_ids so it can do vision/text combine first, then this
+                # helper converts BSHD -> CP-local THD.
                 if attention_mask is None:
                     attention_mask = torch.ones_like(input_ids, dtype=torch.bool, device=input_ids.device)
                 input_ids_thd, _ = preprocess_packed_seqs(
@@ -551,7 +623,7 @@ def forward(
                 )
                 combined_embeddings = combined_embeddings_thd
 
-            if self.config.sequence_parallel:
+            if self.config.sequence_parallel and not qwen3vl_thd_compact_packing:
                 combined_embeddings = tensor_parallel.scatter_to_sequence_parallel_region(combined_embeddings)
                 combined_embeddings = combined_embeddings.contiguous()
 
@@ -559,13 +631,71 @@ def forward(
             combined_embeddings = None
             # On non-pre_process PP stages (e.g. the last stage where MTP runs),
             # convert lm_input_ids to THD format so it matches position_ids.
-            if packed_seq_params is not None:
+            if packed_seq_params is not None and not qwen3vl_thd_compact_packing:
+                # Same legacy conversion is needed on later PP stages because MTP
+                # reads input_ids even though embeddings were produced upstream.
                 if attention_mask is None:
                     attention_mask = torch.ones_like(input_ids, dtype=torch.bool, device=input_ids.device)
                 lm_input_ids, _ = preprocess_packed_seqs(
                     input_ids, attention_mask, pre_process=True, pg_collection=self.pg_collection
                 )
 
+        if qwen3vl_thd_compact_packing:
+            # Compact path: input_ids/combined_embeddings are already global THD.
+            # Keep the old invariant "combine before CP split", but use TE's THD
+            # partition index instead of preprocess_packed_seqs' BSHD conversion.
+            if compact_cp_index is not None:
+                if self.pre_process and combined_embeddings is not None:
+                    full_vision_mask = vision_mask
+                    local_vision_mask = (
+                        full_vision_mask.index_select(1, compact_cp_index) if full_vision_mask is not None else None
+                    )
+                    if deepstack_feature_lists is not None and full_vision_mask is not None:
+                        # Deepstack features are produced only at visual token
+                        # positions. Re-materialize them in full THD layout, CP
+                        # slice, then gather back the local visual positions.
+                        full_embeddings_bsh = combined_embeddings.transpose(0, 1).contiguous()
+                        new_deepstack_feature_lists = []
+                        for deepstack_visual_embed in deepstack_feature_lists:
+                            tmp_embeddings = torch.zeros_like(full_embeddings_bsh)
+                            tmp_embeddings[full_vision_mask] = deepstack_visual_embed
+                            tmp_embeddings = tmp_embeddings.index_select(1, compact_cp_index)
+                            new_deepstack_feature_lists.append(tmp_embeddings[local_vision_mask].contiguous())
+                        deepstack_feature_lists = new_deepstack_feature_lists
+                    combined_embeddings = combined_embeddings.index_select(0, compact_cp_index)
+                    vision_mask = local_vision_mask
+
+                # input_ids is passed to language_model for MTP/future-token
+                # embedding paths; it must match the CP-local decoder_input.
+                input_ids = input_ids.index_select(1, compact_cp_index)
+                lm_input_ids = input_ids
+                if labels is not None and labels.dim() >= 2 and labels.size(1) == compact_total_tokens:
+                    labels = labels.index_select(1, compact_cp_index)
+                if loss_mask is not None and loss_mask.dim() >= 2 and loss_mask.size(1) == compact_total_tokens:
+                    loss_mask = loss_mask.index_select(1, compact_cp_index)
+                if (
+                    moe_padding_mask is not None
+                    and moe_padding_mask.dim() >= 2
+                    and moe_padding_mask.size(1) == compact_total_tokens
+                ):
+                    moe_padding_mask = moe_padding_mask.index_select(1, compact_cp_index)
+
+            if self.pre_process and self.config.sequence_parallel:
+                # Legacy path scatters after BSHD->THD conversion. Compact path
+                # skips that block, so SP scatter has to happen here instead.
+                combined_embeddings = tensor_parallel.scatter_to_sequence_parallel_region(combined_embeddings)
+                combined_embeddings = combined_embeddings.contiguous()
+            if self.config.sequence_parallel and moe_padding_mask is not None:
+                tp_size = self.pg_collection.tp.size()
+                if tp_size > 1:
+                    tp_rank = self.pg_collection.tp.rank()
+                    assert moe_padding_mask.size(1) % tp_size == 0, (
+                        "compact THD MoE padding mask must be divisible by TP size under sequence parallelism"
+                    )
+                    seq_len_per_tp = moe_padding_mask.size(1) // tp_size
+                    start = tp_rank * seq_len_per_tp
+                    moe_padding_mask = moe_padding_mask[:, start : start + seq_len_per_tp].contiguous()
+
         visual_pos_masks = vision_mask
         deepstack_visual_embeds = deepstack_feature_lists
         if self.config.sequence_parallel or cp_size > 1:
@@ -590,7 +720,30 @@ def forward(
                     sequence_parallel=self.config.sequence_parallel,
                 )
 
-        if position_ids is None:
+        if qwen3vl_thd_compact_packing:
+            if position_ids is None:
+                assert qwen3vl_compact_input_ids_bshd is not None, "compact THD packing requires original BSHD input_ids"
+                # MRoPE must see per-sample BSHD order, otherwise [1, T] would
+                # look like one long sample and image/video grids would be
+                # consumed with the wrong boundaries.
+                position_ids, _ = get_rope_index(
+                    self.config.spatial_merge_size,
+                    self.image_token_id,
+                    self.video_token_id,
+                    self.vision_start_token_id,
+                    qwen3vl_compact_input_ids_bshd,
+                    image_grid_thw=image_grid_thw,
+                    video_grid_thw=video_grid_thw,
+                    attention_mask=qwen3vl_compact_attention_mask_bshd,
+                )
+                position_ids = _pack_position_ids_to_compact_thd(position_ids, packed_seq_params)
+            if compact_cp_index is not None:
+                position_ids = position_ids.index_select(2, compact_cp_index)
+            # THD attention receives boundaries via packed_seq_params. A dense
+            # mask would describe BSHD layout and is intentionally disabled.
+            attention_mask = None
+            self.language_model.rotary_pos_emb.is_thd_format = True
+        elif position_ids is None:
             # BSHD
             # Megatron uses 4D bool masks ([B|1,1,S,S], True=masked); HF uses 2D keep masks ([B,S], 1=keep)
             # For simplicity, we set hf_attention_mask to None.
@@ -630,6 +783,7 @@ def forward(
             decoder_input=combined_embeddings,  # only not None in the first decoder PP stage
             labels=labels,  # only not None in the last decoder PP stage
             loss_mask=loss_mask,  # Added for THD training compatibility
+            padding_mask=moe_padding_mask,
             inference_params=inference_params,  # currently always None
             packed_seq_params=packed_seq_params,  # currently always None
             runtime_gather_output=runtime_gather_output,
diff --git a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/text_model.py b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/text_model.py
index d067655a1b..73610fb8d3 100644
--- a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/text_model.py
+++ b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/text_model.py
@@ -120,6 +120,9 @@ def forward(
         *,
         inference_params: Optional[BaseInferenceContext] = None,
         loss_mask: Optional[Tensor] = None,
+        # Compact THD inserts small alignment padding between samples. Thread
+        # this mask to decoder layers so MoE router losses/stats exclude it.
+        padding_mask: Optional[Tensor] = None,
         # args for deepstack
         visual_pos_masks: Optional[torch.Tensor] = None,
         deepstack_visual_embeds: Optional[list[torch.Tensor]] = None,
@@ -170,6 +173,7 @@ def forward(
             # the standard components only.
             packed_seq_params=packed_seq_params,
             sequence_len_offset=sequence_len_offset,
+            padding_mask=padding_mask,
             visual_pos_masks=visual_pos_masks,
             deepstack_visual_embeds=deepstack_visual_embeds,
             **(extra_block_kwargs or {}),
diff --git a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/transformer_block.py b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/transformer_block.py
index 6284617bed..9a720d37ba 100644
--- a/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/transformer_block.py
+++ b/src/megatron/bridge/models/qwen_vl/modelling_qwen3_vl/transformer_block.py
@@ -504,6 +504,9 @@ def _checkpointed_forward(
         attention_bias: Tensor,
         packed_seq_params: PackedSeqParams,
         use_inner_fp8_context: bool,
+        # Must be a checkpoint input, not a closure-only value, so activation
+        # recompute replays the exact same MoE padding mask as the first forward.
+        padding_mask: Optional[Tensor] = None,
         # args for deepstack
         visual_pos_masks: Optional[torch.Tensor] = None,
         deepstack_visual_embeds: Optional[list[torch.Tensor]] = None,
@@ -517,6 +520,7 @@ def custom_forward(
                 context,
                 context_mask,
                 rotary_pos_emb,
+                padding_mask,
                 visual_pos_masks,
                 *deepstack_visual_embeds_args,
             ):
@@ -529,6 +533,9 @@ def custom_forward(
                         else nullcontext()
                     )
                     with inner_fp8_context:
+                        # MCore TransformerLayer passes padding_mask to MoE
+                        # router loss/stat paths. Compact THD uses it so
+                        # alignment padding does not skew load-balancing state.
                         hidden_states, context = layer(
                             hidden_states=hidden_states,
                             attention_mask=attention_mask,
@@ -538,6 +545,7 @@ def custom_forward(
                             attention_bias=attention_bias,
                             inference_context=None,
                             packed_seq_params=packed_seq_params,
+                            padding_mask=padding_mask,
                         )
 
                         if self.pre_process and deepstack_visual_embeds is not None:
@@ -567,6 +575,7 @@ def checkpoint_handler(forward_func):
                     context,
                     context_mask,
                     rotary_pos_emb,
+                    padding_mask,
                     visual_pos_masks,
                     *deepstack_visual_embeds_tuple,
                 )
@@ -579,6 +588,7 @@ def checkpoint_handler(forward_func):
                     context,
                     context_mask,
                     rotary_pos_emb,
+                    padding_mask,
                     visual_pos_masks,
                     *deepstack_visual_embeds_tuple,
                 )
@@ -618,6 +628,7 @@ def checkpoint_handler(forward_func):
                         context,
                         context_mask,
                         rotary_pos_emb,
+                        padding_mask,
                         visual_pos_masks,
                         *deepstack_visual_embeds_tuple,
                     )
@@ -639,6 +650,7 @@ def forward(
         inference_context: Optional[BaseInferenceContext] = None,
         packed_seq_params: Optional[PackedSeqParams] = None,
         sequence_len_offset: Optional[Tensor] = None,
+        padding_mask: Optional[Tensor] = None,
         *,
         inference_params: Optional[BaseInferenceContext] = None,
         # args for deepstack
@@ -731,6 +743,7 @@ def forward(
                     attention_bias=attention_bias,
                     packed_seq_params=packed_seq_params,
                     use_inner_fp8_context=use_inner_fp8_context,
+                    padding_mask=padding_mask,
                     visual_pos_masks=visual_pos_masks,
                     deepstack_visual_embeds=deepstack_visual_embeds,
                 )
@@ -754,6 +767,7 @@ def forward(
                             inference_context=inference_context,
                             packed_seq_params=packed_seq_params,
                             sequence_len_offset=sequence_len_offset,
+                            padding_mask=padding_mask,
                         )
 
                         if self.pre_process and deepstack_visual_embeds is not None:
diff --git a/src/megatron/bridge/models/qwen_vl/qwen3_vl_step.py b/src/megatron/bridge/models/qwen_vl/qwen3_vl_step.py
index c87202262b..3ebd51b2c1 100644
--- a/src/megatron/bridge/models/qwen_vl/qwen3_vl_step.py
+++ b/src/megatron/bridge/models/qwen_vl/qwen3_vl_step.py
@@ -13,6 +13,7 @@
 # limitations under the License.
 import logging
 import math
+import os
 from functools import partial
 from typing import Any, Iterable
 
@@ -27,6 +28,7 @@
     create_masked_next_token_loss_function as _create_loss_function,
 )
 from megatron.bridge.training.state import GlobalState
+from megatron.bridge.training.utils.flop_utils import accumulate_flops_metadata, accumulate_token_throughput_metadata
 from megatron.bridge.training.utils.padding_utils import (
     pad_or_truncate_2d_to_len,
     pad_or_truncate_attn_to_len,
@@ -37,6 +39,11 @@
 
 logger = logging.getLogger(__name__)
 
+try:
+    import transformer_engine_torch as tex
+except ImportError:
+    tex = None
+
 
 def get_batch_from_iterator(
     data_iterator: Iterable,
@@ -94,6 +101,10 @@ def get_batch_from_iterator(
         else:
             _batch_required_keys[key] = None
 
+    raw_attn = batch.get("attention_mask")
+    if isinstance(raw_attn, torch.Tensor) and raw_attn.dim() == 2:
+        _batch_required_keys["_padding_mask"] = raw_attn.cuda(non_blocking=True)
+
     return _batch_required_keys
 
 
@@ -110,6 +121,7 @@ def get_batch(
     torch.Tensor,
     torch.Tensor,
     torch.Tensor,
+    torch.Tensor,
     Any,
 ]:
     """Generate a batch.
@@ -147,6 +159,7 @@ def get_batch(
         batch.get("loss_mask"),
         batch.get("attention_mask"),
         batch.get("position_ids"),
+        batch.get("_padding_mask"),
         multi_modal_inputs,
     )
 
@@ -207,6 +220,115 @@ def pack_or_pad_batch_sequences(
     return tokens, labels, loss_mask, attention_mask, position_ids, packed_seq_params
 
 
+def _get_real_sequence_lengths(
+    tokens: torch.Tensor,
+    padding_mask: torch.Tensor | None,
+) -> torch.Tensor:
+    """Return per-sample non-pad lengths. The 2D mask follows HF convention: 1/True means keep."""
+    if padding_mask is None:
+        return torch.full((tokens.size(0),), tokens.size(1), dtype=torch.int32, device=tokens.device)
+    if padding_mask.dim() != 2:
+        raise ValueError(f"QWEN3VL_THD_COMPACT_PACKING expects a 2D padding mask, got {tuple(padding_mask.shape)}")
+    real_lengths = padding_mask.to(dtype=torch.int32, device=tokens.device).sum(dim=1)
+    return real_lengths.clamp(min=1, max=tokens.size(1))
+
+
+def _build_compact_packed_seq_params(
+    real_lengths: torch.Tensor,
+    pad_to_multiple_of: int,
+) -> tuple[PackedSeqParams, torch.Tensor]:
+    """Build THD metadata where real and physical padded boundaries are both preserved.
+
+    ``cu_seqlens_q`` describes the true sample lengths used for attention/FLOPS
+    semantics. ``cu_seqlens_q_padded`` describes the physical packed layout used
+    by TransformerEngine's THD CP partitioner.
+    """
+    padded_lengths = torch.div(
+        real_lengths + pad_to_multiple_of - 1,
+        pad_to_multiple_of,
+        rounding_mode="floor",
+    ) * pad_to_multiple_of
+    cu_seqlens = torch.zeros(real_lengths.numel() + 1, dtype=torch.int32, device=real_lengths.device)
+    cu_seqlens[1:] = torch.cumsum(real_lengths.to(torch.int32), dim=0)
+    cu_seqlens_padded = torch.zeros_like(cu_seqlens)
+    cu_seqlens_padded[1:] = torch.cumsum(padded_lengths.to(torch.int32), dim=0)
+    max_seqlen = int(padded_lengths.max().item())
+    return (
+        PackedSeqParams(
+            qkv_format="thd",
+            cu_seqlens_q=cu_seqlens,
+            max_seqlen_q=max_seqlen,
+            cu_seqlens_kv=cu_seqlens,
+            max_seqlen_kv=max_seqlen,
+            cu_seqlens_q_padded=cu_seqlens_padded,
+            cu_seqlens_kv_padded=cu_seqlens_padded,
+        ),
+        padded_lengths,
+    )
+
+
+def _compact_pack_batch_sequences(
+    tokens: torch.Tensor,
+    labels: torch.Tensor | None,
+    loss_mask: torch.Tensor | None,
+    real_lengths: torch.Tensor,
+    padded_lengths: torch.Tensor,
+    pad_token_id: int = 0,
+) -> tuple[torch.Tensor, torch.Tensor | None, torch.Tensor | None, torch.Tensor]:
+    """Pack BSHD rows directly into one THD stream without first padding to seq_length.
+
+    The returned ``moe_padding_mask`` marks only the alignment padding inserted
+    between packed samples. It is not an attention mask; MCore MoE uses it to
+    exclude fake padding from router losses and expert-bias token statistics.
+    """
+    total_padded_len = int(padded_lengths.sum().item())
+    packed_tokens = torch.full((1, total_padded_len), pad_token_id, dtype=tokens.dtype, device=tokens.device)
+    packed_labels = (
+        torch.full((1, total_padded_len), -100, dtype=labels.dtype, device=labels.device) if labels is not None else None
+    )
+    packed_loss_mask = (
+        torch.zeros((1, total_padded_len), dtype=loss_mask.dtype, device=loss_mask.device)
+        if loss_mask is not None
+        else None
+    )
+    moe_padding_mask = torch.zeros((1, total_padded_len), dtype=torch.bool, device=tokens.device)
+
+    offset = 0
+    for batch_idx in range(tokens.size(0)):
+        real_len = int(real_lengths[batch_idx].item())
+        padded_len = int(padded_lengths[batch_idx].item())
+        packed_tokens[0, offset : offset + real_len] = tokens[batch_idx, :real_len]
+        if packed_labels is not None:
+            packed_labels[0, offset : offset + real_len] = labels[batch_idx, :real_len]
+        if packed_loss_mask is not None:
+            packed_loss_mask[0, offset : offset + real_len] = loss_mask[batch_idx, :real_len]
+        if padded_len > real_len:
+            moe_padding_mask[0, offset + real_len : offset + padded_len] = True
+        offset += padded_len
+
+    return packed_tokens, packed_labels, packed_loss_mask, moe_padding_mask
+
+
+def _get_compact_thd_cp_index(
+    packed_seq_params: PackedSeqParams,
+    total_tokens: int,
+    cp_size: int,
+    cp_rank: int,
+) -> torch.Tensor | None:
+    if cp_size <= 1:
+        return None
+    if tex is None:
+        raise RuntimeError("QWEN3VL_THD_COMPACT_PACKING with CP>1 requires transformer_engine_torch")
+    # THD CP must split each packed segment independently. TE's helper uses
+    # cu_seqlens_q_padded to apply the standard zigzag CP split per segment.
+    return tex.thd_get_partitioned_indices(
+        packed_seq_params.cu_seqlens_q_padded,
+        total_tokens,
+        cp_size,
+        cp_rank,
+    )
+
+
 def forward_step(
     state: GlobalState,
     data_iterator: Iterable,
@@ -242,6 +364,7 @@ def forward_step(
             loss_mask,
             attention_mask,
             position_ids,
+            padding_mask,
             multi_modal_inputs,
         ) = get_batch(data_iterator, state.cfg, use_mtp, is_first_pp_stage=is_first, is_last_pp_stage=is_last)
     timers("batch-generator").stop()
@@ -250,44 +373,132 @@ def forward_step(
     # Qwen3VL model need the original input and do cp and sp split in model.forward.
     pack_sequences_in_batch = getattr(state.cfg.dataset, "pack_sequences_in_batch", False)
 
-    tokens, labels, loss_mask, attention_mask, position_ids, packed_seq_params = pack_or_pad_batch_sequences(
-        tokens,
-        labels,
-        loss_mask,
-        attention_mask,
-        position_ids,
-        this_pg_collection,
-        use_fp8_padding=True,
-        force_to_pad_to_seq_len=this_pg_collection.pp.size() > 1 or this_pg_collection.ep.size() > 1,
-        seq_length=config.seq_length,
-    )
-    forward_args = {
-        "input_ids": tokens,
-        "labels": labels,
-        "loss_mask": loss_mask,
-        "attention_mask": attention_mask,
-        "position_ids": position_ids,
-    }
-
-    original_tokens = tokens.clone()
-    forward_args = get_batch_on_this_cp_rank(forward_args, cp_group=this_pg_collection.cp)
-    forward_args["packed_seq_params"] = None
-    forward_args["input_ids"] = original_tokens
-    # calculate position_ids in model forward
-    forward_args["position_ids"] = None
-    if pack_sequences_in_batch:
-        if forward_args["labels"] is not None:
-            # When using pp, labels could be None
-            forward_args["labels"] = forward_args["labels"].reshape(1, -1)
-        attention_mask = torch.ones(
-            original_tokens.shape[0], original_tokens.shape[1], dtype=torch.bool, device=original_tokens.device
+    compact_thd_packing = pack_sequences_in_batch and os.getenv("QWEN3VL_THD_COMPACT_PACKING", "0") == "1"
+
+    if compact_thd_packing:
+        # Contract with Qwen3VLModel.forward:
+        # - ``input_ids`` below becomes compact THD [1, total_padded_tokens].
+        # - Qwen3-VL still needs the original BSHD layout to compute MRoPE
+        #   and to interpret image/video placeholder order per sample.
+        compact_input_ids_bshd = tokens
+        compact_attention_mask_bshd = (
+            padding_mask
+            if padding_mask is not None
+            else torch.ones_like(compact_input_ids_bshd, dtype=torch.bool, device=compact_input_ids_bshd.device)
         )
-        forward_args["attention_mask"] = attention_mask
-        if forward_args["loss_mask"] is not None:
-            forward_args["loss_mask"] = forward_args["loss_mask"].reshape(1, -1)
-        # qwen3vl need the original input_ids and position_ids
-        # use split attention mask for calculate loss
-        forward_args["packed_seq_params"] = packed_seq_params
+        real_lengths = _get_real_sequence_lengths(tokens, compact_attention_mask_bshd)
+        cp_size = this_pg_collection.cp.size()
+        tp_size = this_pg_collection.tp.size()
+        # CP requires each packed segment to be divisible by 2*CP for zigzag
+        # splitting. Sequence parallelism further requires the CP-local chunk to
+        # be divisible across TP ranks. Keep 16 for FP8/TE-friendly alignment.
+        cp_multiple = 2 * cp_size if cp_size > 1 else 1
+        sp_multiple = cp_size * tp_size if getattr(config, "sequence_parallel", False) and tp_size > 1 else 1
+        pad_to_multiple_of = math.lcm(16, cp_multiple, sp_multiple)
+        packed_seq_params, padded_lengths = _build_compact_packed_seq_params(real_lengths, pad_to_multiple_of)
+        tokens, labels, loss_mask, moe_padding_mask = _compact_pack_batch_sequences(
+            tokens,
+            labels,
+            loss_mask,
+            real_lengths,
+            padded_lengths,
+        )
+        # Compact THD uses packed_seq_params for sequence boundaries, so there
+        # is no dense attention mask or precomputed BSHD position_ids to pass.
+        attention_mask = None
+        position_ids = None
+
+        accumulate_flops_metadata(
+            state,
+            tokens,
+            cu_seqlens_unpadded=packed_seq_params.cu_seqlens_q,
+            image_grid_thw=multi_modal_inputs.get("image_grid_thw") if isinstance(multi_modal_inputs, dict) else None,
+            video_grid_thw=multi_modal_inputs.get("video_grid_thw") if isinstance(multi_modal_inputs, dict) else None,
+        )
+        accumulate_token_throughput_metadata(
+            state,
+            real_tokens=int(real_lengths.sum().item()),
+            packed_tokens=int(padded_lengths.sum().item()),
+        )
+
+        cp_index = _get_compact_thd_cp_index(
+            packed_seq_params,
+            tokens.size(1),
+            this_pg_collection.cp.size(),
+            this_pg_collection.cp.rank(),
+        )
+        if cp_index is not None:
+            # The model CP-splits embeddings after vision/text combine. The
+            # loss tensors live in the step closure, so they must be split here
+            # with exactly the same THD index.
+            labels = labels.index_select(1, cp_index) if labels is not None else None
+            loss_mask = loss_mask.index_select(1, cp_index) if loss_mask is not None else None
+            moe_padding_mask = moe_padding_mask.index_select(1, cp_index)
+
+        forward_args = {
+            "input_ids": tokens,
+            "labels": labels,
+            "loss_mask": loss_mask,
+            "attention_mask": attention_mask,
+            "position_ids": position_ids,
+            "packed_seq_params": packed_seq_params,
+            "moe_padding_mask": moe_padding_mask,
+            "qwen3vl_thd_compact_packing": True,
+            # Extra internal metadata needed because compact THD removes the
+            # original batch dimension before model.forward runs MRoPE.
+            "qwen3vl_compact_input_ids_bshd": compact_input_ids_bshd,
+            "qwen3vl_compact_attention_mask_bshd": compact_attention_mask_bshd,
+        }
+    else:
+        tokens, labels, loss_mask, attention_mask, position_ids, packed_seq_params = pack_or_pad_batch_sequences(
+            tokens,
+            labels,
+            loss_mask,
+            attention_mask,
+            position_ids,
+            this_pg_collection,
+            use_fp8_padding=True,
+            force_to_pad_to_seq_len=this_pg_collection.pp.size() > 1 or this_pg_collection.ep.size() > 1,
+            seq_length=config.seq_length,
+        )
+
+        # Accumulate FLOPS metadata across micro-batches. When in-batch packing is
+        # active, cu_seqlens_q describes the sub-seq boundaries used by THD attention.
+        accumulate_flops_metadata(
+            state,
+            tokens,
+            cu_seqlens=getattr(packed_seq_params, "cu_seqlens_q", None) if packed_seq_params is not None else None,
+            image_grid_thw=multi_modal_inputs.get("image_grid_thw") if isinstance(multi_modal_inputs, dict) else None,
+            video_grid_thw=multi_modal_inputs.get("video_grid_thw") if isinstance(multi_modal_inputs, dict) else None,
+        )
+
+        forward_args = {
+            "input_ids": tokens,
+            "labels": labels,
+            "loss_mask": loss_mask,
+            "attention_mask": attention_mask,
+            "position_ids": position_ids,
+        }
+
+        original_tokens = tokens.clone()
+        forward_args = get_batch_on_this_cp_rank(forward_args, cp_group=this_pg_collection.cp)
+        forward_args["packed_seq_params"] = None
+        forward_args["input_ids"] = original_tokens
+        # calculate position_ids in model forward
+        forward_args["position_ids"] = None
+        if pack_sequences_in_batch:
+            if forward_args["labels"] is not None:
+                # When using pp, labels could be None
+                forward_args["labels"] = forward_args["labels"].reshape(1, -1)
+            attention_mask = torch.ones(
+                original_tokens.shape[0], original_tokens.shape[1], dtype=torch.bool, device=original_tokens.device
+            )
+            forward_args["attention_mask"] = attention_mask
+            if forward_args["loss_mask"] is not None:
+                forward_args["loss_mask"] = forward_args["loss_mask"].reshape(1, -1)
+            # qwen3vl need the original input_ids and position_ids
+            # use split attention mask for calculate loss
+            forward_args["packed_seq_params"] = packed_seq_params
 
     # use cp split loss mask for calculate loss
     loss_mask = forward_args["loss_mask"]
diff --git a/src/megatron/bridge/training/config.py b/src/megatron/bridge/training/config.py
index 1ed3633162..dd2f1ff6ee 100644
--- a/src/megatron/bridge/training/config.py
+++ b/src/megatron/bridge/training/config.py
@@ -1154,6 +1154,11 @@ def validate(self) -> None:
         if getattr(self.dataset, "pack_sequences_in_batch", False):
             self.model._pack_sequences_in_batch = True
 
+            # Local THD packing does not use Megatron Core's sequence_packing_scheduler,
+            # but HybridEP still needs equal token counts inside each dispatch group.
+            if getattr(self.model, "moe_flex_dispatcher_backend", None) == "hybridep":
+                self.model.moe_hybridep_pad_variable_tokens = True
+
         if hasattr(self.dataset, "finalize"):
             self.dataset.finalize()
         if hasattr(self.ddp, "finalize"):
diff --git a/src/megatron/bridge/training/gpt_step.py b/src/megatron/bridge/training/gpt_step.py
index 9744c2827e..694f6c58a3 100644
--- a/src/megatron/bridge/training/gpt_step.py
+++ b/src/megatron/bridge/training/gpt_step.py
@@ -32,6 +32,7 @@
 from megatron.bridge.training.losses import masked_next_token_loss
 from megatron.bridge.training.post_training.distillation import loss_func_kd
 from megatron.bridge.training.state import GlobalState
+from megatron.bridge.training.utils.flop_utils import accumulate_flops_metadata
 from megatron.bridge.training.utils.packed_seq_utils import get_packed_seq_params
 from megatron.bridge.training.utils.pg_utils import get_pg_collection
 
@@ -238,6 +239,21 @@ def _forward_step_common(
         ) = get_batch(data_iterator, state.cfg, use_mtp, pg_collection=pg_collection)
     timers("batch-generator").stop()
 
+    # Accumulate FLOPS metadata across micro-batches. For offline-packed THD
+    # SFT, ``cu_seqlens`` (and ``cu_seqlens_unpadded`` when ``pad_seq_to_mult
+    # > 1``) describe the real sub-sequence boundaries within the pack, so
+    # the helper computes the THD-correct Σᵢ sᵢ² for the attention term
+    # instead of the pack-length² BSHD approximation. train.py resets these
+    # before each step and reads accumulated values afterwards.
+    accumulate_flops_metadata(
+        state,
+        tokens,
+        cu_seqlens=cu_seqlens,
+        cu_seqlens_argmin=cu_seqlens_argmin,
+        cu_seqlens_unpadded=cu_seqlens_unpadded,
+        cu_seqlens_unpadded_argmin=cu_seqlens_unpadded_argmin,
+    )
+
     forward_args = {
         "input_ids": tokens,
         "position_ids": position_ids,
diff --git a/src/megatron/bridge/training/train.py b/src/megatron/bridge/training/train.py
index 41b246e6e4..b61e637dcb 100644
--- a/src/megatron/bridge/training/train.py
+++ b/src/megatron/bridge/training/train.py
@@ -291,8 +291,16 @@ def train(
                 break
     # Track train step elapsed time for throughput logging
     history_wct = None
+    history_real_tokens = None
+    history_packed_tokens = None
+    # Compact THD can distinguish useful tokens from physical packed tokens.
+    # Keep those histories separate from the legacy nominal tokens/sec counter.
+    total_real_tokens_processed = 0
+    total_packed_tokens_processed = 0
     if config.logger.log_throughput_to_tensorboard:
         history_wct = deque(maxlen=config.logger.throughput_window_size + 1)
+        history_real_tokens = deque(maxlen=config.logger.throughput_window_size + 1)
+        history_packed_tokens = deque(maxlen=config.logger.throughput_window_size + 1)
 
     # Wrap forward_backward_func for Full iteration CUDA graph
     forward_backward_func = get_forward_backward_func(
@@ -436,6 +444,10 @@ def train(
         global_state._flops_seqlen_sum = 0
         global_state._flops_seqlen_sq_sum = 0
         global_state._flops_vision_patches = 0
+        # Optional per-forward-step counters. Step functions that do not know
+        # real/packed token counts simply leave these at zero.
+        global_state._throughput_real_tokens = 0
+        global_state._throughput_packed_tokens = 0
 
         (
             loss_dict,
@@ -551,6 +563,8 @@ def train(
         local_seqlen_sum = getattr(global_state, "_flops_seqlen_sum", 0)
         local_seqlen_sq_sum = getattr(global_state, "_flops_seqlen_sq_sum", 0)
         num_vision_patches = getattr(global_state, "_flops_vision_patches", 0)
+        local_real_tokens = getattr(global_state, "_throughput_real_tokens", 0)
+        local_packed_tokens = getattr(global_state, "_throughput_packed_tokens", 0)
         # Coerce to int — getattr on MagicMock test doubles returns a MagicMock
         # (not the default), which breaks the numeric comparisons below.
         if not isinstance(local_seqlen_sum, int):
@@ -559,6 +573,10 @@ def train(
             local_seqlen_sq_sum = 0
         if not isinstance(num_vision_patches, int):
             num_vision_patches = 0
+        if not isinstance(local_real_tokens, int):
+            local_real_tokens = 0
+        if not isinstance(local_packed_tokens, int):
+            local_packed_tokens = 0
 
         # Correct for VPP over-counting: each microbatch's seqlen is accumulated
         # once per virtual stage, but FLOPS formula already covers all stages.
@@ -567,6 +585,8 @@ def train(
             local_seqlen_sum = local_seqlen_sum // vp_size
             local_seqlen_sq_sum = local_seqlen_sq_sum // vp_size
             num_vision_patches = num_vision_patches // vp_size
+            local_real_tokens = local_real_tokens // vp_size
+            local_packed_tokens = local_packed_tokens // vp_size
 
         if local_seqlen_sum > 0:
             seqlen_sum = local_seqlen_sum * dp_size
@@ -589,6 +609,13 @@ def train(
         global_state.train_state.floating_point_operations_so_far += num_floating_point_operations_in_batch
         num_floating_point_operations_so_far = global_state.train_state.floating_point_operations_so_far
         num_floating_point_operations_since_last_log_event += num_floating_point_operations_in_batch
+        if config.logger.log_throughput_to_tensorboard:
+            # These are rank-local counts scaled by DP, mirroring FLOPS metadata.
+            # They intentionally do not replace the existing nominal token metric.
+            total_real_tokens_processed += local_real_tokens * dp_size
+            total_packed_tokens_processed += local_packed_tokens * dp_size
+            history_real_tokens.append(total_real_tokens_processed)
+            history_packed_tokens.append(total_packed_tokens_processed)
 
         # Logging.
         if not config.logger.skip_train_metrics_log:
@@ -627,6 +654,8 @@ def train(
                 history_wct,
                 model,
                 log_max_attention_logit,
+                history_real_tokens=history_real_tokens,
+                history_packed_tokens=history_packed_tokens,
                 loaded_iteration=start_iteration,
                 seq_length=seqlen_sum // batch_size if seqlen_sum else None,
             )
diff --git a/src/megatron/bridge/training/utils/flop_utils.py b/src/megatron/bridge/training/utils/flop_utils.py
index 8ff3ff242c..ad7fd02d3f 100644
--- a/src/megatron/bridge/training/utils/flop_utils.py
+++ b/src/megatron/bridge/training/utils/flop_utils.py
@@ -15,6 +15,7 @@
 import importlib
 from pathlib import Path
 
+import torch
 import torch.nn.functional as F
 
 from megatron.bridge.data.datasets.packing_utils import calculate_avg_seqlen
@@ -26,6 +27,109 @@
 _lora_seq_stats_cache: dict = {}
 
 
+def _real_subseq_lengths(
+    cu_seqlens: torch.Tensor | None,
+    cu_seqlens_argmin: torch.Tensor | None = None,
+    cu_seqlens_unpadded: torch.Tensor | None = None,
+    cu_seqlens_unpadded_argmin: torch.Tensor | None = None,
+) -> torch.Tensor | None:
+    """Extract real (non-pad) sub-sequence lengths from cu_seqlens metadata.
+
+    Prefers ``cu_seqlens_unpadded`` (true sub-sequence boundaries when
+    ``pad_seq_to_mult > 1``) over the padded ``cu_seqlens``. Truncates by the
+    corresponding ``*_argmin`` when provided. Returns ``None`` when no
+    cu_seqlens info is available.
+    """
+    if cu_seqlens_unpadded is not None:
+        cu = cu_seqlens_unpadded.squeeze()
+        argmin = cu_seqlens_unpadded_argmin
+    elif cu_seqlens is not None:
+        cu = cu_seqlens.squeeze()
+        argmin = cu_seqlens_argmin
+    else:
+        return None
+
+    if argmin is not None:
+        cu = cu[: int(argmin.item())]
+
+    if cu.numel() < 2:
+        return cu.new_empty(0, dtype=torch.long)
+
+    sub_seq_lens = (cu[1:] - cu[:-1]).long()
+    return sub_seq_lens[sub_seq_lens > 0]
+
+
+def accumulate_flops_metadata(
+    state,
+    tokens: torch.Tensor | None,
+    *,
+    cu_seqlens: torch.Tensor | None = None,
+    cu_seqlens_argmin: torch.Tensor | None = None,
+    cu_seqlens_unpadded: torch.Tensor | None = None,
+    cu_seqlens_unpadded_argmin: torch.Tensor | None = None,
+    image_grid_thw: torch.Tensor | None = None,
+    video_grid_thw: torch.Tensor | None = None,
+) -> None:
+    """Accumulate per-microbatch FLOPS metadata onto ``state``.
+
+    Writes three accumulators consumed by ``train.py`` at end of step:
+
+    - ``_flops_seqlen_sum``: ``mbs * tokens.shape[1]`` (padded total tokens
+      this microbatch contributes). Drives the linear MLP/proj/logit terms.
+    - ``_flops_seqlen_sq_sum``: Σᵢ sᵢ² over real sub-sequence lengths derived
+      from ``cu_seqlens`` when available (THD-correct attention work), else
+      ``mbs * seq_len²`` (BSHD fallback, matches legacy behavior).
+    - ``_flops_vision_patches``: Σ patches across the provided image/video
+      grid tensors (each shaped ``[num_images, 3]`` with rows ``(t, h, w)``).
+
+    The BSHD fallback applies when cu_seqlens is not provided (e.g. dense
+    pretraining or non-packed SFT) and reproduces the existing single-pack-as-
+    one-sequence computation.
+
+    For THD packed training (offline packed LLM SFT or VLM in-batch packing),
+    treating the whole pack as one length-``seq_len`` sequence over-counts
+    attention FLOPS by a large factor: actual attention work is Σᵢ sᵢ²,
+    not (Σᵢ sᵢ)². Using ``cu_seqlens`` here closes that gap.
+    """
+    if tokens is None:
+        return
+
+    mbs = tokens.shape[0]
+    seq_len = tokens.shape[1]
+    state._flops_seqlen_sum = getattr(state, "_flops_seqlen_sum", 0) + mbs * seq_len
+
+    sub_seq_lens = _real_subseq_lengths(cu_seqlens, cu_seqlens_argmin, cu_seqlens_unpadded, cu_seqlens_unpadded_argmin)
+    if sub_seq_lens is not None and sub_seq_lens.numel() > 0:
+        sq_delta = int((sub_seq_lens.long() ** 2).sum().item())
+    else:
+        sq_delta = mbs * seq_len**2
+    state._flops_seqlen_sq_sum = getattr(state, "_flops_seqlen_sq_sum", 0) + sq_delta
+
+    for grid in (image_grid_thw, video_grid_thw):
+        if grid is not None and grid.numel() > 0:
+            state._flops_vision_patches = getattr(state, "_flops_vision_patches", 0) + int(
+                grid.prod(dim=-1).sum().item()
+            )
+
+
+def accumulate_token_throughput_metadata(
+    state,
+    *,
+    real_tokens: int | None = None,
+    packed_tokens: int | None = None,
+) -> None:
+    """Accumulate optional token-count metadata for compact-packing throughput logging.
+
+    Existing ``throughput/tokens_per_sec`` remains the legacy nominal counter.
+    These optional counters let compact THD report useful non-padding tokens
+    separately from physical packed tokens without changing legacy metrics.
+    """
+    if real_tokens is not None:
+        state._throughput_real_tokens = getattr(state, "_throughput_real_tokens", 0) + int(real_tokens)
+    if packed_tokens is not None:
+        state._throughput_packed_tokens = getattr(state, "_throughput_packed_tokens", 0) + int(packed_tokens)
+
+
 def vit_flops(
     cfg: ConfigContainer,
     batch_size: int,
diff --git a/src/megatron/bridge/training/utils/train_utils.py b/src/megatron/bridge/training/utils/train_utils.py
index 05e8956ef7..7f857763dd 100644
--- a/src/megatron/bridge/training/utils/train_utils.py
+++ b/src/megatron/bridge/training/utils/train_utils.py
@@ -532,6 +532,8 @@ def training_log(
     model: list[MegatronModule],
     pg_collection: Optional[Any] = None,
     log_max_attention_logit: Optional[float] = None,
+    history_real_tokens: Optional[list] = None,
+    history_packed_tokens: Optional[list] = None,
     loaded_iteration: int = 0,
     seq_length: Optional[int] = None,
 ) -> bool:
@@ -685,6 +687,8 @@ def training_log(
                 train_config=train_config,
                 seq_length=config.dataset.seq_length,
                 history_wct=history_wct,
+                history_real_tokens=history_real_tokens,
+                history_packed_tokens=history_packed_tokens,
                 window_size=logger_config.throughput_window_size,
             )
             if writer:
@@ -1265,6 +1269,8 @@ def report_throughput(
     seq_length: int,
     history_wct: list,
     window_size: int,
+    history_real_tokens: Optional[list] = None,
+    history_packed_tokens: Optional[list] = None,
 ) -> dict:
     """
     Logs the training throughput and utilization.
@@ -1338,6 +1344,20 @@ def report_throughput(
             dev_tokens_per_sec = tokens_per_sec / world_size
             metrics.update({"throughput/tokens_per_sec": tokens_per_sec})
             metrics.update({"throughput/device/tokens_per_sec": dev_tokens_per_sec})
+        # Optional compact-packing counters. They are additive metrics, not a
+        # semantic change to the legacy nominal tokens/sec keys above.
+        if history_real_tokens is not None and len(history_real_tokens) >= window_size:
+            elapsed_real_tokens = int(history_real_tokens[-1]) - int(history_real_tokens[0])
+            if elapsed_real_tokens > 0:
+                real_tokens_per_sec = elapsed_real_tokens / elapsed_wct
+                metrics.update({"throughput/real_tokens_per_sec": real_tokens_per_sec})
+                metrics.update({"throughput/device/real_tokens_per_sec": real_tokens_per_sec / world_size})
+        if history_packed_tokens is not None and len(history_packed_tokens) >= window_size:
+            elapsed_packed_tokens = int(history_packed_tokens[-1]) - int(history_packed_tokens[0])
+            if elapsed_packed_tokens > 0:
+                packed_tokens_per_sec = elapsed_packed_tokens / elapsed_wct
+                metrics.update({"throughput/packed_tokens_per_sec": packed_tokens_per_sec})
+                metrics.update({"throughput/device/packed_tokens_per_sec": packed_tokens_per_sec / world_size})
 
         return metrics
 
diff --git a/src/megatron/bridge/training/vlm_step.py b/src/megatron/bridge/training/vlm_step.py
index ad5ca5f2a8..76ea177911 100644
--- a/src/megatron/bridge/training/vlm_step.py
+++ b/src/megatron/bridge/training/vlm_step.py
@@ -27,6 +27,7 @@
     create_masked_next_token_loss_function as _create_loss_function,
 )
 from megatron.bridge.training.state import GlobalState
+from megatron.bridge.training.utils.flop_utils import accumulate_flops_metadata
 from megatron.bridge.training.utils.packed_seq_utils import get_packed_seq_params
 from megatron.bridge.training.utils.padding_utils import (
     pad_or_truncate_2d_to_len,
@@ -435,21 +436,19 @@ def forward_step(
         ) = get_batch(data_iterator, state.cfg, use_mtp, pg_collection=pg_collection)
     timers("batch-generator").stop()
 
-    # Accumulate FLOPS metadata across micro-batches.
-    # Each micro-batch contributes its actual padded seq_length (not cfg.model.seq_length).
-    # train.py resets these before each step and reads accumulated values afterwards.
-    if tokens is not None:
-        mbs = tokens.shape[0]
-        seq_len = tokens.shape[1]
-        state._flops_seqlen_sum = getattr(state, "_flops_seqlen_sum", 0) + mbs * seq_len
-        state._flops_seqlen_sq_sum = getattr(state, "_flops_seqlen_sq_sum", 0) + mbs * seq_len**2
-    if visual_inputs is not None:
-        for attr in ("image_grid_thw", "video_grid_thw"):
-            grid = getattr(visual_inputs, attr, None)
-            if grid is not None and grid.numel() > 0:
-                state._flops_vision_patches = getattr(state, "_flops_vision_patches", 0) + int(
-                    grid.prod(dim=-1).sum().item()
-                )
+    # Accumulate FLOPS metadata across micro-batches. For in-batch packed
+    # batches, ``cu_seqlens`` describes the real sub-sequence boundaries (no
+    # padding sub-seqs), so the helper computes the THD-correct Σᵢ sᵢ² for
+    # the attention term instead of the pack-length² BSHD approximation.
+    # train.py resets these before each step and reads accumulated values
+    # afterwards.
+    accumulate_flops_metadata(
+        state,
+        tokens,
+        cu_seqlens=cu_seqlens,
+        image_grid_thw=getattr(visual_inputs, "image_grid_thw", None) if visual_inputs is not None else None,
+        video_grid_thw=getattr(visual_inputs, "video_grid_thw", None) if visual_inputs is not None else None,
+    )
 
     forward_args = {
         "input_ids": tokens,
diff --git a/tests/unit_tests/training/utils/test_flop_utils.py b/tests/unit_tests/training/utils/test_flop_utils.py
index b0a0814670..4e3e3eed31 100644
--- a/tests/unit_tests/training/utils/test_flop_utils.py
+++ b/tests/unit_tests/training/utils/test_flop_utils.py
@@ -19,8 +19,13 @@
 from unittest.mock import MagicMock, patch
 
 import pytest
+import torch
 
-from megatron.bridge.training.utils.flop_utils import num_floating_point_operations, vit_flops
+from megatron.bridge.training.utils.flop_utils import (
+    accumulate_flops_metadata,
+    num_floating_point_operations,
+    vit_flops,
+)
 
 
 @dataclass
@@ -1662,3 +1667,114 @@ def custom(batch_size):
         assert num_floating_point_operations(cfg, batch_size=4) == sentinel * 4
         # Override must have been invoked twice with the right batch_size args.
         assert captured == [1, 4], f"Override call log mismatch: {captured}"
+
+
+class _State:
+    """Minimal stand-in for GlobalState — just an attribute bag."""
+
+
+class TestAccumulateFlopsMetadata:
+    """Unit tests for ``accumulate_flops_metadata``."""
+
+    def test_bshd_no_cu_seqlens_uses_pack_length_squared(self):
+        # Without cu_seqlens, the accumulator falls back to BSHD math —
+        # mbs * seq_len² — matching the pre-existing behavior on dense
+        # pretraining / non-packed paths.
+        state = _State()
+        tokens = torch.zeros(2, 512)
+        accumulate_flops_metadata(state, tokens)
+        assert state._flops_seqlen_sum == 2 * 512
+        assert state._flops_seqlen_sq_sum == 2 * 512**2
+
+    def test_thd_cu_seqlens_uses_sum_of_squares(self):
+        # cu_seqlens = [0, 256, 512, 4096] → sub-seq lengths [256, 256, 3584].
+        # THD attention work = 256² + 256² + 3584² = 12,975,488; the BSHD
+        # approximation (1 × 4096²) would be 16,777,216 — much larger.
+        state = _State()
+        tokens = torch.zeros(1, 4096)
+        cu_seqlens = torch.tensor([0, 256, 512, 4096])
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_seqlens)
+        assert state._flops_seqlen_sum == 1 * 4096
+        assert state._flops_seqlen_sq_sum == 256**2 + 256**2 + 3584**2
+
+    def test_thd_padded_cu_seqlens_with_argmin(self):
+        # Offline packed SFT pads cu_seqlens for CUDA graphs; the real
+        # entries end at cu_seqlens_argmin. Pad entries past argmin must be
+        # ignored (here they would otherwise contribute zero-length, but we
+        # exercise the truncation explicitly).
+        state = _State()
+        tokens = torch.zeros(1, 8192)
+        cu_seqlens = torch.tensor([0, 1024, 4096, 8192, 8192, 8192, 8192])
+        argmin = torch.tensor(4)  # real entries [0, 1024, 4096, 8192]
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_seqlens, cu_seqlens_argmin=argmin)
+        assert state._flops_seqlen_sq_sum == 1024**2 + 3072**2 + 4096**2
+
+    def test_thd_unpadded_takes_precedence_over_padded(self):
+        # When both cu_seqlens_unpadded and cu_seqlens are present, the
+        # unpadded variant describes the actual sub-sequence boundaries used
+        # by the attention kernel (cu_seqlens_q in PackedSeqParams) and must
+        # be the source of Σᵢ sᵢ².
+        state = _State()
+        tokens = torch.zeros(1, 4096)
+        cu_seqlens_padded = torch.tensor([0, 4096, 4096, 4096])  # 1 pad-aligned sub-seq
+        cu_seqlens_unpadded = torch.tensor([0, 1000, 3500, 4096])  # 3 real sub-seqs
+        accumulate_flops_metadata(
+            state,
+            tokens,
+            cu_seqlens=cu_seqlens_padded,
+            cu_seqlens_unpadded=cu_seqlens_unpadded,
+        )
+        assert state._flops_seqlen_sq_sum == 1000**2 + 2500**2 + 596**2
+
+    def test_accumulates_additively_across_microbatches(self):
+        # Each call adds to existing accumulators (microbatch loop semantics).
+        state = _State()
+        tokens = torch.zeros(1, 128)
+        cu_a = torch.tensor([0, 32, 128])
+        cu_b = torch.tensor([0, 64, 128])
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_a)
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_b)
+        assert state._flops_seqlen_sum == 2 * 128
+        assert state._flops_seqlen_sq_sum == (32**2 + 96**2) + (64**2 + 64**2)
+
+    def test_tokens_none_is_noop(self):
+        state = _State()
+        accumulate_flops_metadata(state, None)
+        assert not hasattr(state, "_flops_seqlen_sum")
+        assert not hasattr(state, "_flops_seqlen_sq_sum")
+
+    def test_visual_inputs_image_and_video_grids(self):
+        state = _State()
+        tokens = torch.zeros(1, 64)
+        accumulate_flops_metadata(
+            state,
+            tokens,
+            image_grid_thw=torch.tensor([[1, 4, 4], [1, 2, 8]]),  # 16 + 16 = 32 patches
+            video_grid_thw=torch.tensor([[2, 2, 2]]),  # 8 patches
+        )
+        assert state._flops_vision_patches == 32 + 8
+
+    def test_empty_cu_seqlens_falls_back_to_bshd(self):
+        # Degenerate cu_seqlens (only one element after argmin truncation)
+        # yields no sub-seqs, so the helper must fall back to BSHD rather
+        # than report 0 attention work.
+        state = _State()
+        tokens = torch.zeros(1, 256)
+        cu_seqlens = torch.tensor([0])
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_seqlens)
+        assert state._flops_seqlen_sq_sum == 1 * 256**2
+
+    def test_thd_substantially_smaller_than_bshd_for_short_samples(self):
+        # Regression check on the headline claim: a pack containing many
+        # short samples has dramatically less attention work than the BSHD
+        # approximation would suggest.
+        state = _State()
+        tokens = torch.zeros(1, 8192)
+        # 32 sub-seqs of length 256 → pack length 8192.
+        cu_seqlens = torch.tensor([i * 256 for i in range(33)])
+        accumulate_flops_metadata(state, tokens, cu_seqlens=cu_seqlens)
+        thd_sq = state._flops_seqlen_sq_sum
+        bshd_sq = 1 * 8192**2
+        # 32 * 256² = 2,097,152 vs 8192² = 67,108,864 → 32× smaller.
+        assert thd_sq == 32 * 256**2
+        assert bshd_sq // thd_sq == 32