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

Read in ak, bk coefficients #36

Merged
merged 70 commits into from
Jan 30, 2024
Merged

Read in ak, bk coefficients #36

merged 70 commits into from
Jan 30, 2024

Conversation

mlee03
Copy link
Collaborator

@mlee03 mlee03 commented Nov 6, 2023

Purpose

The ak and bk coefficients required to compute the pressure levels are hard-coded in the eta module for km=71, km=72, km=91, and km=139 number of vertical levels. This limits the number of vertical levels a user can choose to run models. By allowing the users to provide the ak and bk values in a simple ASCII file, the users have the freedom to specify any number of vertical levels.

Code changes:

In this PR,

  1. The hard-coded ak and bk array values have been removed. Instead, these values are read in, for example, from the input/eta79.txt file where the input directory is expected to reside in the directory of the main program and 79 represents the value of km. Checks were put in place to ensure that

    • the size of the ak and bk arrays are equal to the specified value of km
    • the resulting values of eta and eta_vare monotonically increasing
  2. Existing unit tests have been modified to read in the coefficients from an eta file and new unit tests were created to test the modified eta module. These tests were modified under the assumption that in the future, all input files will reside in one main input directory where the input directory resides in the directory of the main program. This assumption required creating input directories to all the possible testing directories in tests.

  3. Lastly, the functions vertical_coordinate and compute_eta, originally located in the baroclinic_jablonowski_williamson module have been moved to a newly created util module in the pace/util/grid directory. This change was made because it seemed more fitting for the abovementioned functions to reside with other grid-related functions.

Requirements changes:

  • This PR requires an input directory with eta files named as eta##.txt where ## refers to the value of km

Checklist

Before submitting this PR, please make sure:

  • You have followed the coding standards guidelines established at Code Review Checklist.
  • Docstrings and type hints are added to new and updated routines, as appropriate
  • All relevant documentation has been updated or added (e.g. README, CONTRIBUTING docs)
  • For each public change and fix in pace-util, HISTORY has been updated
  • Unit tests are added or updated for non-stencil code changes

Additionally, if this PR contains code authored by new contributors:

  • The names of all the new contributors have been added to CONTRIBUTORS.md

@lharris4
Copy link

lharris4 commented Nov 6, 2023

Hi @mlee03 . Glad that you could clean up the predefined vertical levels, which is getting to be a real mess. I agree that the best place for creating the vertical levels is in the preprocessor.

One concern (not related to this commit) is that there is no foolproof way to automatically generate a good set of ak/bk levels, and when creating an entirely novel set of levels some degree of hand-tuning is necessary to avoid stability problems. Xi Chen had developed a tool which can help, and interpolating from predefined level sets is usually pretty good, but some user discretion is needed when making new level sets.

NB: In the future when I have broader comments like this, is there somewhere more appropriate to bring them up than in the comments of a semi-related PR? (Adding @bensonr @oelbert for their opinion.) Thanks.

@FlorianDeconinck
Copy link
Collaborator

@lharris4: There's no way to auto-generate stable levels, but is there a way to check that they are unstable? We could eventually code that in.

@lharris4
Copy link

lharris4 commented Nov 7, 2023

@FlorianDeconinck there is no a priori way to ensure stability but typically we can examine the differences in delp and delz in a variety of situations, which is what Xi's tool does. We could run tests to ensure that the differences do not exceed some tolerance. This might be a good way forward; furthermore smoothing can be applied to improve the smoothness and stability of the selected coordinate.

@FlorianDeconinck
Copy link
Collaborator

Logged for posterity: #38

On the other topic, though obviously I am but a stranger here, I shall remind that there's a "Discussion" system on Github, which is a wiki-liky structure and could accommodate ongoing discussions and/or thinking on feature sets.

@lharris4
Copy link

lharris4 commented Nov 7, 2023

@FlorianDeconinck thanks for logging that. Also the discussions could potentially be a great resource. There are a precious few discussions on the FV3 repo and I would like to see this functionality used more.

Copy link
Contributor

@oelbert oelbert left a comment

Choose a reason for hiding this comment

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

This is a good start but there are a few things still do do. I think we talked about moving from .txt to .netcdf format for the pressure coefficient files, but also it would be good to have the yaml config file set the location of the ak and bk input file, and then the driver can pass that location as necessary during initialization. That way we also don't need so many copies of the files in the repo as well.

Comment on lines 9 to 11
km = 79
eta_file = "./input/eta79.txt"
ak_79, bk_79 = np.loadtxt(eta_file, unpack=True)
Copy link
Contributor

Choose a reason for hiding this comment

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

I think this could go inside of test_set_hybrid_pressure_coefficients_correct

pressure_data = set_hybrid_pressure_coefficients(km=90)


def test_set_hybrid_pressure_coefficients_correct():
Copy link
Contributor

Choose a reason for hiding this comment

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

It might be worthwhile to parameterize this over the 79 and 91 k-level versions

raise ValueError("Unexpected ks value")

if ptop != pressure_data.ptop:
raise ValueError("Unexpected ptopt value")
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
raise ValueError("Unexpected ptopt value")
raise ValueError("Unexpected ptop value")

pressure_data = set_hybrid_pressure_coefficients(km)

if not np.array_equal(ak_79, pressure_data.ak):
raise ValueError("Unexpected values in ak array")
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
raise ValueError("Unexpected values in ak array")
raise ValueError("ak values do not match input file")

raise ValueError("Unexpected values in ak array")

if not np.array_equal(bk_79, pressure_data.bk):
raise ValueError("Unexpected values in bk array")
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
raise ValueError("Unexpected values in bk array")
raise ValueError("bk values do not match input file")

Comment on lines 4 to 5
eta_0 = 0.252
surface_pressure = 1.0e5 # units of (Pa), from Table VI of DCMIP2016
Copy link
Contributor

Choose a reason for hiding this comment

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

Since eta_0 is only used in vertical_coordinate maybe it should live there? I'm also starting to think that surface_pressure belongs in the util.constants at this point, could we move it there?

"""
eta = 0.5 * ((ak[:-1] + ak[1:]) / surface_pressure + bk[:-1] + bk[1:])
eta_v = vertical_coordinate(eta)
return eta, eta_v
Copy link
Contributor

Choose a reason for hiding this comment

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

Could these functions live in eta.py?

]
)
# set path where the eta file lives
GRID_DIR = os.path.join(os.path.abspath("./"), "input/")
Copy link
Contributor

Choose a reason for hiding this comment

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

Instead of hardcoding the ak and bk files to path/input/etaX.txt we should pass the location of the file through the yaml config, similar to how Frank is setting up the model to read grid files. I'm happy to chat with you about how to do this

Copy link
Contributor

Choose a reason for hiding this comment

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

Mi Kyung and I have been discussing this and plan to come up with a good location to place all necessary testing/example files.

Copy link
Contributor

Choose a reason for hiding this comment

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

In order to reduce the need for multiple versions of the same data, could the test refer to the file in input/?

Copy link
Contributor

Choose a reason for hiding this comment

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

Same comment as for tests/main/driver/input/eta79.txt

Copy link
Contributor

@fmalatino fmalatino Nov 9, 2023

Choose a reason for hiding this comment

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

Same comment as for tests/main/driver/input/eta79.txt

tests/main/fv3core/input/eta91.txt Outdated Show resolved Hide resolved
@mlee03
Copy link
Collaborator Author

mlee03 commented Jan 24, 2024

@oelbert @bensonr ready for another round of reviews 🤞

@bensonr
Copy link
Contributor

bensonr commented Jan 24, 2024

@mlee03 - can you or one of the other reviewers discuss the necessity of the empty init.py file in the tests/main/grid/ directory

@mlee03
Copy link
Collaborator Author

mlee03 commented Jan 24, 2024

@bensonr, good point, to my limited understanding, I don't think it's needed here since we're not creating a package and the test is not being referenced by any modules. I think I created it here because there was one is the tests/main/driver directory 😁 😁 😁 but I will remove it

@bensonr
Copy link
Contributor

bensonr commented Jan 24, 2024

@mlee03 - is there a procedure for creating the ak/bk netcdf files and is it documented anywhere?

@mlee03
Copy link
Collaborator Author

mlee03 commented Jan 24, 2024

I took the ak/bk coefficients from what was coded in the original version of the eta module and created a python script to generate the netcdf files. I can document how to use the python netcdf module or the xarray module to create these files. Where should it be documented?

@mlee03
Copy link
Collaborator Author

mlee03 commented Jan 25, 2024

Or a Jupyter notebook?

Copy link
Contributor

@oelbert oelbert left a comment

Choose a reason for hiding this comment

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

Looks good for the most part, but I still think the MetricTerms class should be dealing with any external files in the init method and not during potential runtime calls, so I've added some suggestions to the grid generation file.

Comment on lines 126 to 130
metric_terms = pace.util.grid.MetricTerms(
quantity_factory=quantity_factory, communicator=self.communicator
)
metric_terms._eta_file = namelist["grid_config"]["config"]["eta_file"]

Copy link
Contributor

Choose a reason for hiding this comment

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

Shouldn't the eta file go in the metric terms init?

):
self._grid_type = grid_type
self._dx_const = dx_const
self._dy_const = dy_const
self._deglat = deglat
self._halo = N_HALO_DEFAULT
self._eta_file = eta_file
Copy link
Contributor

Choose a reason for hiding this comment

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

The way this is currently structured it can cause reading the eta file at runtime, which I really want to avoid. I think a better way to do it is to set the private variables in the init method and just return them as cached properties:

Suggested change
self._eta_file = eta_file
(
self._ks,
self._ptop,
self._ak,
self._bk,
) = self._set_hybrid_pressure_coefficients(eta_file)

Comment on lines 592 to 600
if not self._eta_file == "None":
(
self._ks,
self._ptop,
self._ak,
self._bk,
) = self._set_hybrid_pressure_coefficients()
else:
raise ValueError("eta file is not specified")
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
if not self._eta_file == "None":
(
self._ks,
self._ptop,
self._ak,
self._bk,
) = self._set_hybrid_pressure_coefficients()
else:
raise ValueError("eta file is not specified")

@@ -587,12 +589,15 @@ def ks(self) -> util.Quantity:
number of levels where the vertical coordinate is purely pressure-based
"""
if self._ks is None:
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
if self._ks is None:

@@ -602,12 +607,15 @@ def ak(self) -> util.Quantity:
pk = ak + (bk * ps)
"""
if self._ak is None:
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
if self._ak is None:

@@ -631,12 +642,15 @@ def ptop(self) -> util.Quantity:
the pressure of the top of atmosphere level
"""
if self._ptop is None:
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
if self._ptop is None:

Comment on lines 645 to 653
if not self._eta_file == "None":
(
self._ks,
self._ptop,
self._ak,
self._bk,
) = self._set_hybrid_pressure_coefficients()
else:
raise ValueError("eta file is not specified")
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
if not self._eta_file == "None":
(
self._ks,
self._ptop,
self._ak,
self._bk,
) = self._set_hybrid_pressure_coefficients()
else:
raise ValueError("eta file is not specified")

@@ -2149,7 +2163,9 @@ def _set_hybrid_pressure_coefficients(self):
"",
dtype=Float,
)
pressure_coefficients = set_hybrid_pressure_coefficients(self._npz)
pressure_coefficients = eta.set_hybrid_pressure_coefficients(
self._npz, self._eta_file
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
self._npz, self._eta_file
self._npz, eta_file

plus adding eta_file as an explicit argument

Copy link
Contributor

Choose a reason for hiding this comment

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

I don't know if this needs to be a notebook but I don't think it has to not be a notebook ¯_(ツ)_/¯

Copy link
Contributor

Choose a reason for hiding this comment

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

Same as the other notebook ¯_(ツ)_/¯

Copy link
Contributor

@oelbert oelbert left a comment

Choose a reason for hiding this comment

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

Looks good!

Copy link
Contributor

@bensonr bensonr left a comment

Choose a reason for hiding this comment

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

ready to go on my end

@fmalatino fmalatino merged commit 34eeea4 into NOAA-GFDL:main Jan 30, 2024
2 checks passed
fmalatino added a commit that referenced this pull request Jan 30, 2024
* Benchmark changes

* Added benchmark configurations

* Removed benchmark configs from configs to be tested in unit tests

* Changed dace submodule to point to Florian fix

* Dace submodule points to gcc_dies_on_dacecpu branch

* Set 'rf_fast' to true in baroclinic_c384_cpu/gpu.yaml

* fix mpi4py version (#51)

* [Feature] Better translate test (#39) (#47)

Translate test: small improvements

- Parametrize perturbation upon failure
- Refactor error folder to be `pwd` based
- Fix GPU translate unable to dump error `nc`
- Fix mixmatch precision on translate test
- Update README.md

Test fix:
- Orchestrate YPPM for translate purposes

Misc:
- Fix bad logger formatting on DaCeProgress

* [NASA] [Feature] Guarding against unimplemented configuration (#40) (#48)

* [Feature] Guarding against unimplemented configuration (#40)

Guarding against unimplemented namelists options:
- a2b_ord4
- d_sw
- fv_dynamics
- fv_subgridz
- neg_adj3
- divergence damping
- xppm
- yppm

Misc:
- Fix `netcdf_monitor` not mkdir the directory
- Add `as_dict` to the dycore state to dump the dycore more easily

* Unused assert

* Update fv3core/pace/fv3core/stencils/yppm.py

Co-authored-by: Oliver Elbert <[email protected]>

* Update fv3core/pace/fv3core/stencils/xppm.py

Co-authored-by: Oliver Elbert <[email protected]>

* Change NotImplemented to ValueError for n_sponge<3

* lint

---------

Co-authored-by: Oliver Elbert <[email protected]>

* Re-try of updating dace submodule to track Florian fix branch

* Reverted gt4py submodule to match main checkout

* Added benchmark README file

* Read in ak, bk coefficients  (#36)

* initial changes to read in ak bk

* read ak/bk

* add xfail

* remove input dir

* further changes to unit tests

* finish up test

* add history

* commit uncommited files

* fix test comment

* add input to top

* read in data

* read in netcdf file in eta mod

* remove txt file

* test

* modify test and fix generate.py

* remove emacs backup file

* driver tests pass

* fix helper.py

* fix fv3core tests

* fix physics test

* fix grid tests

* nullcommconfig

* cleanup input

* remove driver input

* remove top level input

* fix circular import problems

* modify eta_file readin for test_restart_serial

* comment out 91 test

* rm safety checks

* revert diagnostics.py

* restore driver.py

* revert initialization.py

* restore state.py

* restore analytic_init.py

* restore init_utils.py and analytic_init.py

* restore c_sw.py

* d2a2c_vect.py

* restore fv3core/stensils

* restore translate_fvdynamics

* restore physics/stencils

* restore stencils

* remove circular dependency

* use pytest parametrize

* cleanup generation.py

* fstrinngs

* add eta_file to MetricTerm init

* remove eta_file argument in new_from_metric_terms and geos_wrapper

* use pytest parametrize for the xfail tests

* use pytest parametrize for the xfail tests

* fix geos_wrapper and grid

* fix tests

* fstring

* add test comments

* fix util/HISTORY.md

* fix comments

* remove __init__.py from tests/main/grid

* add jupyter notebooks to generate eta files

* generate ak,bk,ptop on metricterm init

* fix tests

* exploit np.all in eta mod

* remove tests/main/grid/input

* update ci

* test

* remove input

* edit ci yaml

* remove push

---------

Co-authored-by: mlee03 <[email protected]>

* Move Active Physics Schemes to Config (#44)

* initial commit, need to adapt and run tests

* revising scheme name

* tests pass

* update history

* linting

* changing typehints for physics schemes to enum instead of str

* driver now works with physics config enum, tests pass

* fixed tests

* missed one

* D-grid data input method (#42)

* Testing changes reflected across branches

* Undoing changes made in build_gaea_c5.sh

* Testing vscode functionality, by adding a change to external_grid branch

* Testing vscode functionality, by adding a change to external_grid branch

* Addition of from_generated method and calc_flag to util/pace/util/grid/generation.py

* Added get_grid method for external grid data to driver/pace/driver/grid.py

* Preliminary xarray netcdf read in method added to driver/pace/driver/grid.py

* Updating util/pace/util/grid/generation.py from_generated method

* Addition of external grid data read in methods for initialization of grid. Current method uses xarray to interact with netcdf tile files. Values for longitutde, latitude, number of points in x an y, grid edge distances read in.

* driver/examples/configs/test_external_C12_1x1.yaml

* Preliminary unit test for external grid data read in

* Current state of unit tests as of 27 Nov 2023

* External grid method and unit tests added

* Re-excluding external grid data yamls from test_example_configs.py

* Update driver/pace/driver/grid.py

Co-authored-by: Florian Deconinck <[email protected]>

* Changed name of grid initializer function to match NetCDF dependency and class descriptor

* Update util/pace/util/grid/generation.py

Moved position of doc string for "from_external" MetricTerms class method

Co-authored-by: Oliver Elbert <[email protected]>

* Fixed indentation error in generation.py from suggestion in PR 42

* Removal of TODO comment in grid.py, changes to method of file accessing in test_analytic_init, test_external_grid_*

* Changed grid data read-in unit tests to compare data directly from file to driver grid data generated from yaml

* Change to reading in lon and lat, other metric terms calculated as needed

* Removed read-in of dx, dy, and area. Changed unit tests to compare calculated area to 'ideal' surface area as given by selected constants type.

* Update tests/mpi_54rank/test_external_grid_1x1.py

Incorrect name of test in test_external_grid_1x1.py changed to match file name

Co-authored-by: Oliver Elbert <[email protected]>

* Added comparisons for read-in vs generated by driver lon, lat, dx, dy, and area data to unit tests

* Added relative error calculations to unit tests for external grid data read-in

* External grid data read in tests changed: relative errors printed by each rank and get_tile_number replacing get_tile_index

* Removing commented out sections in test_external_grid_2x2.py

Co-authored-by: Oliver Elbert <[email protected]>

* Updated external grid data read-in to take configuration and input data locations from command line, updated test description, and added documentation on grid construction to external grid data configuration selection dataclass.

* Updated documentation in grid.py

* Updated external grid data read in unit test to use parametrize functionality of pytest

* Ammended files to reference changes to PR 36

---------

Co-authored-by: Frank Malatino <[email protected]>
Co-authored-by: Florian Deconinck <[email protected]>
Co-authored-by: Oliver Elbert <[email protected]>

---------

Co-authored-by: Oliver Elbert <[email protected]>
Co-authored-by: Florian Deconinck <[email protected]>
Co-authored-by: Oliver Elbert <[email protected]>
Co-authored-by: MiKyung Lee <[email protected]>
Co-authored-by: mlee03 <[email protected]>
Co-authored-by: Frank Malatino <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants