Skip to content
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.

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

Initial condition for jw test #446

Merged
merged 29 commits into from
May 6, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
49e165a
Pull changes for initial conditions from Jablonowski_Williamson_test …
OngChia Apr 15, 2024
7c20ea6
Changes function names testcase_functions.py for better readability.
OngChia Apr 16, 2024
6e02762
Added initial conditions for JW.
OngChia Apr 17, 2024
7a900f7
Add back io_utils which later will become initialization_utils.py
OngChia Apr 17, 2024
8242587
remove existing initialization
OngChia Apr 17, 2024
4951b81
rename io_utils.py as initialization_utils.py
OngChia Apr 17, 2024
17e1d51
Add JW initial conditions in initialization_utils.py. Move non-compul…
OngChia Apr 17, 2024
8be6e12
Removed unecessary savepoints.
OngChia Apr 17, 2024
92c3db0
clean up. Remove unecessary codes that were used to output data for C…
OngChia Apr 18, 2024
994b32e
precommit
OngChia Apr 19, 2024
33abcd2
Merge branch 'main' into initial_condition_for_JW_test
OngChia Apr 22, 2024
37725c1
merge main
OngChia Apr 22, 2024
9445dd5
Remove dummy log file that was accidentally added
OngChia Apr 22, 2024
9ac856c
Removed unwanted diagnostic metric state variables that were used to …
OngChia Apr 23, 2024
ebd03db
merge main
OngChia Apr 23, 2024
617b312
add docstring to variables in horizontal.py
OngChia Apr 26, 2024
b092957
move referecen solution to second argument in allclose.
OngChia Apr 26, 2024
ccc4b05
create diagnostic_stencils and their tests
OngChia Apr 30, 2024
b44895d
merge main
OngChia Apr 30, 2024
b9a3061
Finalizing all stencil and data tests for diagnostic state variables,…
OngChia May 2, 2024
3342960
Removed init_exner_pr and temporary printing statement in test_diagno…
OngChia May 2, 2024
b557e1c
remove comment
OngChia May 2, 2024
21eb5a0
remove unsued rbf numpy stencil in testcase_functions.py
OngChia May 2, 2024
53a0d6c
remove bug in pressure stencil test
OngChia May 2, 2024
00ce346
Resolved pressure problem
OngChia May 2, 2024
4d16b5d
changes according to review
OngChia May 3, 2024
9fdc961
remove the pytest skip mark for test_jabw_initial_condition because j…
OngChia May 3, 2024
d3e1d7a
merge main
OngChia May 3, 2024
1908c52
skip test_diagnose_pressure when roundtrip backend
OngChia May 3, 2024
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
6 changes: 6 additions & 0 deletions model/common/src/icon4py/model/common/constants.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,8 @@
SPECIFIC_HEAT_CONSTANT_VOLUME: Final[wpfloat] = CPD - RD
CVD: Final[wpfloat] = SPECIFIC_HEAT_CONSTANT_VOLUME
CVD_O_RD: Final[wpfloat] = CVD / RD
RD_O_CPD: Final[wpfloat] = RD / CPD
CPD_O_RD: Final[wpfloat] = CPD / RD

#: Gas constant for water vapor [J/K/kg], rv in ICON.
GAS_CONSTANT_WATER_VAPOR: Final[wpfloat] = 461.51
Expand All @@ -38,6 +40,7 @@
#: Av. gravitational acceleration [m/s^2]
GRAVITATIONAL_ACCELERATION: Final[wpfloat] = 9.80665
GRAV: Final[wpfloat] = GRAVITATIONAL_ACCELERATION
GRAV_O_RD: Final[wpfloat] = GRAV / RD

#: reference pressure for Exner function [Pa]
REFERENCE_PRESSURE: Final[wpfloat] = 100000.0
Expand All @@ -50,6 +53,9 @@
# average earth radius in [m]
EARTH_RADIUS: Final[float] = 6.371229e6

#: Earth angular velocity [rad/s]
EARTH_ANGULAR_VELOCITY: Final[wpfloat] = 7.29212e-5

#: sea level temperature for reference atmosphere [K]
SEA_LEVEL_TEMPERATURE: Final[wpfloat] = 288.15
T0SL_BG: Final[wpfloat] = SEA_LEVEL_TEMPERATURE
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# 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 <https://www.gnu.org/licenses/>.
#
# SPDX-License-Identifier: GPL-3.0-or-later

Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
# 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 <https://www.gnu.org/licenses/>.
#
# SPDX-License-Identifier: GPL-3.0-or-later

from gt4py.next.common import GridType
from gt4py.next.ffront.decorator import field_operator, program, scan_operator
from gt4py.next.ffront.fbuiltins import Field, exp, int32, sqrt

from icon4py.model.common.dimension import CellDim, KDim
from icon4py.model.common.settings import backend
from icon4py.model.common.type_alias import vpfloat, wpfloat


@scan_operator(axis=KDim, forward=False, init=(0.0, 0.0, True))
def _scan_pressure(
state: tuple[vpfloat, vpfloat, bool],
ddqz_z_full: vpfloat,
temperature: vpfloat,
pressure_sfc: vpfloat,
grav_o_rd: wpfloat,
):
pressure_interface = (
pressure_sfc * exp(-grav_o_rd * ddqz_z_full / temperature)
if state[2]
else state[1] * exp(-grav_o_rd * ddqz_z_full / temperature)
)
pressure = (
sqrt(pressure_sfc * pressure_interface) if state[2] else sqrt(state[1] * pressure_interface)
)
return pressure, pressure_interface, False
Copy link
Contributor

Choose a reason for hiding this comment

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

So you return here two fields one with shape (cell, num_lev) and one with shape (cell, num_lev + 1) interesting for the horizontal stencils we sometimes had problems with that due to the bounds. But lets keep it that way as long as it works.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Good to know that. Probably, in this case, because the diagnose_pressure stencil should be called with bounds (:, 0:nlev), gt4py did not complain anything?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This stencil does not compile when backend is roundtrip. I added a few lines in test_diagnose_pressure.py to skip it for roundtrip backend.

Copy link
Contributor

Choose a reason for hiding this comment

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

there is a is_roundtrip function for that just in case.



@field_operator
def _diagnose_pressure(
ddqz_z_full: Field[[CellDim, KDim], wpfloat],
temperature: Field[[CellDim, KDim], vpfloat],
pressure_sfc: Field[[CellDim], vpfloat],
grav_o_rd: wpfloat,
) -> tuple[Field[[CellDim, KDim], vpfloat], Field[[CellDim, KDim], vpfloat]]:
pressure, pressure_ifc, _ = _scan_pressure(ddqz_z_full, temperature, pressure_sfc, grav_o_rd)
return pressure, pressure_ifc


@program(grid_type=GridType.UNSTRUCTURED, backend=backend)
def diagnose_pressure(
ddqz_z_full: Field[[CellDim, KDim], wpfloat],
temperature: Field[[CellDim, KDim], vpfloat],
pressure_sfc: Field[[CellDim], vpfloat],
pressure: Field[[CellDim, KDim], vpfloat],
pressure_ifc: Field[[CellDim, KDim], vpfloat],
grav_o_rd: wpfloat,
horizontal_start: int32,
horizontal_end: int32,
vertical_start: int32,
vertical_end: int32,
):
_diagnose_pressure(
ddqz_z_full,
temperature,
pressure_sfc,
grav_o_rd,
out=(pressure, pressure_ifc),
domain={CellDim: (horizontal_start, horizontal_end), KDim: (vertical_start, vertical_end)},
)
Original file line number Diff line number Diff line change
@@ -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 <https://www.gnu.org/licenses/>.
#
# SPDX-License-Identifier: GPL-3.0-or-later

from gt4py.next.common import GridType
from gt4py.next.ffront.decorator import field_operator, program
from gt4py.next.ffront.fbuiltins import Field, exp, int32, log

from icon4py.model.common.dimension import CellDim, KDim, Koff
from icon4py.model.common.settings import backend
from icon4py.model.common.type_alias import vpfloat, wpfloat


@field_operator
def _diagnose_surface_pressure(
exner: Field[[CellDim, KDim], vpfloat],
temperature: Field[[CellDim, KDim], vpfloat],
ddqz_z_full: Field[[CellDim, KDim], wpfloat],
cpd_o_rd: wpfloat,
p0ref: wpfloat,
grav_o_rd: wpfloat,
) -> Field[[CellDim, KDim], vpfloat]:
pressure_sfc = p0ref * exp(
cpd_o_rd * log(exner(Koff[-3]))
+ grav_o_rd
* (
ddqz_z_full(Koff[-1]) / temperature(Koff[-1])
+ ddqz_z_full(Koff[-2]) / temperature(Koff[-2])
+ 0.5 * ddqz_z_full(Koff[-3]) / temperature(Koff[-3])
)
)
return pressure_sfc


@program(grid_type=GridType.UNSTRUCTURED, backend=backend)
def diagnose_surface_pressure(
exner: Field[[CellDim, KDim], vpfloat],
temperature: Field[[CellDim, KDim], vpfloat],
ddqz_z_full: Field[[CellDim, KDim], wpfloat],
pressure_sfc: Field[[CellDim, KDim], vpfloat],
cpd_o_rd: wpfloat,
p0ref: wpfloat,
grav_o_rd: wpfloat,
horizontal_start: int32,
horizontal_end: int32,
vertical_start: int32,
vertical_end: int32,
):
_diagnose_surface_pressure(
exner,
temperature,
ddqz_z_full,
cpd_o_rd,
p0ref,
grav_o_rd,
out=pressure_sfc,
domain={
CellDim: (horizontal_start, horizontal_end),
KDim: (vertical_start, vertical_end),
},
)
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
# 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 <https://www.gnu.org/licenses/>.
#
# SPDX-License-Identifier: GPL-3.0-or-later

from gt4py.next.common import GridType
from gt4py.next.ffront.decorator import field_operator, program
from gt4py.next.ffront.fbuiltins import Field, int32

from icon4py.model.common.dimension import CellDim, KDim
from icon4py.model.common.settings import backend
from icon4py.model.common.type_alias import vpfloat


@field_operator
def _diagnose_temperature(
theta_v: Field[[CellDim, KDim], vpfloat],
exner: Field[[CellDim, KDim], vpfloat],
) -> Field[[CellDim, KDim], vpfloat]:
temperature = theta_v * exner
return temperature


@program(grid_type=GridType.UNSTRUCTURED, backend=backend)
def diagnose_temperature(
theta_v: Field[[CellDim, KDim], vpfloat],
exner: Field[[CellDim, KDim], vpfloat],
temperature: Field[[CellDim, KDim], vpfloat],
horizontal_start: int32,
horizontal_end: int32,
vertical_start: int32,
vertical_end: int32,
):
_diagnose_temperature(
theta_v,
exner,
out=(temperature),
domain={
CellDim: (horizontal_start, horizontal_end),
KDim: (vertical_start, vertical_end),
},
)
4 changes: 2 additions & 2 deletions model/common/src/icon4py/model/common/dimension.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@
CECDim = Dimension("CEC")
ECDim = Dimension("EC")
ECVDim = Dimension("ECV")
CECDim = Dimension("CEC")
CECECDim = Dimension("CECEC")
E2CDim = Dimension("E2C", DimensionKind.LOCAL)
E2VDim = Dimension("E2V", DimensionKind.LOCAL)
Expand All @@ -38,6 +37,7 @@
E2C2EODim = Dimension("E2C2EO", DimensionKind.LOCAL)
E2C2EDim = Dimension("E2C2E", DimensionKind.LOCAL)
C2E2CDim = Dimension("C2E2C", DimensionKind.LOCAL)
C2E2C2EDim = Dimension("C2E2C2E", DimensionKind.LOCAL)
C2E2C2E2CDim = Dimension("C2E2C2E2C", DimensionKind.LOCAL)
E2C = FieldOffset("E2C", source=CellDim, target=(EdgeDim, E2CDim))
C2E = FieldOffset("C2E", source=EdgeDim, target=(CellDim, C2EDim))
Expand All @@ -54,8 +54,8 @@
E2C2EO = FieldOffset("E2C2EO", source=EdgeDim, target=(EdgeDim, E2C2EODim))
E2C2E = FieldOffset("E2C2E", source=EdgeDim, target=(EdgeDim, E2C2EDim))
C2E2C = FieldOffset("C2E2C", source=CellDim, target=(CellDim, C2E2CDim))
C2E2C2E = FieldOffset("C2E2C2E", source=EdgeDim, target=(CellDim, C2E2C2EDim))
C2E2C2E2C = FieldOffset("C2E2C2E2C", source=CellDim, target=(CellDim, C2E2C2E2CDim))
C2CEC = FieldOffset("C2CEC", source=CECDim, target=(CellDim, C2E2CDim))
C2CECEC = FieldOffset("C2CECEC", source=CECECDim, target=(CellDim, C2E2C2E2CDim))

Koff = FieldOffset("Koff", source=KDim, target=(KDim,))
Expand Down
Loading
Loading