From e92896bd5d263a2ddc13811aa8b866e38865cc0c Mon Sep 17 00:00:00 2001 From: Magdalena Luz Date: Tue, 13 Sep 2022 17:19:12 +0200 Subject: [PATCH 01/10] run pre-commit --- advection/setup.py | 1 + advection/src/icon4py/advection/__init__.py | 1 - .../advection/hflx_limiter_pd_stencil_01.py | 54 ++++++++++++------- .../advection/hflx_limiter_pd_stencil_02.py | 24 +++++---- .../advection/rbf_intp_edge_stencil_01.py | 9 ++-- .../advection/step_advection_stencil_04.py | 8 ++- .../tests/test_hflx_limiter_pd_stencil_01.py | 27 ++++++---- .../tests/test_hflx_limiter_pd_stencil_02.py | 18 ++++--- .../tests/test_rbf_intp_edge_stencil_01.py | 20 +++---- .../tests/test_step_advection_stencil_04.py | 23 ++++---- .../icon4py/atm_dyn_iconam/compute_airmass.py | 1 - atm_dyn_iconam/tests/test_compute_airmass.py | 4 +- 12 files changed, 110 insertions(+), 80 deletions(-) diff --git a/advection/setup.py b/advection/setup.py index ceec936ea..9c9f7b81c 100644 --- a/advection/setup.py +++ b/advection/setup.py @@ -1,3 +1,4 @@ +# ICON4Py - ICON inspired code in Python and GT4Py # # Copyright (c) 2022, ETH Zurich and MeteoSwiss # All rights reserved. diff --git a/advection/src/icon4py/advection/__init__.py b/advection/src/icon4py/advection/__init__.py index a6d2d236c..15dfdb009 100644 --- a/advection/src/icon4py/advection/__init__.py +++ b/advection/src/icon4py/advection/__init__.py @@ -10,4 +10,3 @@ # distribution for a copy of the license or check . # # SPDX-License-Identifier: GPL-3.0-or-later - diff --git a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py index 94c13720b..e157efeaf 100644 --- a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py +++ b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py @@ -12,8 +12,16 @@ # SPDX-License-Identifier: GPL-3.0-or-later from functional.ffront.decorator import field_operator, program -from functional.ffront.fbuiltins import Field, maximum, minimum, neighbor_sum, broadcast -from icon4py.common.dimension import CellDim, EdgeDim, KDim, C2EDim, C2E +from functional.ffront.fbuiltins import ( + Field, + broadcast, + maximum, + minimum, + neighbor_sum, +) + +from icon4py.common.dimension import C2E, C2EDim, CellDim, EdgeDim, KDim + @field_operator def _hflx_limiter_pd_stencil_01( @@ -21,20 +29,28 @@ def _hflx_limiter_pd_stencil_01( p_cc: Field[[CellDim, KDim], float], p_rhodz_now: Field[[CellDim, KDim], float], p_mflx_tracer_h: Field[[EdgeDim, KDim], float], - p_dtime:float, + p_dtime: float, dbl_eps: float, ) -> Field[[CellDim, KDim], float]: -# p_m: Field[[CellDim, KDim], float] -# p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) -# p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[1])*geofac_div(1)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[2])*geofac_div(2)*p_dtime) - #p_m = neighbor_sum(a_a, axis=C2EDim) - #r_m = minimum(1.0, (p_cc*p_rhodz_now)/(p_m+dbl_eps)) -# r_m = minimum(1.0, (p_cc*p_rhodz_now)/(1.0+dbl_eps)) + # p_m: Field[[CellDim, KDim], float] + # p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) + # p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[1])*geofac_div(1)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[2])*geofac_div(2)*p_dtime) + # p_m = neighbor_sum(a_a, axis=C2EDim) + # r_m = minimum(1.0, (p_cc*p_rhodz_now)/(p_m+dbl_eps)) + # r_m = minimum(1.0, (p_cc*p_rhodz_now)/(1.0+dbl_eps)) - p_m = neighbor_sum(maximum(broadcast(0.0, (CellDim, KDim)), p_mflx_tracer_h(C2E)*geofac_div*p_dtime), axis=C2EDim) - r_m = minimum(broadcast(1.0, (CellDim, KDim)), (p_cc*p_rhodz_now)/(p_m+dbl_eps)) + p_m = neighbor_sum( + maximum( + broadcast(0.0, (CellDim, KDim)), p_mflx_tracer_h(C2E) * geofac_div * p_dtime + ), + axis=C2EDim, + ) + r_m = minimum( + broadcast(1.0, (CellDim, KDim)), (p_cc * p_rhodz_now) / (p_m + dbl_eps) + ) return r_m + @program def hflx_limiter_pd_stencil_01( geofac_div: Field[[CellDim, C2EDim], float], @@ -42,15 +58,15 @@ def hflx_limiter_pd_stencil_01( p_rhodz_now: Field[[CellDim, KDim], float], p_mflx_tracer_h: Field[[EdgeDim, KDim], float], r_m: Field[[CellDim, KDim], float], - p_dtime:float, + p_dtime: float, dbl_eps: float, ): _hflx_limiter_pd_stencil_01( - geofac_div, - p_cc, - p_rhodz_now, - p_mflx_tracer_h, - p_dtime, - dbl_eps, - out=r_m, + geofac_div, + p_cc, + p_rhodz_now, + p_mflx_tracer_h, + p_dtime, + dbl_eps, + out=r_m, ) diff --git a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_02.py b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_02.py index 2285800f8..ac59f1690 100644 --- a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_02.py +++ b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_02.py @@ -14,26 +14,28 @@ from functional.ffront.decorator import field_operator, program from functional.ffront.fbuiltins import Field, where -from icon4py.common.dimension import CellDim, EdgeDim, KDim, E2C +from icon4py.common.dimension import E2C, CellDim, EdgeDim, KDim + @field_operator def _hflx_limiter_pd_stencil_02( - refin_ctrl: Field[[EdgeDim], float], + refin_ctrl: Field[[EdgeDim], float], r_m: Field[[CellDim, KDim], float], p_mflx_tracer_h_in: Field[[EdgeDim, KDim], float], bound: float, ) -> Field[[EdgeDim, KDim], float]: - p_mflx_tracer_h_out = where ( + p_mflx_tracer_h_out = where( refin_ctrl == bound, p_mflx_tracer_h_in, - where ( - (p_mflx_tracer_h_in > 0.0) | (p_mflx_tracer_h_in==0.0), + where( + (p_mflx_tracer_h_in > 0.0) | (p_mflx_tracer_h_in == 0.0), p_mflx_tracer_h_in * r_m(E2C[0]), p_mflx_tracer_h_in * r_m(E2C[1]), ), ) return p_mflx_tracer_h_out + @program def hflx_limiter_pd_stencil_02( refin_ctrl: Field[[EdgeDim], float], @@ -43,9 +45,9 @@ def hflx_limiter_pd_stencil_02( bound: float, ): _hflx_limiter_pd_stencil_02( - refin_ctrl, - r_m, - p_mflx_tracer_h_in, - bound, - out=p_mflx_tracer_h_out, -) + refin_ctrl, + r_m, + p_mflx_tracer_h_in, + bound, + out=p_mflx_tracer_h_out, + ) diff --git a/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py b/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py index 76a97b58c..b879d031b 100644 --- a/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py +++ b/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py @@ -19,17 +19,18 @@ @field_operator def _rbf_intp_edge_stencil_01( - p_vn_in: Field[[EdgeDim, KDim], float], - # ptr_coeff: Field[[EdgeDim, E2CDim, E2C2EDim], float] + p_vn_in: Field[[EdgeDim, KDim], float], + # ptr_coeff: Field[[EdgeDim, E2CDim, E2C2EDim], float] ptr_coeff: Field[[EdgeDim, E2C2EDim], float], ) -> Field[[EdgeDim, KDim], float]: - p_vt_out = neighbor_sum(p_vn_in(E2C2E) * ptr_coeff, axis=E2C2EDim) + p_vt_out = neighbor_sum(p_vn_in(E2C2E) * ptr_coeff, axis=E2C2EDim) return p_vt_out + @program def rbf_intp_edge_stencil_01( p_vn_in: Field[[EdgeDim, KDim], float], ptr_coeff: Field[[EdgeDim, E2C2EDim], float], p_vt_out: Field[[EdgeDim, KDim], float], ): - _rbf_intp_edge_stencil_01(p_vn_in, ptr_coeff, out=p_vt_out) + _rbf_intp_edge_stencil_01(p_vn_in, ptr_coeff, out=p_vt_out) diff --git a/advection/src/icon4py/advection/step_advection_stencil_04.py b/advection/src/icon4py/advection/step_advection_stencil_04.py index f299d11d9..f0177fabe 100644 --- a/advection/src/icon4py/advection/step_advection_stencil_04.py +++ b/advection/src/icon4py/advection/step_advection_stencil_04.py @@ -16,15 +16,17 @@ from icon4py.common.dimension import CellDim, KDim + @field_operator def _step_advection_stencil_04( p_tracer_now: Field[[CellDim, KDim], float], p_tracer_new: Field[[CellDim, KDim], float], p_dtime: float, ) -> Field[[CellDim, KDim], float]: - opt_ddt_tracer_adv=(p_tracer_new-p_tracer_now)/p_dtime + opt_ddt_tracer_adv = (p_tracer_new - p_tracer_now) / p_dtime return opt_ddt_tracer_adv + @program def step_advection_stencil_04( p_tracer_now: Field[[CellDim, KDim], float], @@ -32,4 +34,6 @@ def step_advection_stencil_04( opt_ddt_tracer_adv: Field[[CellDim, KDim], float], p_dtime: float, ): - _step_advection_stencil_04(p_tracer_now, p_tracer_new, p_dtime, out=opt_ddt_tracer_adv) + _step_advection_stencil_04( + p_tracer_now, p_tracer_new, p_dtime, out=opt_ddt_tracer_adv + ) diff --git a/advection/tests/test_hflx_limiter_pd_stencil_01.py b/advection/tests/test_hflx_limiter_pd_stencil_01.py index a75ffe6bb..6cd13adde 100644 --- a/advection/tests/test_hflx_limiter_pd_stencil_01.py +++ b/advection/tests/test_hflx_limiter_pd_stencil_01.py @@ -11,15 +11,16 @@ # # SPDX-License-Identifier: GPL-3.0-or-later -import numpy as np +import numpy as np from icon4py.advection.hflx_limiter_pd_stencil_01 import ( hflx_limiter_pd_stencil_01, ) -from icon4py.common.dimension import CellDim, EdgeDim, KDim, C2EDim +from icon4py.common.dimension import C2EDim, CellDim, EdgeDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh from icon4py.testutils.utils import random_field, zero_field + def hflx_limiter_pd_stencil_01_numpy( c2e: np.array, geofac_div: np.array, @@ -29,16 +30,24 @@ def hflx_limiter_pd_stencil_01_numpy( p_dtime, dbl_eps, ): -# p_m: np.array - p_mflx_tracer_h_c2e = p_mflx_tracer_h[c2e] - #p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) + max(0.0, geofac_div[:, 1]*p_mflx_tracer_h_c2e[:, 1]*p_dtime) + max(0.0, geofac_div[:, 2]*p_mflx_tracer_h_c2e[:, 2]*p_dtime) - geofac_div=np.expand_dims(geofac_div,axis=-1) -# p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) - p_m = np.sum(np.maximum(np.zeros(p_mflx_tracer_h_c2e.shape,dtype=float), p_mflx_tracer_h_c2e*geofac_div*p_dtime)) - r_m = np.minimum(np.ones(p_cc.shape,dtype=float),(p_cc*p_rhodz_now)/(p_m+dbl_eps)) + # p_m: np.array + p_mflx_tracer_h_c2e = p_mflx_tracer_h[c2e] + # p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) + max(0.0, geofac_div[:, 1]*p_mflx_tracer_h_c2e[:, 1]*p_dtime) + max(0.0, geofac_div[:, 2]*p_mflx_tracer_h_c2e[:, 2]*p_dtime) + geofac_div = np.expand_dims(geofac_div, axis=-1) + # p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) + p_m = np.sum( + np.maximum( + np.zeros(p_mflx_tracer_h_c2e.shape, dtype=float), + p_mflx_tracer_h_c2e * geofac_div * p_dtime, + ) + ) + r_m = np.minimum( + np.ones(p_cc.shape, dtype=float), (p_cc * p_rhodz_now) / (p_m + dbl_eps) + ) return r_m + def test_hflx_limiter_pd_stencil_01(): mesh = SimpleMesh() geofac_div = random_field(mesh, CellDim, C2EDim) diff --git a/advection/tests/test_hflx_limiter_pd_stencil_02.py b/advection/tests/test_hflx_limiter_pd_stencil_02.py index 2079beb60..fdb93587c 100644 --- a/advection/tests/test_hflx_limiter_pd_stencil_02.py +++ b/advection/tests/test_hflx_limiter_pd_stencil_02.py @@ -11,28 +11,29 @@ # # SPDX-License-Identifier: GPL-3.0-or-later -import numpy as np +import numpy as np from icon4py.advection.hflx_limiter_pd_stencil_02 import ( hflx_limiter_pd_stencil_02, ) from icon4py.common.dimension import CellDim, EdgeDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh -from icon4py.testutils.utils import random_field, zero_field +from icon4py.testutils.utils import random_field + def hflx_limiter_pd_stencil_02_numpy( e2c: np.array, - refin_ctrl: np.array, + refin_ctrl: np.array, r_m: np.array, p_mflx_tracer_h_in: np.array, bound, ): r_m_e2c = r_m[e2c] refin_ctrl = np.expand_dims(refin_ctrl, axis=-1) - p_mflx_tracer_h_out = np.where ( + p_mflx_tracer_h_out = np.where( refin_ctrl != bound, - np.where ( - p_mflx_tracer_h_in >=0, + np.where( + p_mflx_tracer_h_in >= 0, p_mflx_tracer_h_in * r_m_e2c[:, 0], p_mflx_tracer_h_in * r_m_e2c[:, 1], ), @@ -40,6 +41,7 @@ def hflx_limiter_pd_stencil_02_numpy( ) return p_mflx_tracer_h_out + def test_hflx_limiter_pd_stencil_02(): mesh = SimpleMesh() @@ -47,7 +49,7 @@ def test_hflx_limiter_pd_stencil_02(): r_m = random_field(mesh, CellDim, KDim) p_mflx_tracer_h_in = random_field(mesh, EdgeDim, KDim) p_mflx_tracer_h_out = random_field(mesh, EdgeDim, KDim) - bound=np.float64(7.0) + bound = np.float64(7.0) ref = hflx_limiter_pd_stencil_02_numpy( mesh.e2c, @@ -58,7 +60,7 @@ def test_hflx_limiter_pd_stencil_02(): ) hflx_limiter_pd_stencil_02( - refin_ctrl, + refin_ctrl, r_m, p_mflx_tracer_h_in, p_mflx_tracer_h_out, diff --git a/advection/tests/test_rbf_intp_edge_stencil_01.py b/advection/tests/test_rbf_intp_edge_stencil_01.py index f9a20110e..96cd338d1 100644 --- a/advection/tests/test_rbf_intp_edge_stencil_01.py +++ b/advection/tests/test_rbf_intp_edge_stencil_01.py @@ -11,34 +11,34 @@ # # SPDX-License-Identifier: GPL-3.0-or-later -import numpy as np +import numpy as np -from icon4py.advection.rbf_intp_edge_stencil_01 import ( - rbf_intp_edge_stencil_01, -) -from icon4py.common.dimension import EdgeDim, E2C2EDim, KDim +from icon4py.advection.rbf_intp_edge_stencil_01 import rbf_intp_edge_stencil_01 +from icon4py.common.dimension import E2C2EDim, EdgeDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh from icon4py.testutils.utils import random_field, zero_field + def rbf_intp_edge_stencil_01_numpy( e2c2e: np.array, - p_vn_in: np.array, + p_vn_in: np.array, ptr_coeff: np.array, ) -> np.array: ptr_coeff = np.expand_dims(ptr_coeff, axis=-1) p_vt_out = np.sum(p_vn_in[e2c2e] * ptr_coeff, axis=1) -# p_vn_ptr = p_vn_in[e2c2e] * ptr_coeff -# p_vt_out = np.sum(p_vn_ptr, axis=1) + # p_vn_ptr = p_vn_in[e2c2e] * ptr_coeff + # p_vt_out = np.sum(p_vn_ptr, axis=1) return p_vt_out + def test_rbf_intp_edge_stencil_01(): - mesh= SimpleMesh() + mesh = SimpleMesh() p_vn_in = random_field(mesh, EdgeDim, KDim) ptr_coeff = random_field(mesh, EdgeDim, E2C2EDim) - p_vt_out = zero_field(mesh, EdgeDim, KDim) + p_vt_out = zero_field(mesh, EdgeDim, KDim) ref = rbf_intp_edge_stencil_01_numpy( mesh.e2c2e, diff --git a/advection/tests/test_step_advection_stencil_04.py b/advection/tests/test_step_advection_stencil_04.py index 463293458..6da3ec9b5 100644 --- a/advection/tests/test_step_advection_stencil_04.py +++ b/advection/tests/test_step_advection_stencil_04.py @@ -13,34 +13,33 @@ import numpy as np -from icon4py.advection.step_advection_stencil_04 import ( - step_advection_stencil_04, -) - +from icon4py.advection.step_advection_stencil_04 import step_advection_stencil_04 from icon4py.common.dimension import CellDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh from icon4py.testutils.utils import random_field, zero_field + def step_advection_stencil_04_numpy( - p_tracer_now: np.array, - p_tracer_new: np.array, + p_tracer_now: np.array, + p_tracer_new: np.array, p_dtime, ) -> np.array: - opt_ddt_tracer_adv = (p_tracer_new-p_tracer_now)/p_dtime + opt_ddt_tracer_adv = (p_tracer_new - p_tracer_now) / p_dtime return opt_ddt_tracer_adv + def test_step_advection_stencil_04(): mesh = SimpleMesh() p_tracer_now = random_field(mesh, CellDim, KDim) p_tracer_new = random_field(mesh, CellDim, KDim) - opt_ddt_tracer_adv = zero_field(mesh, CellDim, KDim) - p_dtime = np.float64(5.0) + opt_ddt_tracer_adv = zero_field(mesh, CellDim, KDim) + p_dtime = np.float64(5.0) ref = step_advection_stencil_04_numpy( np.asarray(p_tracer_now), - np.asarray(p_tracer_new), - p_dtime, + np.asarray(p_tracer_new), + p_dtime, ) step_advection_stencil_04( p_tracer_now, @@ -48,5 +47,5 @@ def test_step_advection_stencil_04(): opt_ddt_tracer_adv, p_dtime, offset_provider={}, -) + ) assert np.allclose(opt_ddt_tracer_adv, ref) diff --git a/atm_dyn_iconam/src/icon4py/atm_dyn_iconam/compute_airmass.py b/atm_dyn_iconam/src/icon4py/atm_dyn_iconam/compute_airmass.py index c259bea68..bac9061ce 100644 --- a/atm_dyn_iconam/src/icon4py/atm_dyn_iconam/compute_airmass.py +++ b/atm_dyn_iconam/src/icon4py/atm_dyn_iconam/compute_airmass.py @@ -26,7 +26,6 @@ def _compute_airmass( return rho_in * ddqz_z_full_in * deepatmo_t1mc_in - @program def compute_airmass( rho_in: Field[[CellDim, KDim], float], diff --git a/atm_dyn_iconam/tests/test_compute_airmass.py b/atm_dyn_iconam/tests/test_compute_airmass.py index 0ca950c89..881838f8b 100644 --- a/atm_dyn_iconam/tests/test_compute_airmass.py +++ b/atm_dyn_iconam/tests/test_compute_airmass.py @@ -13,9 +13,7 @@ import numpy as np -from icon4py.atm_dyn_iconam.compute_airmass import ( - compute_airmass, -) +from icon4py.atm_dyn_iconam.compute_airmass import compute_airmass from icon4py.common.dimension import CellDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh from icon4py.testutils.utils import random_field From e862f05db4a7dbdab47b806ca25b346b2243f369 Mon Sep 17 00:00:00 2001 From: Magdalena Luz Date: Tue, 13 Sep 2022 21:56:44 +0200 Subject: [PATCH 02/10] pre-commit: remove commented out code --- advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py | 1 - advection/tests/test_rbf_intp_edge_stencil_01.py | 3 --- 2 files changed, 4 deletions(-) diff --git a/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py b/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py index b879d031b..46e106942 100644 --- a/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py +++ b/advection/src/icon4py/advection/rbf_intp_edge_stencil_01.py @@ -20,7 +20,6 @@ @field_operator def _rbf_intp_edge_stencil_01( p_vn_in: Field[[EdgeDim, KDim], float], - # ptr_coeff: Field[[EdgeDim, E2CDim, E2C2EDim], float] ptr_coeff: Field[[EdgeDim, E2C2EDim], float], ) -> Field[[EdgeDim, KDim], float]: p_vt_out = neighbor_sum(p_vn_in(E2C2E) * ptr_coeff, axis=E2C2EDim) diff --git a/advection/tests/test_rbf_intp_edge_stencil_01.py b/advection/tests/test_rbf_intp_edge_stencil_01.py index 96cd338d1..84bd29aa4 100644 --- a/advection/tests/test_rbf_intp_edge_stencil_01.py +++ b/advection/tests/test_rbf_intp_edge_stencil_01.py @@ -27,9 +27,6 @@ def rbf_intp_edge_stencil_01_numpy( ptr_coeff = np.expand_dims(ptr_coeff, axis=-1) p_vt_out = np.sum(p_vn_in[e2c2e] * ptr_coeff, axis=1) - # p_vn_ptr = p_vn_in[e2c2e] * ptr_coeff - # p_vt_out = np.sum(p_vn_ptr, axis=1) - return p_vt_out From 916a1f7e6c76b97ccaa9ec132ffc1a1d8acf9ead Mon Sep 17 00:00:00 2001 From: Magdalena Luz Date: Tue, 13 Sep 2022 21:57:57 +0200 Subject: [PATCH 03/10] fix hflx_limiter_pd_stencil_01 --- .../advection/hflx_limiter_pd_stencil_01.py | 32 ++++++------------- .../tests/test_hflx_limiter_pd_stencil_01.py | 27 +++++++--------- 2 files changed, 21 insertions(+), 38 deletions(-) diff --git a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py index e157efeaf..55d84ef3e 100644 --- a/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py +++ b/advection/src/icon4py/advection/hflx_limiter_pd_stencil_01.py @@ -12,48 +12,36 @@ # SPDX-License-Identifier: GPL-3.0-or-later from functional.ffront.decorator import field_operator, program -from functional.ffront.fbuiltins import ( - Field, - broadcast, - maximum, - minimum, - neighbor_sum, -) +from functional.ffront.fbuiltins import Field, broadcast, maximum, minimum -from icon4py.common.dimension import C2E, C2EDim, CellDim, EdgeDim, KDim +from icon4py.common.dimension import C2CE, C2E, CEDim, CellDim, EdgeDim, KDim @field_operator def _hflx_limiter_pd_stencil_01( - geofac_div: Field[[CellDim, C2EDim], float], + geofac_div: Field[[CEDim], float], p_cc: Field[[CellDim, KDim], float], p_rhodz_now: Field[[CellDim, KDim], float], p_mflx_tracer_h: Field[[EdgeDim, KDim], float], p_dtime: float, dbl_eps: float, ) -> Field[[CellDim, KDim], float]: - # p_m: Field[[CellDim, KDim], float] - # p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) - # p_m = maximum(0.0, p_mflx_tracer_h(C2E[0])*geofac_div(0)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[1])*geofac_div(1)*p_dtime) + maximum(0.0, p_mflx_tracer_h(C2E[2])*geofac_div(2)*p_dtime) - # p_m = neighbor_sum(a_a, axis=C2EDim) - # r_m = minimum(1.0, (p_cc*p_rhodz_now)/(p_m+dbl_eps)) - # r_m = minimum(1.0, (p_cc*p_rhodz_now)/(1.0+dbl_eps)) + zero = broadcast(0.0, (CellDim, KDim)) - p_m = neighbor_sum( - maximum( - broadcast(0.0, (CellDim, KDim)), p_mflx_tracer_h(C2E) * geofac_div * p_dtime - ), - axis=C2EDim, - ) + pm_0 = maximum(zero, p_mflx_tracer_h(C2E[0]) * geofac_div(C2CE[0]) * p_dtime) + pm_1 = maximum(zero, p_mflx_tracer_h(C2E[1]) * geofac_div(C2CE[1]) * p_dtime) + pm_2 = maximum(zero, p_mflx_tracer_h(C2E[2]) * geofac_div(C2CE[2]) * p_dtime) + p_m = pm_0 + pm_1 + pm_2 r_m = minimum( broadcast(1.0, (CellDim, KDim)), (p_cc * p_rhodz_now) / (p_m + dbl_eps) ) + return r_m @program def hflx_limiter_pd_stencil_01( - geofac_div: Field[[CellDim, C2EDim], float], + geofac_div: Field[[CEDim], float], p_cc: Field[[CellDim, KDim], float], p_rhodz_now: Field[[CellDim, KDim], float], p_mflx_tracer_h: Field[[EdgeDim, KDim], float], diff --git a/advection/tests/test_hflx_limiter_pd_stencil_01.py b/advection/tests/test_hflx_limiter_pd_stencil_01.py index 6cd13adde..6e7d85da9 100644 --- a/advection/tests/test_hflx_limiter_pd_stencil_01.py +++ b/advection/tests/test_hflx_limiter_pd_stencil_01.py @@ -12,13 +12,14 @@ # SPDX-License-Identifier: GPL-3.0-or-later import numpy as np +from functional.iterator.embedded import StridedNeighborOffsetProvider from icon4py.advection.hflx_limiter_pd_stencil_01 import ( hflx_limiter_pd_stencil_01, ) -from icon4py.common.dimension import C2EDim, CellDim, EdgeDim, KDim +from icon4py.common.dimension import C2EDim, CEDim, CellDim, EdgeDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh -from icon4py.testutils.utils import random_field, zero_field +from icon4py.testutils.utils import as_1D_sparse_field, random_field, zero_field def hflx_limiter_pd_stencil_01_numpy( @@ -30,20 +31,13 @@ def hflx_limiter_pd_stencil_01_numpy( p_dtime, dbl_eps, ): - # p_m: np.array - p_mflx_tracer_h_c2e = p_mflx_tracer_h[c2e] - # p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) + max(0.0, geofac_div[:, 1]*p_mflx_tracer_h_c2e[:, 1]*p_dtime) + max(0.0, geofac_div[:, 2]*p_mflx_tracer_h_c2e[:, 2]*p_dtime) geofac_div = np.expand_dims(geofac_div, axis=-1) - # p_m = max(0.0, geofac_div[:, 0]*p_mflx_tracer_h_c2e[:, 0]*p_dtime) - p_m = np.sum( - np.maximum( - np.zeros(p_mflx_tracer_h_c2e.shape, dtype=float), - p_mflx_tracer_h_c2e * geofac_div * p_dtime, - ) - ) - r_m = np.minimum( - np.ones(p_cc.shape, dtype=float), (p_cc * p_rhodz_now) / (p_m + dbl_eps) - ) + p_m_0 = np.maximum(0.0, p_mflx_tracer_h[c2e[:, 0]] * geofac_div[:, 0] * p_dtime) + p_m_1 = np.maximum(0.0, p_mflx_tracer_h[c2e[:, 1]] * geofac_div[:, 1] * p_dtime) + p_m_2 = np.maximum(0.0, p_mflx_tracer_h[c2e[:, 2]] * geofac_div[:, 2] * p_dtime) + + p_m = p_m_0 + p_m_1 + p_m_2 + r_m = np.minimum(1.0, p_cc * p_rhodz_now / (p_m + dbl_eps)) return r_m @@ -69,7 +63,7 @@ def test_hflx_limiter_pd_stencil_01(): ) hflx_limiter_pd_stencil_01( - geofac_div, + as_1D_sparse_field(geofac_div, CEDim), p_cc, p_rhodz_now, p_mflx_tracer_h, @@ -77,6 +71,7 @@ def test_hflx_limiter_pd_stencil_01(): p_dtime, dbl_eps, offset_provider={ + "C2CE": StridedNeighborOffsetProvider(CellDim, CEDim, mesh.n_c2e), "C2E": mesh.get_c2e_offset_provider(), }, ) From ffda8359b2f578ebf4035575dcad1b120014876f Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Tue, 20 Sep 2022 14:22:36 +0200 Subject: [PATCH 04/10] hflx_limiter_mo_stencil_01 part1 and vert_adv_stencil_01 --- .../advection/hflx_limiter_mo_stencil_01.py | 48 +++++++++++++ .../icon4py/advection/vert_adv_stencil_01.py | 55 ++++++++++++++ .../tests/test_hflx_limiter_mo_stencil_01.py | 66 +++++++++++++++++ advection/tests/test_vert_adv_stencil_01.py | 72 +++++++++++++++++++ 4 files changed, 241 insertions(+) create mode 100644 advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py create mode 100644 advection/src/icon4py/advection/vert_adv_stencil_01.py create mode 100644 advection/tests/test_hflx_limiter_mo_stencil_01.py create mode 100644 advection/tests/test_vert_adv_stencil_01.py diff --git a/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py b/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py new file mode 100644 index 000000000..1832bf6ed --- /dev/null +++ b/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py @@ -0,0 +1,48 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field, abs + +from icon4py.common.dimension import C2CE, E2C, CEDim, CellDim, EdgeDim, KDim + + +@field_operator +def _hflx_limiter_mo_stencil_01( + p_mflx_tracer_h: Field[[EdgeDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], + p_mass_flx_e: Field[[EdgeDim, KDim], float], +) -> tuple[Field[[EdgeDim, KDim], float],Field[[EdgeDim, KDim], float]]: + + z_mflx_low = 0.5 * ( p_mass_flx_e * (p_cc(E2C[0]) + p_cc(E2C[1])) + - abs(p_mass_flx_e) * (-p_cc(E2C[0]) + p_cc(E2C[1]))) + + z_anti = p_mflx_tracer_h - z_mflx_low + + return z_mflx_low, z_anti + + +@program +def hflx_limiter_mo_stencil_01( + p_mflx_tracer_h: Field[[EdgeDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], + p_mass_flx_e: Field[[EdgeDim, KDim], float], + z_mflx_low: Field[[EdgeDim, KDim], float], + z_anti: Field[[EdgeDim, KDim], float], +): + _hflx_limiter_mo_stencil_01( + p_mflx_tracer_h, + p_cc, + p_mass_flx_e, + out = (z_mflx_low, z_anti), + ) diff --git a/advection/src/icon4py/advection/vert_adv_stencil_01.py b/advection/src/icon4py/advection/vert_adv_stencil_01.py new file mode 100644 index 000000000..da9be0234 --- /dev/null +++ b/advection/src/icon4py/advection/vert_adv_stencil_01.py @@ -0,0 +1,55 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field + +from icon4py.common.dimension import CellDim, Koff, KDim + + +@field_operator +def _vert_adv_stencil_01( + tracer_now: Field[[CellDim, KDim], float], + rhodz_now: Field[[CellDim, KDim], float], + p_mflx_tracer_v: Field[[CellDim, KDim], float], + deepatmo_divzl: Field[[KDim], float], + deepatmo_divzu: Field[[KDim], float], + rhodz_new: Field[[CellDim, KDim], float], + p_dtime: float, +) -> Field[[CellDim, KDim], float]: + tracer_new = ( tracer_now * rhodz_now + p_dtime * ( p_mflx_tracer_v(Koff[1]) * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu ) ) / rhodz_new + + return tracer_new + + +@program +def vert_adv_stencil_01( + tracer_now: Field[[CellDim, KDim], float], + rhodz_now: Field[[CellDim, KDim], float], + p_mflx_tracer_v: Field[[CellDim, KDim], float], + deepatmo_divzl: Field[[KDim], float], + deepatmo_divzu: Field[[KDim], float], + rhodz_new: Field[[CellDim, KDim], float], + tracer_new: Field[[CellDim, KDim], float], + p_dtime: float, +): + _vert_adv_stencil_01( + tracer_now, + rhodz_now, + p_mflx_tracer_v, + deepatmo_divzl, + deepatmo_divzu, + rhodz_new, + p_dtime, + out=tracer_new, + ) diff --git a/advection/tests/test_hflx_limiter_mo_stencil_01.py b/advection/tests/test_hflx_limiter_mo_stencil_01.py new file mode 100644 index 000000000..c3745132c --- /dev/null +++ b/advection/tests/test_hflx_limiter_mo_stencil_01.py @@ -0,0 +1,66 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np +from functional.iterator.embedded import StridedNeighborOffsetProvider + +from icon4py.advection.hflx_limiter_mo_stencil_01 import ( + hflx_limiter_mo_stencil_01, +) +from icon4py.common.dimension import C2EDim, CEDim, CellDim, EdgeDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field + + +def hflx_limiter_mo_stencil_01_numpy( + e2c: np.array, + p_mflx_tracer_h: np.array, + p_cc: np.array, + p_mass_flx_e: np.array, +): + p_cc_e2c = p_cc[e2c] + z_mflx_low = 0.5*(p_mass_flx_e*(p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) + - abs(p_mass_flx_e) *(-p_cc_e2c[:, 0] + p_cc_e2c[:, 1])) + + z_anti = p_mflx_tracer_h - z_mflx_low + + return z_mflx_low, z_anti + + +def test_hflx_limiter_mo_stencil_01(): + mesh = SimpleMesh() + p_mflx_tracer_h = random_field(mesh, EdgeDim, KDim) + p_cc = random_field(mesh, CellDim, KDim) + p_mass_flx_e = random_field(mesh, EdgeDim, KDim) + z_mflx_low = random_field(mesh, EdgeDim, KDim) + z_anti = random_field(mesh, EdgeDim, KDim) + + z_mflx_low_ref, z_anti_ref = hflx_limiter_mo_stencil_01_numpy( + mesh.e2c, + np.asarray(p_mflx_tracer_h), + np.asarray(p_cc), + np.asarray(p_mass_flx_e), + ) + + hflx_limiter_mo_stencil_01( + p_mflx_tracer_h, + p_cc, + p_mass_flx_e, + z_mflx_low, + z_anti, + offset_provider={ + "E2C": mesh.get_e2c_offset_provider(), + }, + ) + assert np.allclose(z_mflx_low, z_mflx_low_ref) + assert np.allclose(z_anti, z_anti_ref) diff --git a/advection/tests/test_vert_adv_stencil_01.py b/advection/tests/test_vert_adv_stencil_01.py new file mode 100644 index 000000000..3761be12f --- /dev/null +++ b/advection/tests/test_vert_adv_stencil_01.py @@ -0,0 +1,72 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.advection.vert_adv_stencil_01 import vert_adv_stencil_01 +from icon4py.common.dimension import CellDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field, zero_field + + +def vert_adv_stencil_01_numpy( + tracer_now: np.array, + rhodz_now: np.array, + p_mflx_tracer_v: np.array, + deepatmo_divzl: np.array, + deepatmo_divzu: np.array, + rhodz_new: np.array, + p_dtime, +) -> np.array: + p_mflx_tracer_v_offset_1 = np.roll(p_mflx_tracer_v, shift=-1, axis=1) + + tracer_new = (tracer_now * rhodz_now + p_dtime * ( + p_mflx_tracer_v_offset_1 * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu)) / rhodz_new + + + return tracer_new + + +def test_vert_adv_stencil_01(): + mesh = SimpleMesh() + + tracer_now = random_field(mesh, CellDim, KDim) + rhodz_now = random_field(mesh, CellDim, KDim) + p_mflx_tracer_v = random_field(mesh, CellDim, KDim) + deepatmo_divzl = random_field(mesh, KDim) + deepatmo_divzu = random_field(mesh, KDim) + rhodz_new = random_field(mesh, CellDim, KDim) + tracer_new = random_field(mesh, CellDim, KDim) + p_dtime = np.float64(5.0) + + ref = vert_adv_stencil_01_numpy( + np.asarray(tracer_now), + np.asarray(rhodz_now), + np.asarray(p_mflx_tracer_v), + np.asarray(deepatmo_divzl), + np.asarray(deepatmo_divzu), + np.asarray(rhodz_new), + p_dtime, + ) + vert_adv_stencil_01( + tracer_now, + rhodz_now, + p_mflx_tracer_v, + deepatmo_divzl, + deepatmo_divzu, + rhodz_new, + tracer_new, + p_dtime, + offset_provider={"Koff": KDim}, + ) + assert np.allclose(tracer_new[:,:-1], ref[:,:-1]) From 4505c2683b6168e82d526e35bdf6c2fef28be313 Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Tue, 20 Sep 2022 16:12:23 +0200 Subject: [PATCH 05/10] adding stencils --- .../advection/hflx_limiter_mo_stencil_01.py | 12 ++- .../icon4py/advection/hor_adv_stencil_01.py | 68 ++++++++++++++ .../advection/step_advection_stencil_03.py | 39 ++++++++ .../upwind_hflux_miura_stencil_01.py | 70 ++++++++++++++ .../icon4py/advection/vert_adv_stencil_01.py | 8 +- .../tests/test_hflx_limiter_mo_stencil_01.py | 9 +- advection/tests/test_hor_adv_stencil_01.py | 80 ++++++++++++++++ .../tests/test_step_advection_stencil_03.py | 55 +++++++++++ .../test_upwind_hflux_miura_stencil_01.py | 94 +++++++++++++++++++ advection/tests/test_vert_adv_stencil_01.py | 12 ++- 10 files changed, 431 insertions(+), 16 deletions(-) create mode 100644 advection/src/icon4py/advection/hor_adv_stencil_01.py create mode 100644 advection/src/icon4py/advection/step_advection_stencil_03.py create mode 100644 advection/src/icon4py/advection/upwind_hflux_miura_stencil_01.py create mode 100644 advection/tests/test_hor_adv_stencil_01.py create mode 100644 advection/tests/test_step_advection_stencil_03.py create mode 100644 advection/tests/test_upwind_hflux_miura_stencil_01.py diff --git a/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py b/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py index 1832bf6ed..c3309e2ce 100644 --- a/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py +++ b/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py @@ -14,7 +14,7 @@ from functional.ffront.decorator import field_operator, program from functional.ffront.fbuiltins import Field, abs -from icon4py.common.dimension import C2CE, E2C, CEDim, CellDim, EdgeDim, KDim +from icon4py.common.dimension import E2C, CellDim, EdgeDim, KDim @field_operator @@ -22,10 +22,12 @@ def _hflx_limiter_mo_stencil_01( p_mflx_tracer_h: Field[[EdgeDim, KDim], float], p_cc: Field[[CellDim, KDim], float], p_mass_flx_e: Field[[EdgeDim, KDim], float], -) -> tuple[Field[[EdgeDim, KDim], float],Field[[EdgeDim, KDim], float]]: +) -> tuple[Field[[EdgeDim, KDim], float], Field[[EdgeDim, KDim], float]]: - z_mflx_low = 0.5 * ( p_mass_flx_e * (p_cc(E2C[0]) + p_cc(E2C[1])) - - abs(p_mass_flx_e) * (-p_cc(E2C[0]) + p_cc(E2C[1]))) + z_mflx_low = 0.5 * ( + p_mass_flx_e * (p_cc(E2C[0]) + p_cc(E2C[1])) + - abs(p_mass_flx_e) * (-p_cc(E2C[0]) + p_cc(E2C[1])) + ) z_anti = p_mflx_tracer_h - z_mflx_low @@ -44,5 +46,5 @@ def hflx_limiter_mo_stencil_01( p_mflx_tracer_h, p_cc, p_mass_flx_e, - out = (z_mflx_low, z_anti), + out=(z_mflx_low, z_anti), ) diff --git a/advection/src/icon4py/advection/hor_adv_stencil_01.py b/advection/src/icon4py/advection/hor_adv_stencil_01.py new file mode 100644 index 000000000..2cb692a10 --- /dev/null +++ b/advection/src/icon4py/advection/hor_adv_stencil_01.py @@ -0,0 +1,68 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field, neighbor_sum + +from icon4py.common.dimension import ( + C2CE, + C2E, + C2EDim, + CEDim, + CellDim, + EdgeDim, + KDim, +) + + +@field_operator +def _hor_adv_stencil_01( + p_mflx_tracer_h: Field[[EdgeDim, KDim], float], + deepatmo_divh: Field[[KDim], float], + tracer_now: Field[[CellDim, KDim], float], + rhodz_now: Field[[CellDim, KDim], float], + rhodz_new: Field[[CellDim, KDim], float], + geofac_div: Field[[CEDim], float], + p_dtime: float, +) -> Field[[CellDim, KDim], float]: + tracer_new = ( + tracer_now * rhodz_now + - p_dtime + * deepatmo_divh + * neighbor_sum(p_mflx_tracer_h(C2E) * geofac_div(C2CE), axis=C2EDim) + ) / rhodz_new + + return tracer_new + + +@program +def hor_adv_stencil_01( + p_mflx_tracer_h: Field[[EdgeDim, KDim], float], + deepatmo_divh: Field[[KDim], float], + tracer_now: Field[[CellDim, KDim], float], + rhodz_now: Field[[CellDim, KDim], float], + rhodz_new: Field[[CellDim, KDim], float], + geofac_div: Field[[CEDim], float], + tracer_new: Field[[CellDim, KDim], float], + p_dtime: float, +): + _hor_adv_stencil_01( + p_mflx_tracer_h, + deepatmo_divh, + tracer_now, + rhodz_now, + rhodz_new, + geofac_div, + p_dtime, + out=tracer_new, + ) diff --git a/advection/src/icon4py/advection/step_advection_stencil_03.py b/advection/src/icon4py/advection/step_advection_stencil_03.py new file mode 100644 index 000000000..2c3728ca0 --- /dev/null +++ b/advection/src/icon4py/advection/step_advection_stencil_03.py @@ -0,0 +1,39 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field, maximum + +from icon4py.common.dimension import CellDim, KDim + + +@field_operator +def _step_advection_stencil_03( + p_tracer_now: Field[[CellDim, KDim], float], + p_grf_tend_tracer: Field[[CellDim, KDim], float], + p_dtime: float, +) -> Field[[CellDim, KDim], float]: + p_tracer_new = maximum(0.0, p_tracer_now + p_dtime * p_grf_tend_tracer) + return p_tracer_new + + +@program +def step_advection_stencil_03( + p_tracer_now: Field[[CellDim, KDim], float], + p_grf_tend_tracer: Field[[CellDim, KDim], float], + p_tracer_new: Field[[CellDim, KDim], float], + p_dtime: float, +): + _step_advection_stencil_03( + p_tracer_now, p_grf_tend_tracer, p_dtime, out=p_tracer_new + ) diff --git a/advection/src/icon4py/advection/upwind_hflux_miura_stencil_01.py b/advection/src/icon4py/advection/upwind_hflux_miura_stencil_01.py new file mode 100644 index 000000000..31c44a66c --- /dev/null +++ b/advection/src/icon4py/advection/upwind_hflux_miura_stencil_01.py @@ -0,0 +1,70 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field, where + +from icon4py.common.dimension import E2C, CellDim, EdgeDim, KDim + + +@field_operator +def _upwind_hflux_miura_stencil_01( + p_vn: Field[[EdgeDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], + distv_bary_1: Field[[EdgeDim, KDim], float], + distv_bary_2: Field[[EdgeDim, KDim], float], + z_grad_1: Field[[CellDim, KDim], float], + z_grad_2: Field[[CellDim, KDim], float], + p_mass_flx_e: Field[[EdgeDim, KDim], float], +) -> Field[[EdgeDim, KDim], float]: + + p_out_e = where( + p_vn > 0.0, + ( + p_cc(E2C[0]) + + distv_bary_1 * z_grad_1(E2C[0]) + + distv_bary_2 * z_grad_2(E2C[0]) + ) + * p_mass_flx_e, + ( + p_cc(E2C[1]) + + distv_bary_1 * z_grad_1(E2C[1]) + + distv_bary_2 * z_grad_2(E2C[1]) + ) + * p_mass_flx_e, + ) + + return p_out_e + + +@program +def upwind_hflux_miura_stencil_01( + p_vn: Field[[EdgeDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], + distv_bary_1: Field[[EdgeDim, KDim], float], + distv_bary_2: Field[[EdgeDim, KDim], float], + z_grad_1: Field[[CellDim, KDim], float], + z_grad_2: Field[[CellDim, KDim], float], + p_mass_flx_e: Field[[EdgeDim, KDim], float], + p_out_e: Field[[EdgeDim, KDim], float], +): + _upwind_hflux_miura_stencil_01( + p_vn, + p_cc, + distv_bary_1, + distv_bary_2, + z_grad_1, + z_grad_2, + p_mass_flx_e, + out=p_out_e, + ) diff --git a/advection/src/icon4py/advection/vert_adv_stencil_01.py b/advection/src/icon4py/advection/vert_adv_stencil_01.py index da9be0234..b4a87dd67 100644 --- a/advection/src/icon4py/advection/vert_adv_stencil_01.py +++ b/advection/src/icon4py/advection/vert_adv_stencil_01.py @@ -14,7 +14,7 @@ from functional.ffront.decorator import field_operator, program from functional.ffront.fbuiltins import Field -from icon4py.common.dimension import CellDim, Koff, KDim +from icon4py.common.dimension import CellDim, KDim, Koff @field_operator @@ -27,7 +27,11 @@ def _vert_adv_stencil_01( rhodz_new: Field[[CellDim, KDim], float], p_dtime: float, ) -> Field[[CellDim, KDim], float]: - tracer_new = ( tracer_now * rhodz_now + p_dtime * ( p_mflx_tracer_v(Koff[1]) * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu ) ) / rhodz_new + tracer_new = ( + tracer_now * rhodz_now + + p_dtime + * (p_mflx_tracer_v(Koff[1]) * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu) + ) / rhodz_new return tracer_new diff --git a/advection/tests/test_hflx_limiter_mo_stencil_01.py b/advection/tests/test_hflx_limiter_mo_stencil_01.py index c3745132c..5423f15e0 100644 --- a/advection/tests/test_hflx_limiter_mo_stencil_01.py +++ b/advection/tests/test_hflx_limiter_mo_stencil_01.py @@ -12,12 +12,11 @@ # SPDX-License-Identifier: GPL-3.0-or-later import numpy as np -from functional.iterator.embedded import StridedNeighborOffsetProvider from icon4py.advection.hflx_limiter_mo_stencil_01 import ( hflx_limiter_mo_stencil_01, ) -from icon4py.common.dimension import C2EDim, CEDim, CellDim, EdgeDim, KDim +from icon4py.common.dimension import CellDim, EdgeDim, KDim from icon4py.testutils.simple_mesh import SimpleMesh from icon4py.testutils.utils import random_field @@ -29,8 +28,10 @@ def hflx_limiter_mo_stencil_01_numpy( p_mass_flx_e: np.array, ): p_cc_e2c = p_cc[e2c] - z_mflx_low = 0.5*(p_mass_flx_e*(p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) - - abs(p_mass_flx_e) *(-p_cc_e2c[:, 0] + p_cc_e2c[:, 1])) + z_mflx_low = 0.5 * ( + p_mass_flx_e * (p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) + - abs(p_mass_flx_e) * (-p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) + ) z_anti = p_mflx_tracer_h - z_mflx_low diff --git a/advection/tests/test_hor_adv_stencil_01.py b/advection/tests/test_hor_adv_stencil_01.py new file mode 100644 index 000000000..359840230 --- /dev/null +++ b/advection/tests/test_hor_adv_stencil_01.py @@ -0,0 +1,80 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np +from functional.iterator.embedded import StridedNeighborOffsetProvider + +from icon4py.advection.hor_adv_stencil_01 import hor_adv_stencil_01 +from icon4py.common.dimension import C2EDim, CEDim, CellDim, EdgeDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import as_1D_sparse_field, random_field + + +def hor_adv_stencil_01_numpy( + c2e: np.array, + p_mflx_tracer_h: np.array, + deepatmo_divh: np.array, + tracer_now: np.array, + rhodz_now: np.array, + rhodz_new: np.array, + geofac_div: np.array, + p_dtime, +) -> np.array: + geofac_div = np.expand_dims(geofac_div, axis=-1) + + tracer_new = ( + tracer_now * rhodz_now + - p_dtime * deepatmo_divh * np.sum(p_mflx_tracer_h[c2e] * geofac_div, axis=1) + ) / rhodz_new + + return tracer_new + + +def test_hor_adv_stencil_01(): + mesh = SimpleMesh() + + p_mflx_tracer_h = random_field(mesh, EdgeDim, KDim) + deepatmo_divh = random_field(mesh, KDim) + tracer_now = random_field(mesh, CellDim, KDim) + rhodz_now = random_field(mesh, CellDim, KDim) + rhodz_new = random_field(mesh, CellDim, KDim) + geofac_div = random_field(mesh, CellDim, C2EDim) + geofac_div_new = as_1D_sparse_field(geofac_div, CEDim) + tracer_new = random_field(mesh, CellDim, KDim) + p_dtime = np.float64(5.0) + + ref = hor_adv_stencil_01_numpy( + mesh.c2e, + np.asarray(p_mflx_tracer_h), + np.asarray(deepatmo_divh), + np.asarray(tracer_now), + np.asarray(rhodz_now), + np.asarray(rhodz_new), + np.asarray(geofac_div), + p_dtime, + ) + hor_adv_stencil_01( + p_mflx_tracer_h, + deepatmo_divh, + tracer_now, + rhodz_now, + rhodz_new, + geofac_div_new, + tracer_new, + p_dtime, + offset_provider={ + "C2E": mesh.get_c2e_offset_provider(), + "C2CE": StridedNeighborOffsetProvider(CellDim, CEDim, mesh.n_c2e), + }, + ) + assert np.allclose(tracer_new, ref) diff --git a/advection/tests/test_step_advection_stencil_03.py b/advection/tests/test_step_advection_stencil_03.py new file mode 100644 index 000000000..ade783cb0 --- /dev/null +++ b/advection/tests/test_step_advection_stencil_03.py @@ -0,0 +1,55 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.advection.step_advection_stencil_03 import step_advection_stencil_03 +from icon4py.common.dimension import CellDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field + + +def step_advection_stencil_03_numpy( + p_tracer_now: np.array, + p_grf_tend_tracer: np.array, + p_dtime, +) -> np.array: + + p_tracer_new = p_tracer_now + p_dtime * p_grf_tend_tracer + p_tracer_new = np.where(p_tracer_new < 0.0, 0.0, p_tracer_new) + + return p_tracer_new + + +def test_step_advection_stencil_03(): + mesh = SimpleMesh() + + p_tracer_now = random_field(mesh, CellDim, KDim) + p_grf_tend_tracer = random_field(mesh, CellDim, KDim) + p_tracer_new = random_field(mesh, CellDim, KDim) + + p_dtime = np.float64(5.0) + + ref = step_advection_stencil_03_numpy( + np.asarray(p_tracer_now), + np.asarray(p_grf_tend_tracer), + p_dtime, + ) + step_advection_stencil_03( + p_tracer_now, + p_grf_tend_tracer, + p_tracer_new, + p_dtime, + offset_provider={}, + ) + assert np.allclose(p_tracer_new, ref) diff --git a/advection/tests/test_upwind_hflux_miura_stencil_01.py b/advection/tests/test_upwind_hflux_miura_stencil_01.py new file mode 100644 index 000000000..a4fa164a7 --- /dev/null +++ b/advection/tests/test_upwind_hflux_miura_stencil_01.py @@ -0,0 +1,94 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.advection.upwind_hflux_miura_stencil_01 import ( + upwind_hflux_miura_stencil_01, +) +from icon4py.common.dimension import CellDim, EdgeDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field, zero_field + + +def upwind_hflux_miura_stencil_01_numpy( + e2c: np.array, + p_vn: np.array, + p_cc: np.array, + distv_bary_1: np.array, + distv_bary_2: np.array, + z_grad_1: np.array, + z_grad_2: np.array, + p_mass_flx_e: np.array, +) -> np.array: + + p_cc_e2c = p_cc[e2c] + z_grad_1_e2c = z_grad_1[e2c] + z_grad_2_e2c = z_grad_2[e2c] + + p_out_e = np.where( + p_vn > 0.0, + ( + p_cc_e2c[:, 0] + + distv_bary_1 * z_grad_1_e2c[:, 0] + + distv_bary_2 * z_grad_2_e2c[:, 0] + ) + * p_mass_flx_e, + ( + p_cc_e2c[:, 1] + + distv_bary_1 * z_grad_1_e2c[:, 1] + + distv_bary_2 * z_grad_2_e2c[:, 1] + ) + * p_mass_flx_e, + ) + + return p_out_e + + +def test_upwind_hflux_miura_stencil_01(): + mesh = SimpleMesh() + + p_vn = random_field(mesh, EdgeDim, KDim) + p_cc = random_field(mesh, CellDim, KDim) + distv_bary_1 = random_field(mesh, EdgeDim, KDim) + distv_bary_2 = random_field(mesh, EdgeDim, KDim) + z_grad_1 = random_field(mesh, CellDim, KDim) + z_grad_2 = random_field(mesh, CellDim, KDim) + p_mass_flx_e = random_field(mesh, EdgeDim, KDim) + p_out_e = zero_field(mesh, EdgeDim, KDim) + + ref = upwind_hflux_miura_stencil_01_numpy( + mesh.e2c, + np.asarray(p_vn), + np.asarray(p_cc), + np.asarray(distv_bary_1), + np.asarray(distv_bary_2), + np.asarray(z_grad_1), + np.asarray(z_grad_2), + np.asarray(p_mass_flx_e), + ) + + upwind_hflux_miura_stencil_01( + p_vn, + p_cc, + distv_bary_1, + distv_bary_2, + z_grad_1, + z_grad_2, + p_mass_flx_e, + p_out_e, + offset_provider={ + "E2C": mesh.get_e2c_offset_provider(), + }, + ) + assert np.allclose(p_out_e, ref) diff --git a/advection/tests/test_vert_adv_stencil_01.py b/advection/tests/test_vert_adv_stencil_01.py index 3761be12f..c8de7c3dc 100644 --- a/advection/tests/test_vert_adv_stencil_01.py +++ b/advection/tests/test_vert_adv_stencil_01.py @@ -30,9 +30,11 @@ def vert_adv_stencil_01_numpy( ) -> np.array: p_mflx_tracer_v_offset_1 = np.roll(p_mflx_tracer_v, shift=-1, axis=1) - tracer_new = (tracer_now * rhodz_now + p_dtime * ( - p_mflx_tracer_v_offset_1 * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu)) / rhodz_new - + tracer_new = ( + tracer_now * rhodz_now + + p_dtime + * (p_mflx_tracer_v_offset_1 * deepatmo_divzl - p_mflx_tracer_v * deepatmo_divzu) + ) / rhodz_new return tracer_new @@ -46,7 +48,7 @@ def test_vert_adv_stencil_01(): deepatmo_divzl = random_field(mesh, KDim) deepatmo_divzu = random_field(mesh, KDim) rhodz_new = random_field(mesh, CellDim, KDim) - tracer_new = random_field(mesh, CellDim, KDim) + tracer_new = zero_field(mesh, CellDim, KDim) p_dtime = np.float64(5.0) ref = vert_adv_stencil_01_numpy( @@ -69,4 +71,4 @@ def test_vert_adv_stencil_01(): p_dtime, offset_provider={"Koff": KDim}, ) - assert np.allclose(tracer_new[:,:-1], ref[:,:-1]) + assert np.allclose(tracer_new[:, :-1], ref[:, :-1]) From a1dc59c72dd73214eef636c3f2cb51fd2df8a74c Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Mon, 3 Oct 2022 14:37:19 +0200 Subject: [PATCH 06/10] Adding stencils face_val_ppm_stencil_05 and v_limit_prbl_sm_stencil_01 --- .../advection/face_val_ppm_stencil_05.py | 67 +++++++++++++++ .../advection/v_limit_prbl_sm_stencil_01.py | 70 ++++++++++++++++ .../tests/test_face_val_ppm_stencil_05.py | 83 +++++++++++++++++++ .../tests/test_vlimit_prbl_sm_stencil_01.py | 75 +++++++++++++++++ 4 files changed, 295 insertions(+) create mode 100644 advection/src/icon4py/advection/face_val_ppm_stencil_05.py create mode 100644 advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py create mode 100644 advection/tests/test_face_val_ppm_stencil_05.py create mode 100644 advection/tests/test_vlimit_prbl_sm_stencil_01.py diff --git a/advection/src/icon4py/advection/face_val_ppm_stencil_05.py b/advection/src/icon4py/advection/face_val_ppm_stencil_05.py new file mode 100644 index 000000000..7f2f91cd7 --- /dev/null +++ b/advection/src/icon4py/advection/face_val_ppm_stencil_05.py @@ -0,0 +1,67 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field + +from icon4py.common.dimension import CellDim, KDim, Koff + + +@field_operator +def _face_val_ppm_stencil_05( + p_cc: Field[[CellDim, KDim], float], + p_cellhgt_mc_now: Field[[CellDim, KDim], float], + z_slope: Field[[CellDim, KDim], float], +) -> Field[[CellDim, KDim], float]: + + zgeo1 = p_cellhgt_mc_now(Koff[-1]) / (p_cellhgt_mc_now(Koff[-1]) + p_cellhgt_mc_now) + zgeo2 = 1.0 / ( + p_cellhgt_mc_now(Koff[-2]) + + p_cellhgt_mc_now(Koff[-1]) + + p_cellhgt_mc_now + + p_cellhgt_mc_now(Koff[1]) + ) + zgeo3 = (p_cellhgt_mc_now(Koff[-2]) + p_cellhgt_mc_now(Koff[-1])) / ( + 2.0 * p_cellhgt_mc_now(Koff[-1]) + p_cellhgt_mc_now + ) + zgeo4 = (p_cellhgt_mc_now(Koff[1]) + p_cellhgt_mc_now) / ( + 2.0 * p_cellhgt_mc_now + p_cellhgt_mc_now(Koff[-1]) + ) + + p_face = ( + p_cc(Koff[-1]) + + zgeo1 * (p_cc - p_cc(Koff[-1])) + + zgeo2 + * ( + (2.0 * p_cellhgt_mc_now * zgeo1) * (zgeo3 - zgeo4) * (p_cc - p_cc(Koff[-1])) + - zgeo3 * p_cellhgt_mc_now(Koff[-1]) * z_slope + + zgeo4 * p_cellhgt_mc_now * z_slope(Koff[-1]) + ) + ) + + return p_face + + +@program +def face_val_ppm_stencil_05( + p_cc: Field[[CellDim, KDim], float], + p_cellhgt_mc_now: Field[[CellDim, KDim], float], + z_slope: Field[[CellDim, KDim], float], + p_face: Field[[CellDim, KDim], float], +): + _face_val_ppm_stencil_05( + p_cc, + p_cellhgt_mc_now, + z_slope, + out=p_face, + ) diff --git a/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py b/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py new file mode 100644 index 000000000..ce21024a1 --- /dev/null +++ b/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py @@ -0,0 +1,70 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import ( + Field, + FieldOffset, + abs, + broadcast, + minimum, + where, +) + +from icon4py.common.dimension import CellDim, KDim + + +Koff = FieldOffset("Koff", source=KDim, target=(KDim,)) + + +@field_operator +def _v_limit_prbl_sm_stencil_01( + p_face: Field[[CellDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], +) -> Field[[CellDim, KDim], float]: + + z_delta = p_face - p_face(Koff[1]) + z_a6i = -6.0 * (p_cc - 0.5 * (p_face + p_face(Koff[1]))) + + q_face_up = broadcast(1.0, (CellDim, KDim)) + q_face_low = broadcast(1.0, (CellDim, KDim)) + + q_face_up, q_face_low = where( + abs(z_delta) < z_a6i, + where( + (p_cc < minimum(p_face, p_face(Koff[1]))), + (p_cc, p_cc), + where( + p_face > p_face(Koff[1]), + (3.0 * p_cc - 2.0 * p_face(Koff[1]), p_face(Koff[1])), + (p_face, 3.0 * p_cc - 2.0 * p_face), + ), + ), + (p_face, p_face(Koff[-1])), + ) + + return q_face_up, q_face_low + + +@program +def v_limit_prbl_sm_stencil_01( + p_face: Field[[CellDim, KDim], float], + p_cc: Field[[CellDim, KDim], float], + p_face_up: Field[[CellDim, KDim], float], + p_face_low: Field[[CellDim, KDim], float], +): + _v_limit_prbl_sm_stencil_01( + p_face, + p_cc, + out=(p_face_up, p_face_low), + ) diff --git a/advection/tests/test_face_val_ppm_stencil_05.py b/advection/tests/test_face_val_ppm_stencil_05.py new file mode 100644 index 000000000..c36690a5c --- /dev/null +++ b/advection/tests/test_face_val_ppm_stencil_05.py @@ -0,0 +1,83 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.advection.face_val_ppm_stencil_05 import face_val_ppm_stencil_05 +from icon4py.common.dimension import CellDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field, zero_field + + +def face_val_ppm_stencil_05_numpy( + p_cc: np.array, + p_cellhgt_mc_now: np.array, + z_slope: np.array, +): + p_cellhgt_mc_now_k_minus_1 = np.roll(p_cellhgt_mc_now, shift=1, axis=1) + p_cellhgt_mc_now_k_minus_2 = np.roll(p_cellhgt_mc_now, shift=2, axis=1) + p_cellhgt_mc_now_k_plus_1 = np.roll(p_cellhgt_mc_now, shift=-1, axis=1) + + p_cc_k_minus_1 = np.roll(p_cc, shift=1, axis=1) + z_slope_k_minus_1 = np.roll(z_slope, shift=1, axis=1) + + zgeo1 = p_cellhgt_mc_now_k_minus_1 / (p_cellhgt_mc_now_k_minus_1 + p_cellhgt_mc_now) + zgeo2 = 1.0 / ( + p_cellhgt_mc_now_k_minus_2 + + p_cellhgt_mc_now_k_minus_1 + + p_cellhgt_mc_now + + p_cellhgt_mc_now_k_plus_1 + ) + zgeo3 = (p_cellhgt_mc_now_k_minus_2 + p_cellhgt_mc_now_k_minus_1) / ( + 2.0 * p_cellhgt_mc_now_k_minus_1 + p_cellhgt_mc_now + ) + zgeo4 = (p_cellhgt_mc_now_k_plus_1 + p_cellhgt_mc_now) / ( + 2 * p_cellhgt_mc_now + p_cellhgt_mc_now_k_minus_1 + ) + + p_face = ( + p_cc_k_minus_1 + + zgeo1 * (p_cc - p_cc_k_minus_1) + + zgeo2 + * ( + (2 * p_cellhgt_mc_now * zgeo1) * (zgeo3 - zgeo4) * (p_cc - p_cc_k_minus_1) + - zgeo3 * p_cellhgt_mc_now_k_minus_1 * z_slope + + zgeo4 * p_cellhgt_mc_now * z_slope_k_minus_1 + ) + ) + + return p_face + + +def test_face_val_ppm_stencil_05(): + mesh = SimpleMesh() + p_cc = random_field(mesh, CellDim, KDim) + p_cellhgt_mc_now = random_field(mesh, CellDim, KDim) + z_slope = random_field(mesh, CellDim, KDim) + p_face = zero_field(mesh, CellDim, KDim) + + ref = face_val_ppm_stencil_05_numpy( + np.asarray(p_cc), + np.asarray(p_cellhgt_mc_now), + np.asarray(z_slope), + ) + + face_val_ppm_stencil_05( + p_cc, + p_cellhgt_mc_now, + z_slope, + p_face, + offset_provider={"Koff": KDim}, + ) + + assert np.allclose(ref[:, 1:-1], p_face.__array__()[:, 1:-1]) diff --git a/advection/tests/test_vlimit_prbl_sm_stencil_01.py b/advection/tests/test_vlimit_prbl_sm_stencil_01.py new file mode 100644 index 000000000..e9313f79d --- /dev/null +++ b/advection/tests/test_vlimit_prbl_sm_stencil_01.py @@ -0,0 +1,75 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.advection.v_limit_prbl_sm_stencil_01 import ( + v_limit_prbl_sm_stencil_01, +) +from icon4py.common.dimension import CellDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field, zero_field + + +def v_limit_prbl_sm_stencil_01_numpy( + p_face: np.array, + p_cc: np.array, +): + p_face_k_minus_1 = np.roll(p_face, shift=1, axis=1) + p_face_k_plus_1 = np.roll(p_face, shift=-1, axis=1) + + z_delta = p_face - p_face_k_plus_1 + z_a6i = 6.0 * (p_cc - 0.5 * (p_face + p_face_k_plus_1)) + + q_face_up = 1.0 + q_face_low = 1.0 + + q_face_up, q_face_low = np.where( + abs(z_delta) < -1 * z_a6i, + np.where( + (p_cc < np.minimum(p_face, p_face_k_plus_1)), + (p_cc, p_cc), + np.where( + p_face > p_face_k_plus_1, + (3.0 * p_cc - 2.0 * p_face_k_plus_1, p_face_k_plus_1), + (p_face, 3.0 * p_cc - 2.0 * p_face), + ), + ), + (p_face, p_face_k_minus_1), + ) + + return q_face_up, q_face_low + + +def test_v_limit_prbl_sm_stencil_01(): + mesh = SimpleMesh() + p_cc = random_field(mesh, CellDim, KDim) + p_face = random_field(mesh, CellDim, KDim) + p_face_up = zero_field(mesh, CellDim, KDim) + p_face_low = zero_field(mesh, CellDim, KDim) + + p_face_up_ref, p_face_low_ref = v_limit_prbl_sm_stencil_01_numpy( + np.asarray(p_face), + np.asarray(p_cc), + ) + + v_limit_prbl_sm_stencil_01( + p_face, + p_cc, + p_face_up, + p_face_low, + offset_provider={"Koff": KDim}, + ) + + assert np.allclose(p_face_up_ref[:, :-1], p_face_up[:, :-1]) + assert np.allclose(p_face_low_ref[:, :-1], p_face_low[:, :-1]) From a2263b7fb873f17081404522ff9a1f71e57b8966 Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Tue, 4 Oct 2022 14:14:32 +0200 Subject: [PATCH 07/10] fix to k-offset --- .../src/icon4py/advection/v_limit_prbl_sm_stencil_01.py | 2 +- advection/tests/test_vlimit_prbl_sm_stencil_01.py | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py b/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py index ce21024a1..5d06900ec 100644 --- a/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py +++ b/advection/src/icon4py/advection/v_limit_prbl_sm_stencil_01.py @@ -50,7 +50,7 @@ def _v_limit_prbl_sm_stencil_01( (p_face, 3.0 * p_cc - 2.0 * p_face), ), ), - (p_face, p_face(Koff[-1])), + (p_face, p_face(Koff[1])), ) return q_face_up, q_face_low diff --git a/advection/tests/test_vlimit_prbl_sm_stencil_01.py b/advection/tests/test_vlimit_prbl_sm_stencil_01.py index e9313f79d..7367d9621 100644 --- a/advection/tests/test_vlimit_prbl_sm_stencil_01.py +++ b/advection/tests/test_vlimit_prbl_sm_stencil_01.py @@ -25,15 +25,11 @@ def v_limit_prbl_sm_stencil_01_numpy( p_face: np.array, p_cc: np.array, ): - p_face_k_minus_1 = np.roll(p_face, shift=1, axis=1) p_face_k_plus_1 = np.roll(p_face, shift=-1, axis=1) z_delta = p_face - p_face_k_plus_1 z_a6i = 6.0 * (p_cc - 0.5 * (p_face + p_face_k_plus_1)) - q_face_up = 1.0 - q_face_low = 1.0 - q_face_up, q_face_low = np.where( abs(z_delta) < -1 * z_a6i, np.where( @@ -45,7 +41,7 @@ def v_limit_prbl_sm_stencil_01_numpy( (p_face, 3.0 * p_cc - 2.0 * p_face), ), ), - (p_face, p_face_k_minus_1), + (p_face, p_face_k_plus_1), ) return q_face_up, q_face_low From b4bc1c9cdad13c6aad93fb14b14872d071e76b19 Mon Sep 17 00:00:00 2001 From: Nicoletta Farabullini <41536517+nfarabullini@users.noreply.github.com> Date: Tue, 4 Oct 2022 17:27:08 +0200 Subject: [PATCH 08/10] Tracer adv stencils with differing domain extents (#99) * face_val_ppm_stencil_04 * Redirected repo to functional following gt4py domain PR merge --- .../advection/face_val_ppm_stencil_04.py | 70 ++++++++++++++ .../tests/test_face_val_ppm_stencil_04.py | 91 +++++++++++++++++++ 2 files changed, 161 insertions(+) create mode 100644 advection/src/icon4py/advection/face_val_ppm_stencil_04.py create mode 100644 advection/tests/test_face_val_ppm_stencil_04.py diff --git a/advection/src/icon4py/advection/face_val_ppm_stencil_04.py b/advection/src/icon4py/advection/face_val_ppm_stencil_04.py new file mode 100644 index 000000000..0ddf3d3fb --- /dev/null +++ b/advection/src/icon4py/advection/face_val_ppm_stencil_04.py @@ -0,0 +1,70 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +from functional.ffront.decorator import field_operator, program +from functional.ffront.fbuiltins import Field, FieldOffset, broadcast + +from icon4py.common.dimension import CellDim, KDim + + +Koff = FieldOffset("Koff", source=KDim, target=(KDim,)) + + +@field_operator +def _face_val_ppm_stencil_04_p_face_1( + p_cc: Field[[CellDim, KDim], float], + p_cellhgt_mc_now: Field[[CellDim, KDim], float], +) -> Field[[CellDim, KDim], float]: + + p_face = p_cc * ( + broadcast(1.0, (CellDim, KDim)) + - (p_cellhgt_mc_now / p_cellhgt_mc_now(Koff[-1])) + ) + (p_cellhgt_mc_now / (p_cellhgt_mc_now(Koff[-1]) + p_cellhgt_mc_now)) * ( + (p_cellhgt_mc_now / p_cellhgt_mc_now(Koff[-1])) * p_cc + p_cc(Koff[-1]) + ) + return p_face + + +@field_operator +def _face_val_ppm_stencil_04_p_face_2( + p_cc: Field[[CellDim, KDim], float], +) -> Field[[CellDim, KDim], float]: + p_face = p_cc + return p_face + + +@program +def face_val_ppm_stencil_04( + p_cc: Field[[CellDim, KDim], float], + p_cellhgt_mc_now: Field[[CellDim, KDim], float], + p_face: Field[[CellDim, KDim], float], + nudging: int, + halo: int, + vertical_lower: int, + vertical_upper: int, +): + _face_val_ppm_stencil_04_p_face_1( + p_cc, + p_cellhgt_mc_now, + out=p_face, + domain={CellDim: (nudging, halo), KDim: (1, 2)}, + ) + _face_val_ppm_stencil_04_p_face_1( + p_cc, + p_cellhgt_mc_now, + out=p_face, + domain={CellDim: (nudging, halo), KDim: (2, vertical_upper)}, + ) + _face_val_ppm_stencil_04_p_face_2( + p_cc, out=p_face, domain={CellDim: (nudging, halo), KDim: (vertical_lower, 1)} + ) diff --git a/advection/tests/test_face_val_ppm_stencil_04.py b/advection/tests/test_face_val_ppm_stencil_04.py new file mode 100644 index 000000000..ae02704b3 --- /dev/null +++ b/advection/tests/test_face_val_ppm_stencil_04.py @@ -0,0 +1,91 @@ +# ICON4Py - ICON inspired code in Python and GT4Py +# +# Copyright (c) 2022, ETH Zurich and MeteoSwiss +# All rights reserved. +# +# This file is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or any later +# version. See the LICENSE.txt file at the top-level directory of this +# distribution for a copy of the license or check . +# +# SPDX-License-Identifier: GPL-3.0-or-later + +import numpy as np + +from icon4py.atm_dyn_iconam.face_val_ppm_stencil_04 import face_val_ppm_stencil_04 +from icon4py.common.dimension import CellDim, KDim +from icon4py.testutils.simple_mesh import SimpleMesh +from icon4py.testutils.utils import random_field, zero_field + + +def face_val_ppm_stencil_04_numpy( + p_cc: np.array, + p_cellhgt_mc_now: np.array, + nudging, + halo, + vertical_lower, + vertical_upper, +) -> np.array: + p_cellhgt_mc_now_minus_1 = np.roll(p_cellhgt_mc_now, shift=1, axis=1) + p_cc_minus_1 = np.roll(p_cc, shift=1, axis=1) + p_face_ref = np.zeros_like(p_cc) + + p_face_ref[nudging:halo, 1:vertical_upper] = p_cc[ + nudging:halo, 1:vertical_upper + ] * ( + 1 + - ( + p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] + / p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] + ) + ) + ( + p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] + / ( + p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] + + p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] + ) + ) * ( + ( + p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] + / p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] + ) + * p_cc[nudging:halo, 1:vertical_upper] + + p_cc_minus_1[nudging:halo, 1:vertical_upper] + ) + + p_face_ref[nudging:halo, vertical_lower] = p_cc[nudging:halo, vertical_lower] + return p_face_ref + + +def test_face_val_ppm_stencil_04(): + mesh = SimpleMesh() + + p_cc = random_field(mesh, CellDim, KDim) + p_cellhgt_mc_now = random_field(mesh, CellDim, KDim) + p_face = zero_field(mesh, CellDim, KDim) + nudging = 4 + halo = 16 + vertical_lower = 0 + vertical_upper = np.asarray(p_face).shape[1] + + p_face_ref = face_val_ppm_stencil_04_numpy( + np.asarray(p_cc), + np.asarray(p_cellhgt_mc_now), + nudging, + halo, + vertical_lower, + vertical_upper, + ) + face_val_ppm_stencil_04( + p_cc, + p_cellhgt_mc_now, + p_face, + nudging, + halo, + vertical_lower, + vertical_upper, + offset_provider={"Koff": KDim}, + ) + + assert np.allclose(p_face_ref, p_face) From b9106064c25cffebeed351438c99cbe06dfd5816 Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Mon, 10 Oct 2022 15:00:49 +0200 Subject: [PATCH 09/10] Revert "Tracer adv stencils with differing domain extents (#99)" This reverts commit b4bc1c9cdad13c6aad93fb14b14872d071e76b19. --- .../advection/face_val_ppm_stencil_04.py | 70 -------------- .../tests/test_face_val_ppm_stencil_04.py | 91 ------------------- 2 files changed, 161 deletions(-) delete mode 100644 advection/src/icon4py/advection/face_val_ppm_stencil_04.py delete mode 100644 advection/tests/test_face_val_ppm_stencil_04.py diff --git a/advection/src/icon4py/advection/face_val_ppm_stencil_04.py b/advection/src/icon4py/advection/face_val_ppm_stencil_04.py deleted file mode 100644 index 0ddf3d3fb..000000000 --- a/advection/src/icon4py/advection/face_val_ppm_stencil_04.py +++ /dev/null @@ -1,70 +0,0 @@ -# ICON4Py - ICON inspired code in Python and GT4Py -# -# Copyright (c) 2022, ETH Zurich and MeteoSwiss -# All rights reserved. -# -# This file is free software: you can redistribute it and/or modify it under -# the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or any later -# version. See the LICENSE.txt file at the top-level directory of this -# distribution for a copy of the license or check . -# -# SPDX-License-Identifier: GPL-3.0-or-later - -from functional.ffront.decorator import field_operator, program -from functional.ffront.fbuiltins import Field, FieldOffset, broadcast - -from icon4py.common.dimension import CellDim, KDim - - -Koff = FieldOffset("Koff", source=KDim, target=(KDim,)) - - -@field_operator -def _face_val_ppm_stencil_04_p_face_1( - p_cc: Field[[CellDim, KDim], float], - p_cellhgt_mc_now: Field[[CellDim, KDim], float], -) -> Field[[CellDim, KDim], float]: - - p_face = p_cc * ( - broadcast(1.0, (CellDim, KDim)) - - (p_cellhgt_mc_now / p_cellhgt_mc_now(Koff[-1])) - ) + (p_cellhgt_mc_now / (p_cellhgt_mc_now(Koff[-1]) + p_cellhgt_mc_now)) * ( - (p_cellhgt_mc_now / p_cellhgt_mc_now(Koff[-1])) * p_cc + p_cc(Koff[-1]) - ) - return p_face - - -@field_operator -def _face_val_ppm_stencil_04_p_face_2( - p_cc: Field[[CellDim, KDim], float], -) -> Field[[CellDim, KDim], float]: - p_face = p_cc - return p_face - - -@program -def face_val_ppm_stencil_04( - p_cc: Field[[CellDim, KDim], float], - p_cellhgt_mc_now: Field[[CellDim, KDim], float], - p_face: Field[[CellDim, KDim], float], - nudging: int, - halo: int, - vertical_lower: int, - vertical_upper: int, -): - _face_val_ppm_stencil_04_p_face_1( - p_cc, - p_cellhgt_mc_now, - out=p_face, - domain={CellDim: (nudging, halo), KDim: (1, 2)}, - ) - _face_val_ppm_stencil_04_p_face_1( - p_cc, - p_cellhgt_mc_now, - out=p_face, - domain={CellDim: (nudging, halo), KDim: (2, vertical_upper)}, - ) - _face_val_ppm_stencil_04_p_face_2( - p_cc, out=p_face, domain={CellDim: (nudging, halo), KDim: (vertical_lower, 1)} - ) diff --git a/advection/tests/test_face_val_ppm_stencil_04.py b/advection/tests/test_face_val_ppm_stencil_04.py deleted file mode 100644 index ae02704b3..000000000 --- a/advection/tests/test_face_val_ppm_stencil_04.py +++ /dev/null @@ -1,91 +0,0 @@ -# ICON4Py - ICON inspired code in Python and GT4Py -# -# Copyright (c) 2022, ETH Zurich and MeteoSwiss -# All rights reserved. -# -# This file is free software: you can redistribute it and/or modify it under -# the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or any later -# version. See the LICENSE.txt file at the top-level directory of this -# distribution for a copy of the license or check . -# -# SPDX-License-Identifier: GPL-3.0-or-later - -import numpy as np - -from icon4py.atm_dyn_iconam.face_val_ppm_stencil_04 import face_val_ppm_stencil_04 -from icon4py.common.dimension import CellDim, KDim -from icon4py.testutils.simple_mesh import SimpleMesh -from icon4py.testutils.utils import random_field, zero_field - - -def face_val_ppm_stencil_04_numpy( - p_cc: np.array, - p_cellhgt_mc_now: np.array, - nudging, - halo, - vertical_lower, - vertical_upper, -) -> np.array: - p_cellhgt_mc_now_minus_1 = np.roll(p_cellhgt_mc_now, shift=1, axis=1) - p_cc_minus_1 = np.roll(p_cc, shift=1, axis=1) - p_face_ref = np.zeros_like(p_cc) - - p_face_ref[nudging:halo, 1:vertical_upper] = p_cc[ - nudging:halo, 1:vertical_upper - ] * ( - 1 - - ( - p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] - / p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] - ) - ) + ( - p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] - / ( - p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] - + p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] - ) - ) * ( - ( - p_cellhgt_mc_now[nudging:halo, 1:vertical_upper] - / p_cellhgt_mc_now_minus_1[nudging:halo, 1:vertical_upper] - ) - * p_cc[nudging:halo, 1:vertical_upper] - + p_cc_minus_1[nudging:halo, 1:vertical_upper] - ) - - p_face_ref[nudging:halo, vertical_lower] = p_cc[nudging:halo, vertical_lower] - return p_face_ref - - -def test_face_val_ppm_stencil_04(): - mesh = SimpleMesh() - - p_cc = random_field(mesh, CellDim, KDim) - p_cellhgt_mc_now = random_field(mesh, CellDim, KDim) - p_face = zero_field(mesh, CellDim, KDim) - nudging = 4 - halo = 16 - vertical_lower = 0 - vertical_upper = np.asarray(p_face).shape[1] - - p_face_ref = face_val_ppm_stencil_04_numpy( - np.asarray(p_cc), - np.asarray(p_cellhgt_mc_now), - nudging, - halo, - vertical_lower, - vertical_upper, - ) - face_val_ppm_stencil_04( - p_cc, - p_cellhgt_mc_now, - p_face, - nudging, - halo, - vertical_lower, - vertical_upper, - offset_provider={"Koff": KDim}, - ) - - assert np.allclose(p_face_ref, p_face) From 1c8ca05e9000374e561814137aa76ca9dff00e93 Mon Sep 17 00:00:00 2001 From: abishekg7 Date: Fri, 14 Oct 2022 10:52:00 +0200 Subject: [PATCH 10/10] removing hflx_limiter_mo_stencil_01 for now, until part II is done --- .../advection/hflx_limiter_mo_stencil_01.py | 50 -------------- .../tests/test_hflx_limiter_mo_stencil_01.py | 67 ------------------- 2 files changed, 117 deletions(-) delete mode 100644 advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py delete mode 100644 advection/tests/test_hflx_limiter_mo_stencil_01.py diff --git a/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py b/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py deleted file mode 100644 index c3309e2ce..000000000 --- a/advection/src/icon4py/advection/hflx_limiter_mo_stencil_01.py +++ /dev/null @@ -1,50 +0,0 @@ -# ICON4Py - ICON inspired code in Python and GT4Py -# -# Copyright (c) 2022, ETH Zurich and MeteoSwiss -# All rights reserved. -# -# This file is free software: you can redistribute it and/or modify it under -# the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or any later -# version. See the LICENSE.txt file at the top-level directory of this -# distribution for a copy of the license or check . -# -# SPDX-License-Identifier: GPL-3.0-or-later - -from functional.ffront.decorator import field_operator, program -from functional.ffront.fbuiltins import Field, abs - -from icon4py.common.dimension import E2C, CellDim, EdgeDim, KDim - - -@field_operator -def _hflx_limiter_mo_stencil_01( - p_mflx_tracer_h: Field[[EdgeDim, KDim], float], - p_cc: Field[[CellDim, KDim], float], - p_mass_flx_e: Field[[EdgeDim, KDim], float], -) -> tuple[Field[[EdgeDim, KDim], float], Field[[EdgeDim, KDim], float]]: - - z_mflx_low = 0.5 * ( - p_mass_flx_e * (p_cc(E2C[0]) + p_cc(E2C[1])) - - abs(p_mass_flx_e) * (-p_cc(E2C[0]) + p_cc(E2C[1])) - ) - - z_anti = p_mflx_tracer_h - z_mflx_low - - return z_mflx_low, z_anti - - -@program -def hflx_limiter_mo_stencil_01( - p_mflx_tracer_h: Field[[EdgeDim, KDim], float], - p_cc: Field[[CellDim, KDim], float], - p_mass_flx_e: Field[[EdgeDim, KDim], float], - z_mflx_low: Field[[EdgeDim, KDim], float], - z_anti: Field[[EdgeDim, KDim], float], -): - _hflx_limiter_mo_stencil_01( - p_mflx_tracer_h, - p_cc, - p_mass_flx_e, - out=(z_mflx_low, z_anti), - ) diff --git a/advection/tests/test_hflx_limiter_mo_stencil_01.py b/advection/tests/test_hflx_limiter_mo_stencil_01.py deleted file mode 100644 index 5423f15e0..000000000 --- a/advection/tests/test_hflx_limiter_mo_stencil_01.py +++ /dev/null @@ -1,67 +0,0 @@ -# ICON4Py - ICON inspired code in Python and GT4Py -# -# Copyright (c) 2022, ETH Zurich and MeteoSwiss -# All rights reserved. -# -# This file is free software: you can redistribute it and/or modify it under -# the terms of the GNU General Public License as published by the -# Free Software Foundation, either version 3 of the License, or any later -# version. See the LICENSE.txt file at the top-level directory of this -# distribution for a copy of the license or check . -# -# SPDX-License-Identifier: GPL-3.0-or-later - -import numpy as np - -from icon4py.advection.hflx_limiter_mo_stencil_01 import ( - hflx_limiter_mo_stencil_01, -) -from icon4py.common.dimension import CellDim, EdgeDim, KDim -from icon4py.testutils.simple_mesh import SimpleMesh -from icon4py.testutils.utils import random_field - - -def hflx_limiter_mo_stencil_01_numpy( - e2c: np.array, - p_mflx_tracer_h: np.array, - p_cc: np.array, - p_mass_flx_e: np.array, -): - p_cc_e2c = p_cc[e2c] - z_mflx_low = 0.5 * ( - p_mass_flx_e * (p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) - - abs(p_mass_flx_e) * (-p_cc_e2c[:, 0] + p_cc_e2c[:, 1]) - ) - - z_anti = p_mflx_tracer_h - z_mflx_low - - return z_mflx_low, z_anti - - -def test_hflx_limiter_mo_stencil_01(): - mesh = SimpleMesh() - p_mflx_tracer_h = random_field(mesh, EdgeDim, KDim) - p_cc = random_field(mesh, CellDim, KDim) - p_mass_flx_e = random_field(mesh, EdgeDim, KDim) - z_mflx_low = random_field(mesh, EdgeDim, KDim) - z_anti = random_field(mesh, EdgeDim, KDim) - - z_mflx_low_ref, z_anti_ref = hflx_limiter_mo_stencil_01_numpy( - mesh.e2c, - np.asarray(p_mflx_tracer_h), - np.asarray(p_cc), - np.asarray(p_mass_flx_e), - ) - - hflx_limiter_mo_stencil_01( - p_mflx_tracer_h, - p_cc, - p_mass_flx_e, - z_mflx_low, - z_anti, - offset_provider={ - "E2C": mesh.get_e2c_offset_provider(), - }, - ) - assert np.allclose(z_mflx_low, z_mflx_low_ref) - assert np.allclose(z_anti, z_anti_ref)