Resolved conflicts

model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion.py

@@ -68,6 +68,7 @@
     CPD,
     DEFAULT_PHYSICS_DYNAMICS_TIMESTEP_RATIO,
     GAS_CONSTANT_DRY_AIR,
+    dbl_eps,
 )
 from icon4py.model.common.decomposition.definitions import ExchangeRuntime, SingleNodeExchange
 from icon4py.model.common.dimension import CellDim, EdgeDim, KDim, VertexDim
@@ -771,6 +772,7 @@ def _do_diffusion_step(
             theta_v=prognostic_state.theta_v,
             theta_ref_mc=self.metric_state.theta_ref_mc,
             thresh_tdiff=self.thresh_tdiff,
+            smallest_vpfloat=dbl_eps,
             kh_smag_e=self.kh_smag_e,
             horizontal_start=edge_start_nudging,
             horizontal_end=edge_end_halo,

model/atmosphere/diffusion/src/icon4py/model/atmosphere/diffusion/diffusion_utils.py

@@ -16,9 +16,10 @@
 from gt4py.next import as_field
 from gt4py.next.common import Dimension, Field
 from gt4py.next.ffront.decorator import field_operator, program
-from gt4py.next.ffront.fbuiltins import broadcast, int32, maximum, minimum
+from gt4py.next.ffront.fbuiltins import broadcast, int32, minimum
 
-from icon4py.model.common.dimension import CellDim, EdgeDim, KDim, Koff, VertexDim
+from icon4py.model.common.dimension import CellDim, EdgeDim, KDim, VertexDim
+from icon4py.model.common.math.smagorinsky import en_smag_fac_for_zero_nshift
 
 
 # TODO(Magdalena): fix duplication: duplicated from test testutils/utils.py
@@ -98,35 +99,6 @@ def setup_fields_for_initial_step(
     _setup_fields_for_initial_step(k4, hdiff_efdt_ratio, out=(diff_multfac_vn, smag_limit))
 
 
-@field_operator
-def _en_smag_fac_for_zero_nshift(
-    vect_a: Field[[KDim], float],
-    hdiff_smag_fac: float,
-    hdiff_smag_fac2: float,
-    hdiff_smag_fac3: float,
-    hdiff_smag_fac4: float,
-    hdiff_smag_z: float,
-    hdiff_smag_z2: float,
-    hdiff_smag_z3: float,
-    hdiff_smag_z4: float,
-) -> Field[[KDim], float]:
-    dz21 = hdiff_smag_z2 - hdiff_smag_z
-    alin = (hdiff_smag_fac2 - hdiff_smag_fac) / dz21
-    df32 = hdiff_smag_fac3 - hdiff_smag_fac2
-    df42 = hdiff_smag_fac4 - hdiff_smag_fac2
-    dz32 = hdiff_smag_z3 - hdiff_smag_z2
-    dz42 = hdiff_smag_z4 - hdiff_smag_z2
-
-    bqdr = (df42 * dz32 - df32 * dz42) / (dz32 * dz42 * (dz42 - dz32))
-    aqdr = df32 / dz32 - bqdr * dz32
-    zf = 0.5 * (vect_a + vect_a(Koff[1]))
-
-    dzlin = minimum(dz21, maximum(0.0, zf - hdiff_smag_z))
-    dzqdr = minimum(dz42, maximum(0.0, zf - hdiff_smag_z2))
-    enh_smag_fac = hdiff_smag_fac + (dzlin * alin) + dzqdr * (aqdr + dzqdr * bqdr)
-    return enh_smag_fac
-
-
 @field_operator
 def _init_diffusion_local_fields_for_regular_timestemp(
     k4: float,
@@ -143,7 +115,7 @@ def _init_diffusion_local_fields_for_regular_timestemp(
 ) -> tuple[Field[[KDim], float], Field[[KDim], float], Field[[KDim], float]]:
     diff_multfac_vn = _setup_runtime_diff_multfac_vn(k4, dyn_substeps)
     smag_limit = _setup_smag_limit(diff_multfac_vn)
-    enh_smag_fac = _en_smag_fac_for_zero_nshift(
+    enh_smag_fac = en_smag_fac_for_zero_nshift(
         vect_a,
         hdiff_smag_fac,
         hdiff_smag_fac2,

model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion.py

@@ -19,14 +19,13 @@
 from icon4py.model.atmosphere.diffusion.diffusion_utils import scale_k
 from icon4py.model.common.grid.horizontal import CellParams, EdgeParams
 from icon4py.model.common.grid.vertical import VerticalModelParams
+from icon4py.model.common.test_utils.reference_funcs import enhanced_smagorinski_factor_numpy
 from icon4py.model.common.test_utils.serialbox_utils import IconDiffusionInitSavepoint
 
-from .utils import (
-    diff_multfac_vn_numpy,
-    enhanced_smagorinski_factor_numpy,
-    smag_limit_numpy,
-    verify_diffusion_fields,
-)
+from .utils import diff_multfac_vn_numpy, smag_limit_numpy, verify_diffusion_fields
+
+
+backend = icon4py.model.atmosphere.diffusion.diffusion.backend
 
 
 backend = icon4py.model.atmosphere.diffusion.diffusion.backend

model/atmosphere/diffusion/tests/diffusion_tests/test_diffusion_utils.py

@@ -16,7 +16,6 @@
 
 from icon4py.model.atmosphere.diffusion.diffusion import DiffusionParams
 from icon4py.model.atmosphere.diffusion.diffusion_utils import (
-    _en_smag_fac_for_zero_nshift,
     _setup_runtime_diff_multfac_vn,
     _setup_smag_limit,
     scale_k,
@@ -27,7 +26,7 @@
 from icon4py.model.common.grid.simple import SimpleGrid
 from icon4py.model.common.test_utils.helpers import random_field, zero_field
 
-from .utils import diff_multfac_vn_numpy, enhanced_smagorinski_factor_numpy, smag_limit_numpy
+from .utils import diff_multfac_vn_numpy, smag_limit_numpy
 
 
 def initial_diff_multfac_vn_numpy(shape, k4, hdiff_efdt_ratio):
@@ -104,25 +103,10 @@ def test_diff_multfac_vn_smag_limit_for_loop_run_with_k4_substeps(backend):
     assert np.allclose(expected_smag_limit, smag_limit.asnumpy())
 
 
-def test_init_enh_smag_fac(backend):
-    grid = SimpleGrid()
-    enh_smag_fac = zero_field(grid, KDim)
-    a_vec = random_field(grid, KDim, low=1.0, high=10.0, extend={KDim: 1})
-    fac = (0.67, 0.5, 1.3, 0.8)
-    z = (0.1, 0.2, 0.3, 0.4)
-
-    enhanced_smag_fac_np = enhanced_smagorinski_factor_numpy(fac, z, a_vec.asnumpy())
-
-    _en_smag_fac_for_zero_nshift.with_backend(backend)(
-        a_vec, *fac, *z, out=enh_smag_fac, offset_provider={"Koff": KDim}
-    )
-    assert np.allclose(enhanced_smag_fac_np, enh_smag_fac.asnumpy())
-
-
 def test_set_zero_vertex_k(backend):
     grid = SimpleGrid()
     f = random_field(grid, VertexDim, KDim)
-    set_zero_v_k.with_backend(backend)(f, offset_provider={})
+    set_zero_v_k(f, offset_provider={})
     assert np.allclose(0.0, f.asnumpy())
 
 

model/atmosphere/diffusion/tests/diffusion_tests/utils.py

@@ -58,19 +58,3 @@ def smag_limit_numpy(func, *args):
 def diff_multfac_vn_numpy(shape, k4, substeps):
     factor = min(1.0 / 128.0, k4 * substeps / 3.0)
     return factor * np.ones(shape)
-
-
-def enhanced_smagorinski_factor_numpy(factor_in, heigths_in, a_vec: np.ndarray):
-    alin = (factor_in[1] - factor_in[0]) / (heigths_in[1] - heigths_in[0])
-    df32 = factor_in[2] - factor_in[1]
-    df42 = factor_in[3] - factor_in[1]
-    dz32 = heigths_in[2] - heigths_in[1]
-    dz42 = heigths_in[3] - heigths_in[1]
-    bqdr = (df42 * dz32 - df32 * dz42) / (dz32 * dz42 * (dz42 - dz32))
-    aqdr = df32 / dz32 - bqdr * dz32
-    zf = 0.5 * (a_vec[:-1] + a_vec[1:])
-    max0 = np.maximum(0.0, zf - heigths_in[0])
-    dzlin = np.minimum(heigths_in[1] - heigths_in[0], max0)
-    max1 = np.maximum(0.0, zf - heigths_in[1])
-    dzqdr = np.minimum(heigths_in[3] - heigths_in[1], max1)
-    return factor_in[0] + dzlin * alin + dzqdr * (aqdr + dzqdr * bqdr)