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

+*Thickness-based viscosity and diffusivity units #409

Merged
merged 4 commits into from
Jul 30, 2023

Conversation

Hallberg-NOAA
Copy link
Member

This PR includes a series of 3 commits that revise the units of viscosity in the transparent vert_visc_type to those of a dynamic viscosity in non-Boussinesq more or a kinematic viscosity in Boussinsq mode, with analogous changes to other elements of the transparent vert_visc_type, and of a number of other diapycnal diffusivities and turbulent kinetic energies that are passed between routines or used internally within the routines that set these diffusivities. Dozens of unit conversion factors were also cancelled out with these changes, including using that GV%Z_to_H/GV%Rho_0 = GV%RZ_to_H and GV%Rho_0*GV%H_to_Z = GV%H_to_RZ for both Boussinesq or non-Boussinesq configurations.

Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are bitwise identical in that mode, but in non-Boussinesq mode this conversion involves multiplication and division by GV%Rho_0, so while all answers are mathematically equivalent, this change does change answers at roundoff in non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.

The commits in this PR include:

  • 48c278b14 +*Use [H Z2 T-3 ~> m3 s-3 or W m-2] for TKE units
  • 573f965d7 +Thickness-based diffusivity arguments
  • 3fc64c03e +*Revise the units of 12 vertvisc_type elements

@Hallberg-NOAA Hallberg-NOAA added enhancement New feature or request answer-changing A change in results (actual or potential) labels Jul 19, 2023
@codecov
Copy link

codecov bot commented Jul 19, 2023

Codecov Report

Merging #409 (4fd8ada) into dev/gfdl (84056b1) will decrease coverage by 0.01%.
The diff coverage is 37.34%.

❗ Current head 4fd8ada differs from pull request most recent head 288314a. Consider uploading reports for the commit 288314a to get more accurate results

@@             Coverage Diff              @@
##           dev/gfdl     #409      +/-   ##
============================================
- Coverage     38.15%   38.15%   -0.01%     
============================================
  Files           269      269              
  Lines         76786    76784       -2     
  Branches      14131    14130       -1     
============================================
- Hits          29301    29298       -3     
- Misses        42201    42202       +1     
  Partials       5284     5284              
Files Changed Coverage Δ
src/core/MOM_variables.F90 46.96% <ø> (ø)
src/parameterizations/vertical/MOM_CVMix_KPP.F90 0.77% <0.00%> (ø)
src/parameterizations/vertical/MOM_CVMix_conv.F90 8.13% <0.00%> (ø)
src/parameterizations/vertical/MOM_CVMix_ddiff.F90 8.75% <0.00%> (ø)
src/parameterizations/vertical/MOM_CVMix_shear.F90 11.11% <0.00%> (ø)
...rc/parameterizations/vertical/MOM_diabatic_aux.F90 59.65% <0.00%> (ø)
...rc/parameterizations/vertical/MOM_tidal_mixing.F90 1.40% <0.00%> (+<0.01%) ⬆️
src/user/user_change_diffusivity.F90 0.00% <0.00%> (ø)
...c/parameterizations/vertical/MOM_set_viscosity.F90 60.98% <38.46%> (+0.06%) ⬆️
...rameterizations/vertical/MOM_entrain_diffusive.F90 68.96% <40.00%> (ø)
... and 7 more

... and 1 file with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@marshallward
Copy link
Member

The introduction of GV%Z_to_H into many of the solvers could noticeably affect the runtime. Will these get removed in later pull requests? If not, has anyone looked at timings after this change?

If the time increase is significant, is there any way to pre-compute these H-to-Z corrections? Or must they happen in the expressions?

@Hallberg-NOAA
Copy link
Member Author

This PR is the high-water mark for these GV%Z_to_H conversion factors. They almost all disappear in the follow-on commits with the refactoring and non-Boussinesq translation of the various parts of the code. They are here because I thought that this PR, which deals with the necessary coordinated changes across modules, would be challenging enough to review as is without rolling in even more extensive changes!

@marshallward
Copy link
Member

This PR is the high-water mark for these GV%Z_to_H conversion factors. They almost all disappear in the follow-on commits with the refactoring and non-Boussinesq translation of the various parts of the code. They are here because I thought that this PR, which deals with the necessary coordinated changes across modules, would be challenging enough to review as is without rolling in even more extensive changes!

OK sounds good, I figured they might get unrolled. I'll merge this once the conflicts are resolved (and pending testing).

@marshallward
Copy link
Member

marshallward commented Jul 28, 2023

There seem to be some OpenMP errors with this:

/lustre/f2/scratch/oar.gfdl.ogrp-account/runner/171/work/tc4/symmetric/ocean.stats /lustre/f2/scratch/oar.gfdl.ogrp-account/runner/171/work/tc4/openmp/ocean.stats differ: byte 619, line 4
real	0m12.295s
user	0m0.291s
sys	0m1.138s
2797,2814c2797,2814
< h-point: mean=   2.1147105230955363E-20 min=  -5.9211894646675019E-18 max=   5.9211894646675019E-18 ocean_model-boundary_forcing_saln_tendency
< h-point: c=       800 ocean_model-boundary_forcing_saln_tendency
<  column: mean=   2.1148451835789589E-20 min=  -4.2293921285443795E-20 max=   8.4590824957022972E-20 ocean_model-boundary_forcing_saln_tendency_xyave
<  column: c=        72 ocean_model-boundary_forcing_saln_tendency_xyave
< h-point: mean=   2.1148293340564746E-20 min=  -5.9211894646675019E-18 max=   5.9211894646675019E-18 ocean_model_z-boundary_forcing_saln_tendency
< h-point: c=      1194 ocean_model_z-boundary_forcing_saln_tendency
<  column: mean=   2.1148451835789580E-20 min=  -4.2291776600114608E-20 max=   8.4588680271693768E-20 ocean_model_z-boundary_forcing_saln_tendency_xyave
<  column: c=        78 ocean_model_z-boundary_forcing_saln_tendency_xyave
< h-point: mean=   7.0039212524924156E-19 min=  -1.9610979506978765E-16 max=   1.9610979506978765E-16 ocean_model-boundary_forcing_salt_tendency
< h-point: c=       775 ocean_model-boundary_forcing_salt_tendency
<  column: mean=   7.0039908052431278E-19 min=  -1.4007812635693276E-18 max=   2.8015794246179531E-18 ocean_model-boundary_forcing_salt_tendency_xyave
<  column: c=        58 ocean_model-boundary_forcing_salt_tendency_xyave
< h-point: mean=   7.0039212524924156E-19 min=  -1.9610979506978765E-16 max=   1.9610979506978765E-16 ocean_model_z-boundary_forcing_salt_tendency
< h-point: c=       775 ocean_model_z-boundary_forcing_salt_tendency
<  column: mean=   7.0039908052431278E-19 min=  -1.4007812635693276E-18 max=   2.8015794246179531E-18 ocean_model_z-boundary_forcing_salt_tendency_xyave
<  column: c=        58 ocean_model_z-boundary_forcing_salt_tendency_xyave
< h-point: mean=   1.4007842504984831E-18 min=  -1.9610979506978765E-16 max=   3.9221959013957530E-16 ocean_model-boundary_forcing_salt_tendency_2d
< h-point: c=       745 ocean_model-boundary_forcing_salt_tendency_2d
---
FAIL: Diagnostics tc4.openmp.diag have changed.
make: *** [Makefile:468: tc4.openmp.diag] Error 1
make: *** Waiting for unfinished jobs....
DONE: tc4.dim.h; no runtime errors.
PASS: Diagnostics tc4.nan.diag agree.
4,5c4,5
<      9,       0.125,     0, En 3.4939057301934371E-08, CFL  0.00016, SL  4.2633E-14, M 7.18210E+13, S 35.0000, T  0.3271, Me -1.57E-17, Se -6.74E-15, Te -4.35E-16
<     18,       0.250,     0, En 9.1595372175907653E-08, CFL  0.00026, SL  4.2633E-14, M 7.18210E+13, S 35.0000, T  0.3271, Me -4.25E-19, Se -5.60E-15, Te -4.22E-16
---
>      9,       0.125,     0, En 3.4939057301934371E-08, CFL  0.00016, SL  4.2633E-14, M 7.18210E+13, S 35.0000, T  0.3271, Me -1.57E-17, Se  6.01E-08, Te -4.35E-16
>     18,       0.250,     0, En 9.1595372175907653E-08, CFL  0.00026, SL  4.2633E-14, M 7.18210E+13, S 35.0000, T  0.3271, Me -4.25E-19, Se  7.47E-08, Te -4.22E-16
FAIL: Solutions tc4.openmp have changed.

I can send full output if needed.

BTW this only seems to be happening with Intel, not GNU.

@Hallberg-NOAA
Copy link
Member Author

The openMP bug noted in the pipeline failure above actually came from an intermittent failure introduced with the last PR to be handled (PR #401). This bug will be addressed by PR #427.

  Revised the units of 12 vertvisc_type elements to be based on thicknesses, so
that vertical viscosities (in [H Z T-1 ~> m2 s-1 or Pa s]) are stored as dynamic
viscosites when in non-Boussinesq mode, with analogous changes to the diapycanl
diffusivity (now in [H Z T-1 ~> m2 s-1 or kg m-1 s-1]).  Similarly changed the
units of the 2 Rayleigh drag velocity elements (Ray_u and Ray_v) of the
vertvisc_type from vertical velocity units to thickness flux units and to more
accurately reflect the nature of these fields.  The bottom boundary layer TKE
source element (TKE_BBL) was also revised to [H Z2 T-3 ~> m3 s-3 or W m-2].
This commit also adds required changes to the units of the viscosities or
shear-driven diffusivities returned from KPP_calculate, calculate_CVMix_shear,
calculate_CVMix_conv, Calculate_kappa_shear, Calc_kappa_shear_vertex,
calculate_tidal_mixing and calculate_CVMix_tidal.

 Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
  Rescaled diapycnal diffusivities passed as arguments in non-Boussinesq mode,
to be equivalent to the thermal conductivity divided by the heat capacity,
analogously to the difference between a kinematic viscosity and a dynamic
viscosity, so that the new units are [H Z T-1 ~> m2 s-1 or kg m-1 s-1].

  This includes changing the units of 4 arguments to set_diffusivity;
3 arguments each to calculate_bkgnd_mixing, add_drag_diffusivity,
add_LOTW_BBL_diffusivity, user_change_diff, calculate_tidal_mixing and
add_int_tide_diffusivity; 2 arguments to KPP_calculate, calculate_CVMix_conv,
compute_ddiff_coeffs, differential_diffuse_T_S, entrainment_diffusive,
double_diffusion, add_MLrad_diffusivity, and calculate_CVMix_tidal; and one
argument to energetic_PBL.

  The units of 36 internal variables were also changed, as were a total of
29 elements in various opaque types, including 8 elements in bkgnd_mixing_cs,
2 in diabatic_CC, 3 in tidal_mixing_diags type, 1 in user_change_diff_CS,
9 in set_diffusivity_CS type, and 6 elements in diffusivity_diags.

  Two new internal variables were added, and several redundant GV%H_to_Z
conversion factors were also cancelled out, some using that
GV%H_to_Z*GV%Rho0 = GV%H_to_RZ.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
  Changed the units for TKE arguments to [H Z2 T-3 ~> m3 s-3 or W m-2] for
find_TKE_to_Kd, add_drag_diffusivity, calculate_tidal_mixing and
add_int_tide_diffusivity, with similar changes to the units of 21 diagnostics
or internal variables in the same routines and in add_LOTW_BBL_diffusivity and
add_MLrad_diffusivity.  Dozens of unit conversion factors were also cancelled
out with these changes, including using that GV%Z_to_H/GV%Rho_0 = GV%RZ_to_H and
that GV%Rho_0*GV%H_to_Z = GV%H_to_RZ for both Boussinesq or non-Boussinesq
configurations.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode, but in non-Boussinesq mode this conversion
involves multiplication and division by GV%Rho_0, so while all answers are
mathematically equivalent, this change does change answers at roundoff in
non-Boussinesq mode unless GV%Rho_0 is chosen to be an integer power of 2.
@marshallward
Copy link
Member

@marshallward marshallward merged commit 359bdcb into NOAA-GFDL:dev/gfdl Jul 30, 2023
@Hallberg-NOAA Hallberg-NOAA deleted the visc_type_units branch May 24, 2024 07:41
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
answer-changing A change in results (actual or potential) enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants