added nonhydrostatic capabilities to MPAS-O #21

scalandr · 2022-05-20T21:17:59Z

This is a nonhydrostatic version of MPAS-O. This nonhydrostatic model enables new physics, allowing the correct representation of, for example, high-resolution internal waves dynamics, which are not well represented by the hydrostatic approximation. New vertical momentum and pressure terms are added in the formulation, and a coupled elliptic system is formed for the solution of the nonhydrostatic pressure correction. At every time step, this system is solved with the PETSc library using preconditioned iterative Krylov solvers. The flag config_enable_nonhydrostatic_mode can be used to turn on and off the nonhydrostatic capabilities.

xylar

@scalandr, thank you for this amazing work! This will be such an great new capability for us to have.

I have noticed several things that need a little work. I still have a lot of code to go through but this is all I have time for for a little while so wanted to give you a chance to address it first.

The major one is that most if not all of the variables you added should be in the nonhydrostaticPKG package (I suggest this name change). In the code, it will be important to only access these fields if the package is active. Typically, you can check if the field pointers are associated. See many examples of this throughout the code. If you don't make use of packages, your changes here will add substantially to the memory usage of the model even in hydrostatic mode.

The other, which is preventing the code from compiling and therefore preventing me from verifying that the results in hydrostatic mode are bit-for-bit is the incorrect indentation in src/shared/Makefile. It needs to be with tabs, not spaces.

Please verify that the code compiles for you and then ping me to try testing again.

xylar · 2022-05-21T22:30:16Z

components/mpas-framework/Makefile

-	 ($(FC) pio2.f90 $(FCINCLUDES) $(FFLAGS) $(LDFLAGS) $(LIBS) -o pio2.out &> /dev/null && echo "=> PIO 2 detected") || \
+	 ($(FC) pio2.f90 $(FCINCLUDES) $(FFLAGS) $(LDFLAGS) $(LIBS) -o pio2.out && echo "=> PIO 2 detected") || \


@scalandr, this line needs to go back to what it used to be. I presume this was just to help you debug.

@xylar, done.

components/mpas-ocean/src/Registry.xml

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

xylar · 2022-05-22T18:11:40Z

components/mpas-ocean/src/Registry.xml

+                <var name="dynamicViscosity" type="real" dimensions="nVertLevelsP1 nCells Time" units="Ns/m^2"
+                         description="dynamic viscosity of water as a function of temperature and salinity"
+                         packages="forwardMode;analysisMode"
+                />
+                <var name="uStar" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1}"
+                  	 description="stress velocity for wall bounded flows"
+                  	 packages="forwardMode;analysisMode"
+    		/>
+    		<var name="yPlus" type="real" dimensions="nVertLevels nCells Time" units="m"
+       			 description="boundary normal coordinate"
+      		         packages="forwardMode;analysisMode"
+    		/>


Add to nonhydrostaticPKG

xylar · 2022-05-22T18:12:51Z

components/mpas-ocean/src/Registry.xml

+                <var name="ugrad" type="real" dimensions="R3 nVertLevels nCells Time" units="s^{-1}"
+                        description="three-dimensional gradient of zonal velocity"
+                />
+                <var name="vgrad" type="real" dimensions="R3 nVertLevels nCells Time" units="s^{-1}"
+                        description="three-dimensional gradient of meridional velocity"
+                />
+                <var name="wgrad" type="real" dimensions="R3 nVertLevels nCells Time" units="s^{-1}"
+                        description="three-dimensional gradient of vertical velocity"
+                />
+                <var name="tgrad" type="real" dimensions="R3 nVertLevels nCells Time" units="s^{-1}"
+                        description="three-dimensional gradient of temperature"
+                />
+                <var name="sgrad" type="real" dimensions="R3 nVertLevels nCells Time" units="s^{-1}"
+                        description="three-dimensional gradient of salinity"
+                />
+                <var name="Sh" type="real" dimensions="nVertLevels nCells Time" units="unitless"
+                        description="Structure function for heat and salt -- may want to consider separating at some point"
+                />
+                <var name="Sm" type="real" dimensions="nVertLevels nCells Time" units="unitless"
+                        description="Structure function for momentum"
+                />


Add to nonhydrostaticPKG

xylar · 2022-05-22T18:13:21Z

components/mpas-ocean/src/Registry.xml

+                <var name="Sh" type="real" dimensions="nVertLevels nCells Time" units="unitless"
+                        description="Structure function for heat and salt -- may want to consider separating at some point"
+                />
+                <var name="Sm" type="real" dimensions="nVertLevels nCells Time" units="unitless"
+                        description="Structure function for momentum"
+                />


I'm worried that the names Sh and Sm are not sufficiently descriptive or unique.

xylar · 2022-05-22T18:23:18Z

components/mpas-ocean/src/shared/Makefile

+		 mpas_ocn_vvel_vert_advection.o \
+		 mpas_ocn_vvel_horiz_advection.o \
+		 mpas_ocn_vvel_coriolis.o \
+		 mpas_ocn_vvel_hmix_del2.o \
+		 mpas_ocn_vvel_advection.o \
+		 mpas_ocn_vvel_pgf.o \
+		 mpas_ocn_nonhydrostatic_pressure_solve.o \


The indentation here seems to be incorrect.

xylar · 2022-05-22T18:24:03Z

components/mpas-ocean/src/shared/Makefile

+           mpas_ocn_two_equation_model_driver.o \
+           mpas_ocn_two_equation_structure_functions.o \
+           mpas_ocn_two_equation_model_driver.o \
+           mpas_ocn_tke_glsPsi_transport.o \
+           mpas_ocn_tke_advection_mono.o \
+           mpas_ocn_glsPsi_advection_mono.o \
+           mpas_ocn_shear_production.o \
+           mpas_ocn_buoyancy_production.o


Indentation here needs to be with tabs, not spaces. As a result, I am not able to build the code on Anvil.

scalandr · 2022-07-01T18:48:29Z

Hi @xylar, I am ready to ask for a re-review of the PR. Before I click on the "re-request review" button I wanted to ask you about the two conflicting files. If I click on "Resolve conflicts" and fix them there, then my local branch will not see the changes, correct? What is your recommendation on how to proceed? Thanks!

xylar · 2022-07-01T19:36:14Z

@scalandr, my recommendation for resolving conflicts like this is always to rebase onto E3SM-Project/E3SM/master.

If you use "Resolve conflicts" here, it will do a merge of E3SM-Ocean-Discussion/E3SM/master to this branch, but that's not very helpful since that isn't the same thing as E3SM-Project/E3SM/master and some merge conflicts might remain. Also, I am a strong advocate of avoiding merging at nearly all cost expect for merging pull requests. If you merge E3SM-Project/E3SM/master into your branch instead of rebasing, things can become very messy to untangle later. (Of course, if mistakes happen during a rebase, that can also be tricky to untangle later, but I still very much think it's worth it.)

I will update E3SM-Ocean-Discussion/E3SM/master to match E3SM-Project/E3SM/master right now. That way, if you rebase onto E3SM-Project/E3SM/master, you won't see any extra commits that aren't yours.

xylar · 2022-07-01T19:38:28Z

Okay, I just updated E3SM-Ocean-Discussion/E3SM/master. This shows why it's important to use E3SM-Project/E3SM/master as the base branch instead: there are now 4 conflicting files, not jut 2.

scalandr · 2022-07-01T19:58:17Z

@xylar, when we did the rebase together last time we did

git remote add E3SM-Project/E3SM [email protected]:E3SM-Project/E3SM.git
git fetch --all -p
git rebase -i E3SM-Project/E3SM/master

Shall i just do git fetch --all -p before the git rebase command this time?

xylar · 2022-07-01T20:12:06Z

@scalandr, yes, that's right. You don't need to add E3SM-Porject/E3SM as a remote since it's already been added.

xylar

I see a few basically cosmetic things that I'd like to see fixed. After that, I'm ready for this to move to E3SM.

I'd like to have @dengwirda and @mark-petersen both approve here first, though.

xylar · 2022-07-05T12:55:56Z

components/mpas-ocean/src/Registry.xml

+                <nml_option name="config_use_vertMom_del2" type="logical" default_value=".false." units="unitless"
+                                        description="If true, Laplacian horizontal mixing is used on the vertical momentum equation."


The registry needs to be indented with tabs, not spaces everywhere.

xylar · 2022-07-05T12:56:50Z

components/mpas-ocean/src/Registry.xml

+        <nml_record name="nonhydrostatic_pressure">
+            <nml_option name="config_enable_nonhydrostatic_mode" type="logical" default_value=".false." units="unitless"
+                                        description="Flag to enable the nonhydrostatic pressure calculation"
+                                        possible_values=".true. or .false."
+                />
+            <nml_option name="config_nonhydrostatic_solve_surface_boundary_condition" type="character" default_value="noGradient" units="unitless"
+                                    description="setting for surface boundary condition for the pressure correction solve"
+                                        possible_values="noGradient or no pressure"
+                />
+                <nml_option name="config_enable_vertMomentum_disable_nonHydro" type="logical" default_value=".false." units="unitless"
+                                        description="flag to enable hydrostatic pressure with vertical momentum tendencies"
+                                        possible_values=".true. or .false."
+                />


In the process of replacing spaces with tabs for indentation, please make the indentation consistent

xylar · 2022-07-05T13:00:18Z

components/mpas-ocean/src/Registry.xml

+                <var name="vertVelocityNonhydro" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^{-1}"
+			 description="prognostic vertical velocity, only used when config_enable_nonhydrostatic_mode = .true."
+    		/>
+                <var name="nonhydrostaticPressure" type="real" dimensions="nVertLevels nCells Time" units="N m^{-2}"
+			 description="nonhydrostatic pressure, set to zero if config_enable_nonhydrostatic_mode = .false."
+                />


Can these variables be put in the nonhydrostaticPKG package? If they aren't used at all when the nonhydrostaticPKG is disabled, this should be possible. If not, this will add to the memory used in hydrostatic mode.

xylar · 2022-07-05T13:00:38Z

components/mpas-ocean/src/Registry.xml

+                <var name="tendVertVelocityNonhydro" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^{-2}" name_in_code="vertVelocityNonhydroTend"
+                         description="time tendency of the vertical velocity (Only active for non hydrostatic mode)"
+                />


Same here, can this be in the nonhydrostaticPKG package?

xylar · 2022-07-05T13:01:30Z

components/mpas-ocean/src/Registry.xml

+                <var name="velocityHmixTend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocitydue to horizontal mixing"
+                />
+                <var name="velocityVmixTend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocity due to vertical mixing"
+                />
+                <var name="velocityCoriolisTend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocity due to Coriolis"
+                />
+                <var name="velocityPGFtend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocity due to PGF"
+                />
+                <var name="velocityVadvTend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocity due to vertical advection"
+                />
+                <var name="velocityForcingTend" dimensions="nVertLevels nEdges Time" units="m/s^2" type="real"
+                  description="tendency of normal velocity due to top and bottom forcing"
+                />


I believe this still needs to happen.

xylar · 2022-07-05T13:03:02Z

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

+         print*, dt
+


Please take this out. I assume it's for debugging.

xylar · 2022-07-05T13:43:31Z

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

-
           call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err )
           block => block % next
        end do
-
        call mpas_timer_stop("RK4 vel/thick tendency computations")


Maybe put back these blank lines?

xylar · 2022-07-05T13:43:58Z

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

        call mpas_timer_stop("RK4 vel/thick tendency computations")

-        ! Update halos for prognostic variables.


Maybe put back this comment (below the if)

xylar · 2022-07-05T13:44:29Z

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

+            call ocn_nonhydrostatic_pressure_update_velocity(normalVelocityNew, vertVelocityNonhydroNew, &
+                 layerThicknessNew, layerThickEdgeMean, nonhydrostaticPressure, nonhydrostaticPressureOld, &
+                 nhPressureCorrection, dt)


This is indented 3 spaces too many

xylar · 2022-07-05T13:45:44Z

components/mpas-ocean/src/shared/mpas_ocn_mesh.F

@@ -658,7 +658,7 @@ subroutine ocn_meshCreate(domain) !{{{
         !$acc                   edgeNormalVectors,        &
         !$acc                   localVerticalUnitVectors, &
         !$acc                   cellTangentPlane,         &
-         !$acc                   coeffs_reconstruct)
+         !$acc                   coeffs_reconstruct)       


trailing space should be removed.

xylar · 2022-07-05T13:49:12Z

It looks like one more conflicting commit came in in the last few days. @scalandr, you'll need to rebase one more time.

I ran the nightly test suite from compass with and without PETSc (and NetLib LAPACK). The first one was fine. The second was fine except for tests with different thread decompositions, as we have already discussed. I didn't try rerunning @scalandr's new nonhydro test case yet but will do that as part of an E3SM review.

…uns fine

scalandr · 2022-07-13T02:29:54Z

Hi @vanroekel, I took care of all your comments, so

I got rid of unused flags like config_use_constant_forced_pgf, config_disable_ALE, and config_enable_vertMomentum_disable_nonHydro,
I renamed nonhydrostaticPressure and nonhydrostaticPressureOld to nonhydrostaticPressureNew and nonhydrostaticPressureCur, respectively, in all time-steppers,
I removed viscosity from hmix tend,
I cleaned other files you mentioned.
About the vertical mixing in the vertical momentum, I believe I should inline the tri-diagonal solver. It’s something that we’ll have to do eventually, so let’s just do it now. Beside the tri-diagonal solver, the only thing left is the rebase.

scalandr · 2022-07-13T02:33:07Z

@xylar, If I do the rebase I see

So should I erase the first 7 lines, i.e. commits 88fb to 7355?

xylar · 2022-07-13T05:12:46Z

@scalandr, no, sorry for the confusion. Here's what we're seeing on GitHub:

Those are the commits that we would want to get rid of. But I now think it's just GitHub glitching and not noticing that those commits aren't part of this branch. We've seen that before. I think if I change the target branch to something else and then change it back to E3SM-Ocean-Discussion:master, it should be fine. I'll try that.

xylar · 2022-07-13T05:14:42Z

Indeed, that seems to have taken care of it:

@vanroekel, make sure you agree that the commits that shouldn't have been here are gone.

vanroekel · 2022-07-13T15:06:53Z

@scalandr the new version is looking great! I think there is just one more question about the tridiagonal solve in my mind. It looks like you inlined that code. If that is true, can you remove the code for the solver itself? it shouldn't be needed anymore

scalandr · 2022-07-14T01:46:26Z

Hi @vanroekel, I finished inlining the tridiagonal solve in ocn_vertVel_vmix_tend_implicit and got rid of the subroutine tridiagonal_solve. I also corrected the module description in the horizontal mixing module for the vertical momentum, that is ocn_vvel_hmix_del2.

vanroekel

This is looking great! I'm ready for this to move on to e3sm.

This commit resolves reveiew comments for PR #21 ('Add a regional sea-level projection capability to MALI') - Set up MPI variables outside of routine 'interpolation_init'. - Move reformatting command of 1D SL array to 2D array inside of routine 'interpolation' - Include #ifdef USE_SEALEVELMODEL inside of the check for config_uplift_method - other small fixes

…as_ocn_vel_pressure_grad.F

xylar · 2023-03-10T20:07:38Z

@scalandr, I noticed when building this today that the make targets are out of date compared to master. Could you do a rebase to bring those changes in when you have time? Not super urgent but it would be a good idea.

scalandr · 2023-03-13T18:22:56Z

Hi @xylar, as far as I know @jonbob already rebased this nonhydro branch locally. What would be the best way to proceed? For him to push the rebased version? Thanks!

xylar · 2023-03-13T23:08:11Z

@scalandr and @jonbob, that's something for the 2 of you to coordinate but I think that sounds like the easiest solution. In general, it is very important that no one but the developer force-pushes to a branch on a fork without careful coordination.

…d call

…nhydro

…Dev/develop Previously, when the regional sea-level prediction capability was added to MALI (#21), the restart config option for the sea-level model was not added. This led the sea-level model to get initialized to Timestep zero when coupled MALI-SLM simulations are being restarted, forgetting about the ice loading changes and associated viscoelastic solid earth deformation that happened in the timesteps prior to current model time. This PR fixes the problem by allowing the sea-level model to resume where it was left off. Note in parallel to this PR, the version of the SLM needs to incorporate the changes made in the following accompanying PR (MALI-Dev/1DSeaLevelModel_FWTW#9) * hollyhan/add_restart_functionality_slm: Don't call SLM on init of a restart Add addl info on restart about the calculated time since last SLM call Allow restarts at any interval when using SLM Add missing error flag so model actually dies when error occurs Add missing arguments to log write statement Update restart check to also use time interval division Adjust check if adaptive dt is on or not Update checks using interval division Improve error handling, correct other usage of config_uplift_method Improve synchronization of timesteps between MALI and SLM Add restart option when the SLM is coupled to MALI

scalandr requested review from xylar, dengwirda and vanroekel May 20, 2022 21:18

xylar requested changes May 22, 2022

View reviewed changes

scalandr force-pushed the ocean/nonhydro branch from 38a4930 to 59638ba Compare July 1, 2022 21:07

scalandr requested a review from xylar July 2, 2022 01:58

xylar reviewed Jul 5, 2022

View reviewed changes

scalandr added 17 commits July 6, 2022 19:25

added nonhydro to v2

88fb291

internal wave and overflow run fine

4c0eeaa

added seven modules for turbulence

28d6588

modified files in shared folder for les turbulence

e84e755

added les turbulence changes in the time-stepping schemes, the code r…

91dff6b

…uns fine

correction for RK4

dba0784

removed a couple of print statements

7355644

changes in mpas_ocn_core_interface.F

5ca480c

compiles with OpenMP

a5fd43f

got rid of two equation model (for now)

95295fa

addressed comments 1-5

046ce71

changed name to nonhydrostaticPKG and added variables to package

650c6dd

changed name from verticalVelocity to vertVelocityNonhydro

c2b2412

cleaned nonhydrostatic_pressure_solve.F and vel_hadv_coriolis.F

3eacc9e

cleaned subroutine ocn_vert_transport_velocity_top

c9cf963

cleaned mpas_ocn_vel_pressure_grad.F

e81723c

now it compiles, fixed src/shared/Makefile

2c5bcbf

xylar changed the base branch from master to alternate July 13, 2022 05:13

xylar changed the base branch from alternate to master July 13, 2022 05:13

inlined implicit vert mix in vert momentum

a111c41

vanroekel approved these changes Jul 14, 2022

View reviewed changes

scalandr added 3 commits January 19, 2023 20:17

fixes in ocean.cmake, mpas_ocn_nonhydrostatic_pressure_solve.F and mp…

18a8b6e

…as_ocn_vel_pressure_grad.F

corrected ocn_vertVel_vmix_tend_implicit in mpas_ocn_vmix.F

e209e96

cleaned code

1934c38

scalandr added 8 commits March 21, 2023 20:27

added nonhydro vars to restart file and corrected vvel coriolis

3e9778b

changes to hmix for vertical momentum

bc9cabc

added flags in bld folder and corrected ocn_vel_coriolis_nonhydro_ten…

9b31e2b

…d call

added lines to build the e3sm way on Anvil

b94a4c0

added del4 in vertical momentum

f9a76a1

merged with most recent version of e3sm

2dec6bd

fixed building issue

707cc46

Merge remote-tracking branch 'E3SM-Project/E3SM/master' into ocean/no…

a882fd9

…nhydro

scalandr force-pushed the ocean/nonhydro branch from 57c76f3 to a882fd9 Compare May 19, 2023 19:38

scalandr added 3 commits May 24, 2023 16:48

possible fix for e3sm tests

3759f61

possible fix for problem with ocn_vvel_hmix_del4

85718b4

fixed omp pragma in mpas_ocn_vmix.F

9052bc6

		($(FC) pio2.f90 $(FCINCLUDES) $(FFLAGS) $(LDFLAGS) $(LIBS) -o pio2.out &> /dev/null && echo "=> PIO 2 detected") \|\| \
		($(FC) pio2.f90 $(FCINCLUDES) $(FFLAGS) $(LDFLAGS) $(LIBS) -o pio2.out && echo "=> PIO 2 detected") \|\| \

		<nml_option name="config_use_vertMom_del2" type="logical" default_value=".false." units="unitless"
		description="If true, Laplacian horizontal mixing is used on the vertical momentum equation."

		call mpas_timer_stop("RK4 vel/thick tendency computations")

		! Update halos for prognostic variables.

added nonhydrostatic capabilities to MPAS-O #21

Are you sure you want to change the base?

added nonhydrostatic capabilities to MPAS-O #21

Conversation

scalandr commented May 20, 2022

xylar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scalandr commented Jul 1, 2022

xylar commented Jul 1, 2022

xylar commented Jul 1, 2022

scalandr commented Jul 1, 2022 • edited Loading

xylar commented Jul 1, 2022

xylar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

xylar commented Jul 5, 2022

scalandr commented Jul 13, 2022

scalandr commented Jul 13, 2022 • edited Loading

xylar commented Jul 13, 2022

xylar commented Jul 13, 2022

vanroekel commented Jul 13, 2022

scalandr commented Jul 14, 2022

vanroekel left a comment

Choose a reason for hiding this comment

xylar commented Mar 10, 2023

scalandr commented Mar 13, 2023

xylar commented Mar 13, 2023

scalandr commented Jul 1, 2022 •

edited

Loading

scalandr commented Jul 13, 2022 •

edited

Loading