Remove single-cell pits when initializing global bathymetry #450

mark-petersen · 2020-02-13T23:30:12Z

Before this, single-cell pits were allowed to be one layer deeper than all surrounding cells. These cells can collect tracer extrema and are not able to use horizontal diffusion or advection to clear them out. This PR reduces pits to make them level with the next deepest neighbor cell.

mark-petersen · 2020-02-13T23:31:18Z

Testing now...

xylar · 2020-02-14T05:33:55Z

I was wondering what the logic was in having these values previously. Could you rename the config option following my work in #436 as well?
c4993f3
This PR seems like a more logical place to rename this config option. And I must insist that filling holes isn't smoothing in the sense that a user would expect.

vanroekel

this code needs iterations. I looked at a 60to30 mesh using this initialization and ~700 cells that need alteration of max level cell.

src/core_ocean/mode_init/mpas_ocn_init_global_ocean.F

mark-petersen · 2020-03-09T14:52:07Z

@xylar this change is the same place as #440. It makes more sense to close this PR and make this minor change for single-cell pits in #440, then review it all together. Do you agree? If so, I'll close this PR and add a commit to #440.

xylar · 2020-03-09T15:05:35Z

No I think these need to be separate PRs. I don't want any additional confusion about what we mean by "smoothing" and if we start modifying filling holes in a PR about smoothing, that would do exactly that in my view.

mark-petersen · 2020-03-17T20:53:30Z

rebased. testing now.

mark-petersen · 2020-03-29T15:45:59Z

src/core_ocean/mode_init/mpas_ocn_init_global_ocean.F

+          ! Alter cells for partial bottom cells before filling in bathymetry
+          ! holes.
+          do iCell = 1, nCells
+             call ocn_alter_bottomDepth_for_pbcs(bottomDepth(iCell), refBottomDepth, maxLevelCell(iCell), iErr)


@xylar this is now complete. The critical change was the line above. The bathymetry holes would be filled, and then in a later call the cells were altered. Specifically, if the bottom bathymetry is within 5% of the layer thickness to the top of the cell, we round it up and change maxLevelCell to one less.

My solution here is to call the PBC routine before filling the bathymetry. I left the original call to
ocn_alter_bottomDepth_for_pbcs
because it is later in the sequence, and it is within your iterations for SSH for land ice. I thought there was no harm done leaving the later one in place.

I agree, I think this is safe to do. I believe I had ocn_alter_bottomDepth_for_pbcs in ocn_init_vertical_grid because it currently only gets called if you are not constraining the vertical grid base don the Haney number, which handles partial bottom cells after making substantial alternations to the grid. This earlier call might slightly modify the outcome of the later stage but not in a way that I expect to be important (just probably not bit-for-bit).

mark-petersen · 2020-03-29T15:52:05Z

I tested this with the EC60to30 in COMPASS, and this little script:
https://gist.github.com/mark-petersen/52af011b559786a2d7206e4c225ea992
which showed single-cell holes before, and now shows there are none.

The other changes in this PR are just a little clean-up, and moving init flags to the init_mode Registry.xml file.

mark-petersen · 2020-03-29T15:53:04Z

@xylar and @vanroekel if either of you want to make comments or test, please do.

xylar · 2020-03-29T16:52:08Z

src/core_ocean/mode_init/mpas_ocn_init_global_ocean.F

-             call mpas_allocate_scratch_field(smoothedLevelsField, .true.)
-


I think I agree that the temporary field isn't needed but want to talk it through. If a cell is deeper than its neighbors, it will get filled in. By definition, it is now at the same maxLevelCell as its deepest neighbor. Therefore, no neighbor can be deeper than this cell if it was altered and no neighbors will be filled. So it should be safe to update maxLevelCell as we go.

Yeah, I went through the same logic, and came to the same conclusion.

xylar · 2020-03-29T17:08:44Z

I'd like to do some testing tomorrow.

xylar

In addition to a halo update, I think a similar approach to hole-filling is needed for cases with ice-shelf cavities. I will push two more commits to this branch as soon as I have tested them.

xylar · 2020-03-30T09:13:29Z

src/core_ocean/mode_init/mpas_ocn_init_global_ocean.F

             maxLevelCell(nCells+1) = -1
-             smoothedLevelsField % array = maxLevelCell

             do iCell = 1, nCellsSolve


@mark-petersen, this computation is done over nCellsSolve but I can't find any indication that there is a halo update of bottomDepth, which seems like a potential problem.

xylar · 2020-03-30T13:36:34Z

@mark-petersen and @vanroekel, I tested including my 2 commits with QU240 and QU240wISC on my laptop (with gnu) and everything went fine.

I tested EC60to30 and EC60to30wISC on grizzly with intel, and was able to run init successfully. I had CFL problems in both setups during spin_up but I suspect this is unrelated to this PR. This PR is also not the right place to address those problems.

I used the following script to check that maxLevelCell had the desired property in init/initial_state.nc for each of the 4 setups listed above:

#!/usr/bin/env python

import xarray
import argparse
import numpy


parser = argparse.ArgumentParser()
parser.add_argument("-i", "--infile", dest='infile', required=True,
                    help="MPAS file with maxLevelCell to check")
args = parser.parse_args()

ds = xarray.open_dataset(args.infile)

cellsOnCell = ds.cellsOnCell.values - 1

maxLevelCell = ds.maxLevelCell.values - 1

maxLevelNeighbor = numpy.zeros(maxLevelCell.shape, int)
for index in range(ds.sizes['maxEdges']):
    coc = cellsOnCell[:, index]
    mask = coc >= 0
    maxLevelNeighbor[mask] = numpy.maximum(maxLevelNeighbor[mask],
                                           maxLevelCell[coc[mask]])

mask = maxLevelCell >= 0

assert(numpy.all(maxLevelCell[mask] <= maxLevelNeighbor[mask]))

xylar

I'm happy with this, now that pits are addressed in the case with ice-shelf cavities, too.

xylar · 2020-03-30T14:09:56Z

I suspect that the changes related to reverting Redi d3667cc might have also reverted some changes to the EC60to30 setups that made it possible to spin them up (and aren't related directly to Redi).

vanroekel

ran current head of commit for 60to30 initialization and it works, producing no single cell pits.

mark-petersen · 2020-03-31T02:38:34Z

@xylar I ran the these, with the following results:

EC60to30 on this PR: init and spin-up complete.
EC60to30wISC on this PR: init and SSH adjustment works. Spin-up dies, as you noted.
EC60to30wISC on ocean/develop: init works. SSH adjustment dies on first iteration (go figure).

I think the EC60to30wISC spin-up can work by adjusting Rayleigh damping and time step. I would rather make those adjustments when we have everything in place (64 levels, new mesh, etc), rather than get it to work for this one PR. If you agree, I'll merge it - perhaps after #495

xylar · 2020-03-31T05:01:11Z

@mark-petersen,

EC60to30 on this PR: init and spin-up complete.

That's strange, it didn't work for me. Did you change anything?

I think the EC60to30wISC spin-up can work by adjusting Rayleigh damping and time step. I would rather make those adjustments when we have everything in place (64 levels, new mesh, etc), rather than get it to work for this one PR. If you agree, I'll merge it - perhaps after #495

Our agreement was that we were not going to change to 64 vertical levels for EC60to30wISC for now. So we need to set this up as is, and get spin-up working. This is quite a high priority, since we have to start burning through our time on NERSC before June 30th.

xylar · 2020-03-31T06:07:04Z

I will keep working on spin-up today, so we can hopefully get this IC in place before #495

xylar · 2020-03-31T12:41:00Z

I think this is not working (anymore). maxLevelCell = 60 everywhere and bottomDepth = 5500 everywhere in my EC60to30wISC testing.

xylar · 2020-03-31T12:43:14Z

It's probably the changes I made to try to fill pits in the Haney-number constrained initialization. I'm going to try again with the commit before that one and see if things look normal.

mark-petersen · 2020-03-31T15:23:45Z

This is definitely something I introduced. It occurred without ISC as well. I'm debugging now, but have calls interspersed this morning.

mark-petersen · 2020-03-31T15:39:54Z

@xylar I'm sorry I wasted your time on this. Bug fix is in. I had deleted a print statement and deleted the exit on the next line by mistake. You can proceed with your ISC testing.

xylar · 2020-03-31T16:06:24Z

Okay, great! Thanks for finding that. Then we all need to re-test because all our "tests" passed trivially with a flat bottom... @vanroekel, can you run a quick test when you have time?

xylar · 2020-03-31T16:07:01Z

I'm going to start again with my EC60to30wISC spin-up starting with this branch and modifying the stages of spin-up if needed.

vanroekel · 2020-03-31T16:12:18Z

sure, I should be able to do that today. Would it be enough to run init again? Or would you like me to run the spin up as well?

mark-petersen · 2020-03-31T16:14:10Z

I'm running EC60to30 init and spinup now. I'll post some plots.

xylar · 2020-03-31T16:28:39Z

Would it be enough to run init again? Or would you like me to run the spin up as well?

Just running the test you did before should be enough, given the spin-ups @mark-petersen and I are doing.

vanroekel · 2020-03-31T16:57:17Z

okay, happy to do that. I'll let you know what I find.

xylar · 2020-03-31T17:18:10Z

@mark-petersen, something is not working right with Haney-number constrained initialization at the moment but it's a problem with ocean/develop as well, so not related to this PR. I'll look into it tomorrow but I don't have the energy for it today. Anyway, that is to say that I won't be able to test this PR on EC60to30wISC until I figure out why the Haney Number is comping out as 10572.8534636064!!!

vanroekel · 2020-03-31T17:47:33Z

I just retested and the result is good, no single cell pits.

xylar · 2020-04-01T09:40:45Z

@mark-petersen, I somehow missed linking in (and using) the land-ice mask from cell culling in the EC60to30wISC initial_state step. This is now fixed in #499. After cherry-picking #499 onto this branch, I was able to create an EC60to30wISC mesh that doesn't have any pits and looks sensible in paraview. I'm working on spinning it up now.

xylar · 2020-04-01T10:25:46Z

The EC60to30wISC spin-up was also successful, so we should be good to go (for the maint-1.2 run) as soon as this, #499 and #495 get merged. I don't want to create the initial condition from my current branch because I would like to be able to point to a specific commit on ocean/develop that it was generated from.

Remove single-cell pits from ocean bathymetry Add PBC call just before filling bathymetry holes Move init flags into mode_init

mark-petersen · 2020-04-02T04:51:27Z

rebased. retesting now.

mark-petersen · 2020-04-02T05:15:28Z

Passed nightly regression with debug, and also QU240wISC.

mark-petersen · 2020-04-02T05:34:19Z

@xylar ocean/develop is ready to go. I am running the EC60to30wISC right now here:

/lustre/scratch4/turquoise/mpeterse/runs/200401_EC60to30wISC/ocean/global_ocean/EC60to30wISC/init

remember to use the directory initial_state_60_layer. I renamed it to initial_state so all the links work, i.e.

mv initial_state initial_state_renamed_64
mv initial_state_60_layer initial_state

…evelop Before this, single-cell pits were allowed to be one layer deeper than all surrounding cells. These cells can collect tracer extrema and are not able to use horizontal diffusion or advection to clear them out. This PR reduces pits to make them level with the next deepest neighbor cell.

mark-petersen added Ocean COMPASS labels Feb 13, 2020

mark-petersen requested a review from vanroekel February 13, 2020 23:30

mark-petersen self-assigned this Feb 13, 2020

mark-petersen changed the base branch from master to ocean/develop February 13, 2020 23:30

vanroekel suggested changes Feb 16, 2020

View reviewed changes

src/core_ocean/mode_init/mpas_ocn_init_global_ocean.F Outdated Show resolved Hide resolved

mark-petersen force-pushed the ocean/remove_pits_in_global_init branch from 0769284 to 105fdef Compare March 17, 2020 20:53

mark-petersen added this to the Bathymetry, mesh, IC for E3SM V2 milestone Mar 19, 2020

mark-petersen force-pushed the ocean/remove_pits_in_global_init branch from 105fdef to 2707bcf Compare March 29, 2020 15:31

mark-petersen commented Mar 29, 2020

View reviewed changes

xylar reviewed Mar 29, 2020

View reviewed changes

xylar reviewed Mar 30, 2020

View reviewed changes

xylar approved these changes Mar 30, 2020

View reviewed changes

vanroekel approved these changes Mar 30, 2020

View reviewed changes

mark-petersen force-pushed the ocean/remove_pits_in_global_init branch from 0c0d5a7 to a5513d1 Compare March 31, 2020 15:37

mark-petersen mentioned this pull request Mar 31, 2020

Add script to compute vertical coordinate for initial condition #495

Merged

mark-petersen and others added 4 commits April 1, 2020 22:49

Fill bathymetry holes

78c8ac6

Remove single-cell pits from ocean bathymetry Add PBC call just before filling bathymetry holes Move init flags into mode_init

Bug fix - flat bottom bug.

ad880d0

Add halo update for bottomDepth

97733bf

Remove holes in Haney-number init

b0969b0

mark-petersen force-pushed the ocean/remove_pits_in_global_init branch from a5513d1 to b0969b0 Compare April 2, 2020 04:50

mark-petersen merged this pull request into MPAS-Dev:ocean/develop Apr 2, 2020

mark-petersen deleted the ocean/remove_pits_in_global_init branch April 2, 2020 05:16

matthewhoffman added COMPASS Ocean labels Mar 17, 2021

		call mpas_allocate_scratch_field(smoothedLevelsField, .true.)

Remove single-cell pits when initializing global bathymetry #450

Remove single-cell pits when initializing global bathymetry #450

Uh oh!

Conversation

mark-petersen commented Feb 13, 2020

Uh oh!

mark-petersen commented Feb 13, 2020

Uh oh!

xylar commented Feb 14, 2020

Uh oh!

vanroekel left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

mark-petersen commented Mar 9, 2020

Uh oh!

xylar commented Mar 9, 2020

Uh oh!

mark-petersen commented Mar 17, 2020

Uh oh!

mark-petersen Mar 29, 2020

Choose a reason for hiding this comment

Uh oh!

xylar Mar 29, 2020

Choose a reason for hiding this comment

Uh oh!

mark-petersen commented Mar 29, 2020

Uh oh!

mark-petersen commented Mar 29, 2020

Uh oh!

xylar Mar 29, 2020

Choose a reason for hiding this comment

Uh oh!

mark-petersen Mar 30, 2020

Choose a reason for hiding this comment

Uh oh!

xylar commented Mar 29, 2020

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

xylar Mar 30, 2020

Choose a reason for hiding this comment

Uh oh!

xylar commented Mar 30, 2020

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

xylar commented Mar 30, 2020

Uh oh!

vanroekel left a comment

Choose a reason for hiding this comment

Uh oh!

mark-petersen commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

mark-petersen commented Mar 31, 2020

Uh oh!

mark-petersen commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

vanroekel commented Mar 31, 2020

Uh oh!

mark-petersen commented Mar 31, 2020

Uh oh!

xylar commented Mar 31, 2020

Uh oh!

vanroekel commented Mar 31, 2020

Uh oh!