Skip to content

GOTM in MPAS-Ocean #704

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

Merged
merged 18 commits into from
Apr 23, 2021
Merged

Conversation

qingli411
Copy link

This PR adds the General Ocean Turbulence Model (GOTM) in MPAS-Ocean as an additional option for the vertical turbulence closure, especially for coastal simulations in which the present implementation of KPP in CVMix does not work well. GOTM is a turbulence closure library that includes a set of two equation models such as k-epsilon and the generic length scale model (GLS; Umlauf and Burchard, 2003).

@qingli411
Copy link
Author

qingli411 commented Sep 24, 2020

It has been tested in a single column test case following Section 5.1 of Kärnä, 2020, which is an updated version of the surface gradient forcing test case in #615. The test case is in compass test number 138.
image

The simulated velocity and viscosity profiles are compared with the analytical solution of Equations (67) and (69), and a simulation using constant viscosity of 0.01 m^2/s in CVMix.

  • The velocity profile
    image

  • The viscosity profile
    image

The results depends on the bottom drag coefficient, here set to 0.0173, corresponding to the value inferred from the analytical solution with a layer thickness of 0.06 m.

@qingli411
Copy link
Author

Currently the GOTM module reads in its own namelist file ('gotmturb.nml', included in the compass test) at the beginning of the simulation. The path of the GOTM namelist file can be set in the MPAS-Ocean namelist. This is coded in the GOTM library so we will need to customize GOTM or have a way to generate the GOTM namelist if we want to control the behavior of GOTM from the MPAS-Ocean namelist.

GOTM requires input of bottom drag coefficient (to compute the bottom friction velocity) and the bottom roughness length, both set to constant in the namelist right now. But later may use the values from the bottom drag parameters in other part of the MPAS-Ocean code once they are ready.

@qingli411
Copy link
Author

@sbrus89 Could you please take a look at this when you have time? Looks like I don't have the access to assign a reviewer...

@qingli411
Copy link
Author

The implementation of GOTM in MPAS-Ocean has also been tested in a JRA55do forced G-case before rebasing to the ocean/develop branch. I only ran it for two years to make sure it does not break anything in a global case. It gives shallower mixed layer depth than CVMix in most of the regions and is quite sensitive to the time step.

Here are some MPAS-Analysis results:

@mark-petersen
Copy link
Contributor

@qingli411 please post a gotmturb.nml file so we can test this PR. Also, tell us what flags need to be changed to exercise all parts of the code you wrote. Is it simply config_use_gotm = .true., or do we need to test with various other settings as well?

@qingli411
Copy link
Author

@mark-petersen The gotmturb.nml file is in the default_gotm test case (testing_and_setup/compass/ocean/periodic_planar/2.5km/default_gotm/gotmturb.nml). The default_gotm test case is the single column test case with GOTM I mentioned earlier. To activate GOTM, use config_use_gotm = .true.. Since GOTM replaces CVMix here, we also need to turn off CVMIx by config_use_cvmix = .false..

@mark-petersen
Copy link
Contributor

@qingli411 I rebased on the head of ocean/develop and then tested with gnu 5 debug and intel 19 debug on the nightly regression suite, with cvmix off and gotm on, and gotmturb.nml in every run directory. All the tests run, which is a good sign, as debug catches any out of bounds array index.

This is the output from the nightly regression suite:

 ** Running case Baroclinic Channel 10km - Default Test
      PASS
 ** Running case Global Ocean 240km - Init Test
      PASS
 ** Running case Global Ocean 240km - Performance Test
      PASS
 ** Running case Global Ocean 240km - Restart Test
   ** FAIL (See case_outputs/Global_Ocean_240km_-_Restart_Test for more information)
 ** Running case Global Ocean 240km - Partition Test
      PASS
 ** Running case Global Ocean 240km - RK4 Partition Test
      PASS
 ** Running case Global Ocean 240km - Thread Test
      PASS
 ** Running case Global Ocean 240km - RK4 Thread Test
      PASS
 ** Running case Global Ocean 240km - Analysis Test
      PASS
 ** Running case Global Ocean 240km - BGC Ecosys Test
      PASS
 ** Running case ZISO 20km - Smoke Test
      PASS
 ** Running case ZISO 20km - Smoke Test with frazil
      PASS
 ** Running case Baroclinic Channel 10km - Thread Test
      PASS
 ** Running case Baroclinic Channel 10km - Decomp Test
      PASS
 ** Running case Baroclinic Channel 10km - Restart Test
   ** FAIL (See case_outputs/Baroclinic_Channel_10km_-_Restart_Test for more information)
 ** Running case sub-ice-shelf 2D - restart test
   ** FAIL (See case_outputs/sub-ice-shelf_2D_-_restart_test for more information)

for both gnu and intel, the restart tests all fail. This means one of these:

  1. upon restart, the proper fields are not computed on init, or are in the wrong order.
  2. there is a field in local memory that is needed from step to step, but not known at the beginning of a restart.

either way, the 'state' of the system is not identical just after beginning a restart, as it would be if the model proceeded without restarting. Try to fix (1) first. The solution to (2) is to add variables to the restart file. But we only do that if they are absolutely needed. If they can be computed on init instead, that is better.

You can copy my restart test here, for your convenience:
/lustre/scratch3/turquoise/mpeterse/runs/nightly/pr704_qingli_gotm/ocean/global_ocean/QU240/restart_test
or the whole regression suite:
/lustre/scratch3/turquoise/mpeterse/runs/nightly/pr704_qingli_gotm

you can tar up that whole directory, copy it to your own space, and run nightly_ocean_test_suite.py. You will need to add the link to your own executable in the forward subdirectories.

@qingli411
Copy link
Author

@mark-petersen OK. I think it's probably due to the second reason. GOTM solves prognostic equations for TKE and a length scale and requires the state of some variables from the previous time step. I added a list of variables that GOTM needs in the statePool for that reason (See this line in src/core_ocean/Registry.xml). These arrays are used in src/core_ocean/shared/mpas_ocn_vmix_gotm.F to assign values from the previous time step to GOTM. So maybe we should save these state variables in the restart file? If so is there a way to save these variables to the restart file only when GOTM is used?

Related to this I have a question on the way to initialize these arrays in MPAS-O. Here I initialize the arrays only once using a saved variable first. Now I realized that I should probably do this in the init step and check config_do_restart so that these arrays do not get initialized (to zero) again when doing restart. What do you think?

<!-- FIELDS FOR GOTM -->
<var name="gotmVertViscTopOfCell" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2 s^{-1}"
description="GOTM: vertical viscosity defined at the cell center (horizontally) and top (vertically)"
packages="forwardMode;analysisMode"
Copy link
Contributor

Choose a reason for hiding this comment

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

Here you are adding 3D fields that are allocated whether gotm is in use or not. To fix that, add a new package in Registry.xml called gotmPKG and put all of your new variables in that package (not forwardMode;analysisMode). You will need to add code to turn on the package if the first gotm flag is true. Follow the example here:

src/core_ocean/driver/mpas_ocn_core_interface.F
238       !
239       ! test for variable bottom drag of momentum, variableBottomDragPKG
240       call mpas_pool_get_package(packagePool, 'variableBottomDragPKGActive', variableBottomDragPKGActive)
241       call mpas_pool_get_config(configPool, 'config_use_variable_drag', config_use_variable_drag)
242       if ( config_use_variable_drag ) then
243         variableBottomDragPKGActive = .true.
244       end if

For bfb restarts, add the prognostic variables that are needed for restart to the restart stream in the Registry.xml file here:

1550         <stream name="restart"

Also, think carefully about which variables go into the state var_struct. Those should really only be the variables we have a prognostic equation for (i.e. use a time-stepping method to advance). All the other variables that are secondary should go into the diagnostic var_struct.

The difference is important because the state variables have 2 time levels for the timestepping, and the diagnostics only have one timelevel. From Registry:

1883     <var_struct name="state" time_levs="2">
...
2220     <var_struct name="diagnostics" time_levs="1">

so think carefully about what what variables belong in state. Those get double the memory, and that is needed for the time-stepping routines. The diagnostics are computed from the state variables, but only need a single instance in memory, for one time slice.

Copy link
Author

Choose a reason for hiding this comment

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

This is done in 2f63289 and 6693f2b. I think all the GOTM variables here need to be state variables rather than diagnostics variables as values from the previous time step are required to step GOTM forward.

@mark-petersen
Copy link
Contributor

Related to this I have a question on the way to initialize these arrays in MPAS-O. Here I initialize the arrays only once using a saved variable first. Now I realized that I should probably do this in the init step and check config_do_restart so that these arrays do not get initialized (to zero) again when doing restart. What do you think?

Correct, you do not want to use a first logical that you flip. All modules have at init routine at the very bottom, and all the init routines are called on start-up. That is where you conduct any process that you do on start-up or restart, but not on every time step. This includes initialization of variables.

You can use the flag config_do_restart which is false at the very beginning of a simulation, and true for restart runs. So for prognostic variables you may initialize if config_do_restart=.false. and then depend on the restart read of a prognostic variable after that.

! initialize GOTM arrays at the first step
! MPAS-O (1:nlev+1) <-> GOTM (nlev:0)
if (first) then
!$omp do schedule(runtime)
Copy link
Contributor

@mark-petersen mark-petersen Sep 29, 2020

Choose a reason for hiding this comment

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

@qingli411 all our omp calls have changed format in #513. Please search for all !$omp lines in your PR and change to this:

      !$omp parallel
      !$omp do schedule(runtime) &
      !$omp private(k, ki)
      do iCell = 1, nCells
...

for all cell, edge, and vertex loops. The private variables are any index or scalar value assigned within the loop.

Copy link
Contributor

Choose a reason for hiding this comment

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

Note that the first index (here iCell) is not in the private line on purpose. That is because openMP knows to run thread-parallel on that first loop.

Copy link
Author

Choose a reason for hiding this comment

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

Initialization of the GOTM state variables are now in the init subroutine (see 7b00c93).

I also changed the format for the OpenMP directives. But since my knowledge of OpenMP is very limited I'm not sure whether it work with OpenMP... Here GOTM has its own initialization step (called in the subroutine ocn_vmix_gotm_init), which allocates memory for a list of 1D arrays (vertical profiles of TKE, length scales, viscosity etc). Then I loop over all the MPAS-O cells to call GOTM to update these variables for each cell (the GOTM state variables allocated in MPAS-O).

gotmEpsbTopOfCellNew(k, iCell) = gotm_epsb(ki)
end do
end do
!$omp end do
Copy link
Contributor

Choose a reason for hiding this comment

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

At the end of a cell, edge, or vertex loop, change all omp lines to:

      !$omp end do
      !$omp end parallel

Copy link
Author

Choose a reason for hiding this comment

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

Done

@qingli411
Copy link
Author

qingli411 commented Sep 30, 2020

@mark-petersen I have made the requested changes and run nightly_ocean_test_suite.py with gnu debug on badger. Looks like the tests all pass. For some reason the script stopped at the end when reporting the runtime (see the output below)...

Here is my nightly test suite directory
/lustre/scratch4/turquoise/qingli/pr704_qingli_gotm

Output from the nightly regression suite

 ** Running case Baroclinic Channel 10km - Default Test
      PASS
 ** Running case Global Ocean 240km - Init Test
      PASS
 ** Running case Global Ocean 240km - Performance Test
      PASS
 ** Running case Global Ocean 240km - Restart Test
      PASS
 ** Running case Global Ocean 240km - Partition Test
      PASS
 ** Running case Global Ocean 240km - RK4 Partition Test
      PASS
 ** Running case Global Ocean 240km - Thread Test
      PASS
 ** Running case Global Ocean 240km - RK4 Thread Test
      PASS
 ** Running case Global Ocean 240km - Analysis Test
      PASS
 ** Running case Global Ocean 240km - BGC Ecosys Test
      PASS
 ** Running case ZISO 20km - Smoke Test
      PASS
 ** Running case ZISO 20km - Smoke Test with frazil
      PASS
 ** Running case Baroclinic Channel 10km - Thread Test
      PASS
 ** Running case Baroclinic Channel 10km - Decomp Test
      PASS
 ** Running case Baroclinic Channel 10km - Restart Test
      PASS
 ** Running case sub-ice-shelf 2D - restart test
      PASS
 ** Running case Periodic Planar 20km - LIGHT particle region reset test
      PASS
 ** Running case Periodic Planar 20km - LIGHT particle time reset test
      PASS
TEST RUNTIMES:
Traceback (most recent call last):
  File "./nightly_ocean_test_suite.py", line 240, in <module>
    runtime = np.ceil(float(subprocess.check_output(
ValueError: could not convert string to float: 'file'

@mark-petersen
Copy link
Contributor

@qingli411 please compile with openmp on, i.e.

make $COMPILER CORE=ocean USE_PIO2=true OPENMP=true DEBUG=true GEN_F90=true

and you will see some mismatching !$omp statements.

@qingli411
Copy link
Author

@mark-petersen That should be fixed now.

@mark-petersen
Copy link
Contributor

@qingli411 try to run the nightly regression suite with openmp on and debug. I get a floating invalid error:

206 ocean_model        000000000168B22C  ocn_vmix_gotm_mp_         313  mpas_ocn_vmix_gotm.f90

note the line number is for *.f90, after the C-pre processor (use GEN_F90=true on compile to see it).

All this is printing to the screen.

  init_turbulence
    init_turbulence
        reading turbulence namelists..
        done.
        allocation memory..
        allocation memory..
        allocation memory..
        done.
        reading turbulence namelists..
        done.
        allocation memory..
        allocation memory..
        allocation memory..
        done.

        --------------------------------------------------------
        You are using the generic two-equation model
        with the following properties:

            exponent of k in psi-equation,    m  =   1.50000000000000
            exponent of l in psi-equation,    n  =  -1.00000000000000
            exponent of cm0 in psi-equation,  p  =   3.00000000000000
            cpsi1                                =   1.44000000000000
            cpsi2                                =   1.92000000000000
            cpsi3minus                           = -0.630000000000000
            cpsi3plus                            =   1.00000000000000
            sig_k                                =   1.00000000000000
            sig_psi                              =   1.30000000000000

            Value of the stability function
            in the log-law,                   cm0 =  0.526464710835324
            in shear-free turbulence,        cmsf =  0.731006013490010

            von Karman constant,           kappa =  0.415873799673737
            homogeneous decay rate,            d =  -1.08695652173913
            spatial decay rate (no shear), alpha =  -4.97449646698836
            length-scale slope (no shear),     L =  7.333286748097885E-002
            steady-state Richardson-number, Ri_st=  0.248302382206140
        --------------------------------------------------------

    init_tridiagonal

You should remove all printing to the screen (print statements). If you want to keep verbose debug printing, use call mpas_log_write( so it goes into the output files. In addition, you should have a local CPP flag for verbose output, and have it default off for the PR. For example, see

src/framework/
mpas_forcing.F:#define FORCING_DEBUG_WRITE(M) !call mpas_log_write( M )

where you can remove the ! when you want verbose output.

@mark-petersen
Copy link
Contributor

You may need to set

export OMP_NUM_THREADS=1

in your window before running to test the omp lines properly.

@qingli411
Copy link
Author

@mark-petersen OK I will do that. The printing to the screen is in the GOTM library... We will need to modify the GOTM repository to write that information to the log file.

@qingli411
Copy link
Author

@mark-petersen With the above fix (7c8fc6e) the nightly regression suite with openmp on and debug fails on the Thread and Decomp tests with

SIGFPE: Floating-point exception - erroneous arithmetic operation.
See the output of the nightly regression suite here.
 ** Running case Baroclinic Channel 10km - Default Test
      PASS
 ** Running case Global Ocean 240km - Init Test
      PASS
 ** Running case Global Ocean 240km - Performance Test
      PASS
 ** Running case Global Ocean 240km - Restart Test
      PASS
 ** Running case Global Ocean 240km - Partition Test
      PASS
 ** Running case Global Ocean 240km - RK4 Partition Test
      PASS
 ** Running case Global Ocean 240km - Thread Test
   ** FAIL (See case_outputs/Global_Ocean_240km_-_Thread_Test for more information)
 ** Running case Global Ocean 240km - RK4 Thread Test
      PASS
 ** Running case Global Ocean 240km - Analysis Test
      PASS
 ** Running case Global Ocean 240km - BGC Ecosys Test
      PASS
 ** Running case ZISO 20km - Smoke Test
      PASS
 ** Running case ZISO 20km - Smoke Test with frazil
      PASS
 ** Running case Baroclinic Channel 10km - Thread Test
   ** FAIL (See case_outputs/Baroclinic_Channel_10km_-_Thread_Test for more information)
 ** Running case Baroclinic Channel 10km - Decomp Test
   ** FAIL (See case_outputs/Baroclinic_Channel_10km_-_Decomp_Test for more information)
 ** Running case Baroclinic Channel 10km - Restart Test
      PASS
 ** Running case sub-ice-shelf 2D - restart test
      PASS
 ** Running case Periodic Planar 20km - LIGHT particle region reset test
      PASS
 ** Running case Periodic Planar 20km - LIGHT particle time reset test
      PASS
TEST RUNTIMES:
Traceback (most recent call last):
  File "./nightly_ocean_test_suite.py", line 240, in <module>
    runtime = np.ceil(float(subprocess.check_output(
ValueError: could not convert string to float: 'file'

I guess some variables were not correctly initialized in the Thread tests. I have set export OMP_NUM_THREADS=1 before running the test suite. Any suggestions?

@mark-petersen
Copy link
Contributor

@qingli411 thanks for your work on this. We don't want excessive output written, even to log files. Lines like allocation memory.. done. should really just be on when a verbose setting is on.

@qingli411
Copy link
Author

@mark-petersen I understand that. But I don't think we have an easy way to control how GOTM writes out these messages as it is an external library. Same is for the way it reads in its own namelist. We will need to have a customized version of GOTM in MPAS-O if we want to have better control.

@mark-petersen
Copy link
Contributor

Yes, I understand now. It is best to not customize GOTM, just use it straight. Are there actual Fortran print statements in GOTM? They don't have an internal verbose flag to shut them off? If every processor uses print, and we are using 1k to 10k processors, then GOTM would be completely unusable in E3SM. If so, please try to get GOTM to put the print statements within a CPP verbose flag. Even if it takes several months to go through, it's worth making an issue to them now. Also, better to work with them than split off your own version of the code. It's a really reasonable request.

@qingli411
Copy link
Author

@mark-petersen I'm still in contact with the GOTM developers on possible better solutions. But, with their help, here (a603b7e) is a workaround to suppress the GOTM messages which only involves changing one line of code in GOTM. So I added the change in the Makefile. What do you think?

@mark-petersen
Copy link
Contributor

rebased. Passes nightly regression suite with gnu and intel, both debug, gotm off. Also tested QU240 with gotm on with both.

mark-petersen added a commit that referenced this pull request Apr 22, 2021
jonbob added a commit to E3SM-Project/E3SM that referenced this pull request Apr 22, 2021
Update mpas-source: add GOTM vertical mixing library

Brings in a new mpas-source submodule with changes only to the ocean
core, as well as bld file updates corresponding to Registry changes. It
adds the General Ocean Turbulence Model (GOTM) in MPAS-Ocean as an
additional option for the vertical turbulence closure, especially for
coastal simulations in which the present implementation of KPP in CVMix
does not work well. GOTM is a turbulence closure library that includes a
set of two equation models such as k-epsilon and the generic length
scale model (GLS; Umlauf and Burchard, 2003). See
MPAS-Dev/MPAS-Model#704

GOTM is default off, but is compiled by E3SM.

[NML]
[BFB]
jonbob added a commit to E3SM-Project/E3SM that referenced this pull request Apr 23, 2021
Update mpas-source: add GOTM vertical mixing library

Brings in a new mpas-source submodule with changes only to the ocean
core, as well as bld file updates corresponding to Registry changes. It
adds the General Ocean Turbulence Model (GOTM) in MPAS-Ocean as an
additional option for the vertical turbulence closure, especially for
coastal simulations in which the present implementation of KPP in CVMix
does not work well. GOTM is a turbulence closure library that includes a
set of two equation models such as k-epsilon and the generic length
scale model (GLS; Umlauf and Burchard, 2003). See
MPAS-Dev/MPAS-Model#704

GOTM is default off, but is compiled by E3SM.

[NML]
[BFB]
@mark-petersen mark-petersen merged commit 8674574 into MPAS-Dev:ocean/develop Apr 23, 2021
qingli411 added a commit to qingli411/compass that referenced this pull request Apr 28, 2021
caozd999 pushed a commit to caozd999/compass that referenced this pull request Apr 29, 2021
mark-petersen added a commit to mark-petersen/MPAS-Model that referenced this pull request May 7, 2021
GOTM in MPAS-Ocean MPAS-Dev#704

This PR adds the General Ocean Turbulence Model (GOTM) in MPAS-Ocean as
an additional option for the vertical turbulence closure, especially for
coastal simulations in which the present implementation of KPP in CVMix
does not work well. GOTM is a turbulence closure library that includes a
set of two equation models such as k-epsilon and the generic length
scale model (GLS; Umlauf and Burchard, 2003).
mark-petersen added a commit to mark-petersen/MPAS-Model that referenced this pull request May 7, 2021
vanroekel pushed a commit to vanroekel/MPAS-Model that referenced this pull request Sep 20, 2021
GOTM in MPAS-Ocean MPAS-Dev#704

This PR adds the General Ocean Turbulence Model (GOTM) in MPAS-Ocean as
an additional option for the vertical turbulence closure, especially for
coastal simulations in which the present implementation of KPP in CVMix
does not work well. GOTM is a turbulence closure library that includes a
set of two equation models such as k-epsilon and the generic length
scale model (GLS; Umlauf and Burchard, 2003).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants