Skip to content

bindings: move swap from Fortran to Python #1115

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
havogt wants to merge 6 commits into from
Closed

bindings: move swap from Fortran to Python #1115

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
dc52a12
bindings: move swap from Fortran to Python
havogt Mar 18, 2026
073f425
Apply suggestions from code review
havogt Mar 18, 2026
37c21b2
update references
havogt Mar 18, 2026
491344c
fix test
havogt Mar 20, 2026
1a54450
Merge branch 'main' into blueline_swap_rbf_vec_coeff_e_and_v
havogt Mar 20, 2026
cd4b75c
fix tests
havogt Mar 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 17 additions & 4 deletions model/testing/src/icon4py/model/testing/serialbox.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
# SPDX-License-Identifier: BSD-3-Clause
import functools
import logging
from collections.abc import Sequence
from typing import Final, Literal, TypeAlias

import gt4py.next as gtx
Expand Down Expand Up @@ -84,8 +85,10 @@ def wrapper(self, *args, **kwargs):
def log_meta_info(self):
self.log.info(self.savepoint.metainfo)

def _get_field(self, name, *dimensions, dtype=float):
def _get_field(self, name, *dimensions, dtype=float, transpose: None | Sequence[int] = None):
buffer = np.squeeze(self.serializer.read(name, self.savepoint).astype(dtype))
if transpose is not None:
buffer = np.transpose(buffer, axes=transpose)
buffer = self._reduce_to_dim_size(buffer, dimensions)

self.log.debug(f"{name} {buffer.shape}")
Expand Down Expand Up @@ -673,7 +676,7 @@ def pos_on_tplane_e_y(self):
return self._get_field("pos_on_tplane_e_y", dims.EdgeDim, dims.E2CDim)[:, 0:2]

def rbf_vec_coeff_e(self):
return self._get_field("rbf_vec_coeff_e", dims.EdgeDim, dims.E2C2EDim)
return self._get_field("rbf_vec_coeff_e", dims.EdgeDim, dims.E2C2EDim, transpose=(1, 0))

@IconSavepoint.optionally_registered()
def rbf_vec_coeff_c1(self):
Expand All @@ -694,10 +697,20 @@ def rbf_vec_coeff_c2(self):
return gtx.as_field(dimensions, buffer, allocator=self.backend)

def rbf_vec_coeff_v1(self):
return self._get_field("rbf_vec_coeff_v1", dims.VertexDim, dims.V2EDim)
dimensions = (dims.VertexDim, dims.V2EDim)
buffer = np.squeeze(
self.serializer.read("rbf_vec_coeff_v", self.savepoint).astype(float)[:, 0, :]
).transpose()
buffer = self._reduce_to_dim_size(buffer, dimensions)
return gtx.as_field(dimensions, buffer, allocator=self.backend)

def rbf_vec_coeff_v2(self):
return self._get_field("rbf_vec_coeff_v2", dims.VertexDim, dims.V2EDim)
dimensions = (dims.VertexDim, dims.V2EDim)
buffer = np.squeeze(
self.serializer.read("rbf_vec_coeff_v", self.savepoint).astype(float)[:, 1, :]
).transpose()
buffer = self._reduce_to_dim_size(buffer, dimensions)
return gtx.as_field(dimensions, buffer, allocator=self.backend)

def rbf_vec_idx_v(self):
return self._get_field("rbf_vec_idx_v", dims.VertexDim, dims.V2EDim)
Expand Down
10 changes: 10 additions & 0 deletions tools/src/icon4py/tools/py2fgen/wrappers/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,16 @@
),
]

Float64Array3D: TypeAlias = Annotated[
data_alloc.NDArray,
py2fgen.ArrayParamDescriptor(
rank=3,
dtype=ts.ScalarKind.FLOAT64,
memory_space=py2fgen.MemorySpace.MAYBE_DEVICE,
is_optional=False,
),
]

OptionalFloat64Array1D: TypeAlias = Annotated[
data_alloc.NDArray | None,
py2fgen.ArrayParamDescriptor(
Expand Down
16 changes: 12 additions & 4 deletions tools/src/icon4py/tools/py2fgen/wrappers/diffusion_wrapper.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,7 @@ def diffusion_init(
geofac_grg_y: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
geofac_n2s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
nudgecoeff_e: fa.EdgeField[wpfloat],
rbf_coeff_1: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
rbf_coeff_2: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
rbf_vec_coeff_v: wrapper_common.Float64Array3D,
zd_cellidx: wrapper_common.OptionalInt32Array2D,
zd_vertidx: wrapper_common.OptionalInt32Array2D,
zd_intcoef: wrapper_common.OptionalFloat64Array2D,
Expand Down Expand Up @@ -92,7 +91,8 @@ def diffusion_init(
"Need to initialise grid using 'grid_init' before running 'diffusion_init'."
)

on_gpu = theta_ref_mc.array_ns != np # TODO(havogt): expose `on_gpu` from py2fgen
xp = theta_ref_mc.array_ns
on_gpu = xp != np # TODO(havogt): expose `on_gpu` from py2fgen
actual_backend = wrapper_common.select_backend(
wrapper_common.BackendIntEnum(backend), on_gpu=on_gpu
)
Expand Down Expand Up @@ -134,7 +134,6 @@ def diffusion_init(
dims.KDim: nlev,
}
)
xp = wgtfac_c.array_ns

if zd_cellidx is None:
# then zdiffu_t is False or the list on that rank is empty, then all of the following are not initialized
Expand Down Expand Up @@ -194,6 +193,15 @@ def diffusion_init(
zd_diffcoef=zd_diffcoef,
)

# Create separate fields for the two components of the RBF vector coefficients.
# TODO(havogt): we could use GT4Py's named collections.
rbf_coeff_1 = gtx.as_field(
Copy link
Copy Markdown
Contributor

@muellch muellch Mar 18, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess we should rename those from rbf_coeff_1 to rbf_vec_coeff_v_1 or rbf_coeff_v_1 so it is clear it is the vertex coefficients. There are also edge and cell coeffs ...

Copy link
Copy Markdown
Contributor

@muellch muellch Mar 18, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could also do this later, now these names are icon4py internal?

[dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 0, :]), allocator=allocator
)
rbf_coeff_2 = gtx.as_field(
[dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 1, :]), allocator=allocator
)

# Interpolation state
interpolation_state = DiffusionInterpolationState(
e_bln_c_s=e_bln_c_s,
Expand Down
23 changes: 18 additions & 5 deletions tools/src/icon4py/tools/py2fgen/wrappers/dycore_wrapper.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,10 +56,9 @@ def solve_nh_init(
geofac_rot: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
pos_on_tplane_e_1: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim], gtx.float64],
pos_on_tplane_e_2: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2CDim], gtx.float64],
rbf_vec_coeff_e: gtx.Field[gtx.Dims[dims.EdgeDim, dims.E2C2EDim], gtx.float64],
rbf_vec_coeff_e: wrapper_common.Float64Array2D,
e_bln_c_s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2EDim], gtx.float64],
rbf_coeff_1: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
rbf_coeff_2: gtx.Field[gtx.Dims[dims.VertexDim, dims.V2EDim], gtx.float64],
rbf_vec_coeff_v: wrapper_common.Float64Array3D,
geofac_div: gtx.Field[gtx.Dims[dims.CellDim, dims.C2EDim], gtx.float64],
geofac_n2s: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
geofac_grg_x: gtx.Field[gtx.Dims[dims.CellDim, dims.C2E2CODim], gtx.float64],
Expand Down Expand Up @@ -129,7 +128,8 @@ def solve_nh_init(
if grid_wrapper.grid_state is None:
raise Exception("Need to initialise grid using 'grid_init' before running 'solve_nh_init'.")

on_gpu = c_lin_e.array_ns != np # TODO(havogt): expose `on_gpu` from py2fgen
xp = c_lin_e.array_ns
on_gpu = xp != np # TODO(havogt): expose `on_gpu` from py2fgen
actual_backend = wrapper_common.select_backend(
wrapper_common.BackendIntEnum(backend), on_gpu=on_gpu
)
Expand Down Expand Up @@ -181,6 +181,19 @@ def solve_nh_init(
)
nonhydro_params = solve_nonhydro.NonHydrostaticParams(config)

# Create separate fields for the two components of the RBF vector coefficients.
# TODO(havogt): we could use GT4Py's named collections.
rbf_coeff_1 = gtx.as_field(
[dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 0, :]), allocator=allocator
)
rbf_coeff_2 = gtx.as_field(
[dims.VertexDim, dims.V2EDim], xp.transpose(rbf_vec_coeff_v[:, 1, :]), allocator=allocator
)

# Swap indices in rbf_vec_coeff_e. TODO(havogt): Should eventually be done on the Fortran side.
rbf_vec_coeff_e_transposed = gtx.as_field(
[dims.EdgeDim, dims.E2C2EDim], xp.transpose(rbf_vec_coeff_e), allocator=allocator
)
interpolation_state = dycore_states.InterpolationState(
c_lin_e=c_lin_e,
c_intp=c_intp,
Expand All @@ -189,7 +202,7 @@ def solve_nh_init(
geofac_rot=geofac_rot,
pos_on_tplane_e_1=pos_on_tplane_e_1[:, 0:2],
pos_on_tplane_e_2=pos_on_tplane_e_2[:, 0:2],
rbf_vec_coeff_e=rbf_vec_coeff_e,
rbf_vec_coeff_e=rbf_vec_coeff_e_transposed,
e_bln_c_s=e_bln_c_s,
rbf_coeff_1=rbf_coeff_1,
rbf_coeff_2=rbf_coeff_2,
Expand Down
Loading