Bug fixes in Redi and new tapering options #445

vanroekel · 2020-02-06T21:30:56Z

This PR fixes a few issues with the Redi mixing implementation and also adds tapering based on Danabasoglu and McWilliams 95 and Large et al 1997.

vanroekel · 2020-02-06T21:32:18Z

Here is a plot of minimum temperature from the tests (black is the latest from this PR). This is an E3SM G-case. We are evaluating after 30 years of simulation.

vanroekel · 2020-02-06T21:33:23Z

src/core_ocean/shared/mpas_ocn_gm.F

+            dzTop(maxLevelCell(sVal)+1) = -1e-15_RKIND
+            dTdzTop(maxLevelCell(sVal)+1) = -1e-15_RKIND
+            dSdzTop(maxLevelCell(sVal)+1) = -1e-15_RKIND
+


I think we need to compute dtdz and dsdz each time. There are times when the stratification is computed for iCell but this is not consistent with iCellSelf. This makes most sense to me comparing to the picture in Griffies et al 1998

vanroekel · 2020-02-06T21:34:01Z

src/core_ocean/shared/mpas_ocn_tendency.F

+
+               modifiedGroupName = trim(groupItr % memberName) // "_redi_term3"
+               call mpas_pool_get_array(redi3Pool, trim(modifiedGroupName), redi_term3)
+


these output the individual redi terms (except vertical mixing). we can revert these when we want to merge.

vanroekel · 2020-02-06T21:55:21Z

src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F

@@ -216,12 +218,13 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, scratchPool, layerThicknessEdge,

               ! For top layer, only use triads pointing down:
               flux_term2 = coef * RediKappa(iCell) * layerThicknessEdge(k, iEdge) * &
-                  0.5_RKIND * &
+                  0.25_RKIND * &


I think we need to treat the top and bottom as an average of 4 things even though there are just two. I think of it as an average of 4 triads still but 2 are zero.

vanroekel · 2020-02-06T21:56:30Z

src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F

@@ -317,7 +323,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, scratchPool, layerThicknessEdge,
               flux_term3 = RediKappa(iCell) * 2.0_RKIND * &
                   (fluxRediZTop(iTr,k)/areaSum(k) - fluxRediZTop(iTr,k+1)/areaSum(k+1))

-               tend(iTr,k,iCell) = tend(iTr,k,iCell) + flux_term3
+               redi_term3(iTr,k,iCell) = redi_term3(iTr,k,iCell) - flux_term3
+               tend(iTr,k,iCell) = tend(iTr,k,iCell) - flux_term3


This seems to be the most critical change. It was not clear to me why flux_term3 had a positive on it and was the dominant term of 1-3 in my testing that seemed to cause the large negative temperatures. This seems to fit best with the equation for the Redi term you had on confluence.

Term 3 is d/dz ( h S dot grad phi) = ( S dot grad phi) where grad phi is computed on the edge and middle vertically, and then multiplied by the slope triad that goes through that edge. I don't expect a negative sign here. But just in case, I made this change in the current head of ocean/develop and reran this idealized SO test case.

Compare to identical repo with positive sign on term 3:
#280 (comment)

Also, the analytic polynomial test cases match with a positive sign but not a negative sign.

Thanks @mark-petersen, I agree the sign is correct as positive, I still don't fully understand why though.

vanroekel · 2020-02-12T15:57:41Z

@mark-petersen I've been thinking more about redi term3 and had a question about the fluxRediZtop code (here

MPAS-Model/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F

Lines 303 to 313 in 600fa0a

    
                          do iTr = 1, nTracers 
        
                             ! Gradient of tracer at edge. 
        
                             ! dvEdge on following line is equal to areaEdge/dcEdge 
        
                             dTracerDx = (tracers(iTr,k,cell2) - tracers(iTr,k,cell1)) * dvEdge(iEdge) 
        
                             fluxRediZTop(iTr,k) = fluxRediZTop(iTr,k) + & 
        
                                slopeTriadUp(k,iCellSelf,iEdge) * dTracerDx 
        
                             fluxRediZTop(iTr,k+1) = fluxRediZTop(iTr,k+1) + & 
        
                                slopeTriadDown(k,iCellSelf,iEdge) * dTracerDx 
        
                          end do 
        
                       end do 
        
                    end do

). When I look at this it seems there are only two triads being considered for each edge. The previous code (redi term2) has an explicit loop over iCellSelf (including 1 and 2) but here only seems to have one value per edge. Shouldn't this be an average of the 4? I'm sure I'm not understanding something in the implementation.

vanroekel · 2020-02-19T05:12:25Z

@mark-petersen I think this PR is ready now. Here is the latest plot of temperature min from @jonbob's E3SM run. There is something weird at year 20, but overall this looks good.

mark-petersen · 2020-02-19T05:18:52Z

That plot is with the single-cell pits, right? When the single-cell pits were removed, did the min remain at -1.8, even if it was just MPAS-Ocean stand alone?

milenaveneziani · 2020-02-19T05:47:52Z

Aren't you concerned about all the values with T<-2?

vanroekel · 2020-02-19T17:39:08Z

The min temperature remained above -2 for most of the run. Something weird is indeed happening at year 20. It does recover to be above -2 again quickly. I will take a look at that time period. Other than those few years things look good.

This is the case with no single cell pits. So I'm not sure why we are getting that period. Jon saved extra output so I should be able to see what is happening.

jonbob · 2020-02-19T17:57:56Z

@vanroekel - let me know if I should rerun that time. And I'll check the output to make sure nothing suspicious has happened....

jonbob · 2020-02-19T23:09:39Z

I looked more closely at the output and the log files -- there's an almost linear decrease in the minimum temperature starting very close to the beginning of year 20. Nothing in the logs jumps out to indicate something going off the rails

vanroekel · 2020-02-20T06:31:19Z

@mark-petersen I've looked through the 30 years of simulation from @jonbob. The negatives are most often right near the bottom, but there are a few further up in the column (layer 14 when max level cell = 21). So it looks like bottom tapering could help, but may not completely fix the problem.

mark-petersen · 2020-02-20T17:43:12Z

OK. It sounds like bottom tapering is still worth trying.

mark-petersen · 2020-02-24T20:52:26Z

src/core_ocean/shared/mpas_ocn_gm.F

+               ! Griffies 1998 eqn 34
+               k33(k  ,iCell) = k33(k  ,iCell) + slopeTaperUp * sfcTaperUp * &
+                              areaEdge*dzTop(k)*slopeTriadUp(k,iCellSelf,iEdge)**2
+               k33(k+1,iCell) = k33(k+1,iCell) + slopeTaperDown * sfcTaperDown * &


@vanroekel, these k33 lines were moved up in the last commit, but they still have the taper coefficients, same as slopeTriadUp and slopeTraidDown. The plan was to not taper k33, right? So shouldn't the slope taper variables not appear in the above lines?

vanroekel · 2020-02-25T21:15:09Z

update, after the new commit, here are the results of min temperature after 1 year. Blue is before the latest commit and orange is the latest

milenaveneziani · 2020-02-25T21:24:41Z

Looks great!

vanroekel · 2020-02-28T18:53:56Z

Okay, I've gone through and made a number of new changes, including tapering based on N2 (like Danabasoglu and Marshall 1997), I've also added a convective trigger. I've tested briefly in stand alone, things are looking good so far!

vanroekel · 2020-02-28T18:54:32Z

I should mention the orange line is the commit with the r value fix, and blue is the latest commit.

milenaveneziani · 2020-02-28T18:56:35Z

Great. I can run a G-case whit E3SM when we are ready.

vanroekel · 2020-02-28T21:50:24Z

Thanks @milenaveneziani for offering to do the run. I had asked @jonbob to pull the code together for the E3SM run and he already submitted the test. 4 months in, things look good!

Also, through 6 months of stand alone results look good as well.

vanroekel · 2020-03-02T22:56:43Z

@mark-petersen this branch now has options to disable all tapering. To remove tapering on everything do

config_Redi_diabatic_limit = 1e6
config_Redi_use_slope_taper = .false.
config_Redi_use_surface_taper = .false.
config_Redi_N2_limit_term1 = .false.
config_GM_closure = 'constant'

Jon is testing the new limiting of term1 in a g-case now.

vanroekel · 2020-03-23T17:37:18Z

@mark-petersen I just pushed my latest code with some changes to tapering and disabling terms 2 and 3 in the bottom cell. Could you test your convergence on this commit and let me know how it goes.

mark-petersen · 2020-03-25T20:59:34Z

@vanroekel I find these descriptions confusing:

 398     <nml_option name="config_Redi_half_slope" type="real" default_value="0.004" units="non-dimensional"
 399         description="value of slope where set to half max"
...
 402     <nml_option name="config_Redi_half_width" type="real" default_value="0.001" units="non-dimensional"
 403         description="value of half width of tanh function limiter"

Could you replace the description with something like

nml_option name="config_Redi_half_slope"
"The isopycnal slope used in the Redi calculation is limited by the 
function $equation here$ (see Authors et al eqn. XX). 
config_Redi_half_slope is $\variable$ and is the slope value where the limiting is 1/2 
(i.e. the center of the tanh function)"

and similar for config_Redi_half_width. Thanks!

mark-petersen · 2020-04-15T21:31:57Z

@vanroekel The head is still bfb with your last commit 2f74e37 on the MPAS-Ocean nightly regression suite. Please check in your standard E3SM set-up. Should be bfb there as well, as changes were all omp and formatting, and removing unused variables.

This is almost ready to merge. I'm rechecking convergence and idealized Redi test now.

vanroekel · 2020-04-15T21:47:38Z

@mark-petersen I'm building a case with your latest changes now. I should have something for you by tomorrow morning at the latest.

vanroekel · 2020-04-16T03:22:22Z

@mark-petersen I can confirm that the latest commit is BFB with my previous e3sm tests.

mark-petersen · 2020-04-16T03:46:18Z

Automated tests all look good.

mark-petersen

This now passes all the tests and is read to merge.

vanroekel · 2020-04-16T03:47:46Z

Great! Thanks for all the help getting this across the finish line. I agree this is ready to merge.

Matches bfb with previous

mark-petersen · 2020-04-16T13:58:16Z

src/core_ocean/Registry.xml

@@ -1647,7 +1626,8 @@
 				mode="forward;analysis">

 			<stream name="mesh"/>
-			<var_struct name="tracers"/>
+      <var_struct name="tracers"/>
+      <var name="indMLD"/>


@vanroekel I'm taking a last careful look through. We added indMLD to the restart file here when we were struggling with the restart bug, but I don't think it is needed. The restart test passes without it, with Redi and GM on. Do you agree it can be removed?

I think so. I believe rediGM is now called before the timestep but after the analysis compute, but this might require that MLD be computed on start up. Is that flag set in your testing? If we don't compute on start up we may want to keep indMLD in the restart file.

I guess it depends on the timing of the mixedLayerDepths AM. We can only safely remove the variable if the AM is computed just before the restart is written, to get exact restarts. For example, if we use either of these settings:

config_AM_mixedLayerDepths_compute_interval = 'dt' config_AM_mixedLayerDepths_compute_interval = '01_00:00:00'

for a global run, we can remove indMLD from the restart file. The safest option is to leave it in, at the cost of a slightly larger restart file.

I agree for nearly all E3SM cases this won't be an issue, but it might be safest to just leave it in.

mark-petersen · 2020-04-16T14:07:31Z

src/core_ocean/Registry.xml

@@ -1657,6 +1637,7 @@
 			<var name="xtime"/>
 			<var name="simulationStartTime"/>
 			<var name="boundaryLayerDepth"/>
+			<var name="indexBoundaryLayerDepth"/>


Same thing here. My restart test passes using QU240 and cvmix on, using kpp, with both boundaryLayerDepth and indexBoundaryLayerDepth removed. Are there conditions where those variables are used from the previous timestep, rather than recomputed?

agreed. I think this was from your first tapering implementation when it was based on BLD. I'm sure it can be removed.

What about boundaryLayerDepth in the restart file? That was added a long time ago. Probably before the cvmix:

cb8195861b (Douglas Jacobsen 2015-02-02 13:45:52 -0700 1659) <var name="boundaryLayerDepth"/>

I can't think of any reason to include boundaryLayerDepth in the restart file. It is always computed before mixing and is only used within the mixing routine. So I don't think that should affect restarts.

removed: boundaryLayerDepth indexBoundaryLayerDepth

xylar · 2020-04-16T18:54:21Z

Congratulations! This was clearly a lot of work!

mark-petersen · 2020-06-26T17:40:47Z

src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F

-           call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
-        end if
+ !       if (config_use_GM.or.config_use_Redi) then
+ !          call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &


@vanroekel I saw that this line is commented out in ocean/develop, and looks like we were in a rush and left some confusing code. This was commented out in this PR in:
66a13d9 fixes for restart test
We changed E3SM to call this routine from the driver:

components/mpas-ocean/driver/ocn_comp_mct.F:968: call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, meshPool, scratchPool, timeLevelIn=1)

in E3SM-Project/E3SM@b29cfe5b370 (Luke Van Roekel 2020-05-05 16:01:43 -0500 968)
but it appears that it is not called from the driver in stand-alone. I can make an issue to clean this up, but wanted to ping you on these lines.

yes @mark-petersen I removed those calls and moved them to the driver, I just never removed them. I can clean this up if you want. It should be in stand alone as you note below

mark-petersen · 2020-06-26T17:41:48Z

src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F

-         end if
+!         ! Compute normalGMBolusVelocity; it will be added to normalVelocity in Stage 2 of the next cycle.
+!         if (config_use_GM.or.config_use_Redi) then
+!            call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &


This line was also commented out, and needs to be removed or cleaned up.

yes those can be removed

mark-petersen · 2020-06-26T17:43:08Z

src/core_ocean/mode_forward/mpas_ocn_forward_mode.F

           call ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, diagnosticsPool, statePool, err)

+           ! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on
+           if (config_use_Redi.or.config_use_GM) then
+               call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, meshPool, scratchPool, timeLevelIn=1)


Looks like the stand-alone call is here, so maybe this is just a matter of removing the commented out lines for clarity.

mark-petersen · 2020-06-26T17:45:08Z

src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F

-         if (config_use_standardGM) then
-             call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
+         if (config_use_GM) then
+             call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &


In rk4, the corresponding lines here should also be removed, since it is called from the driver.

…/develop Bug fixes in Redi and new tapering options MPAS-Dev#445 This PR fixes a few issues with the Redi mixing implementation and also adds tapering based on Danabasoglu and McWilliams 95 and Large et al 1997.

The mixed layer depth analysis member was enabled and writing out analysis at every time step for all global ocean forward model runs. This appears to have been turned on for debugging in a pull request from more than a year ago (MPAS-Dev/MPAS-Model#445). During dynamic adjustment, this can lead to files that are many GB in size and presumably slows down all model runs significantly.

The mixed layer depth analysis member was writing out analysis from every time step for all global ocean forward model runs. The MLD AM was turned on for every time step in MPAS-Dev/MPAS-Model#445 and output in legacy compass was set to every day. But in the compass package, output is every time step. During dynamic adjustment, this can lead to files that are many GB in size and presumably slows down all model runs significantly.

vanroekel requested a review from mark-petersen February 6, 2020 21:30

vanroekel commented Feb 6, 2020

View reviewed changes

xylar added bug Ocean labels Feb 9, 2020

xylar assigned mark-petersen Feb 9, 2020

mark-petersen mentioned this pull request Feb 13, 2020

Fix Redi mixing in ocean E3SM-Project/E3SM#3391

Merged

mark-petersen reviewed Feb 24, 2020

View reviewed changes

mark-petersen force-pushed the ocean/rediFixesAndNewTapering branch from 7b96a58 to b70416f Compare March 3, 2020 19:59

mark-petersen force-pushed the ocean/rediFixesAndNewTapering branch from f751f53 to ef2be48 Compare April 15, 2020 21:27

Update Redi convergence and idealized cases

2342524

mark-petersen self-requested a review April 16, 2020 03:46

mark-petersen approved these changes Apr 16, 2020

View reviewed changes

Comments, performance change.

3a113f2

Matches bfb with previous

mark-petersen reviewed Apr 16, 2020

View reviewed changes

Remove BLD variables from restart stream

ba03b7c

removed: boundaryLayerDepth indexBoundaryLayerDepth

mark-petersen merged commit 8326fbb into MPAS-Dev:ocean/develop Apr 16, 2020

This was referenced Apr 17, 2020

Missing stream(s) in spin-up configs #521

Closed

Errors in CUSP 12 spin-up #525

Closed

mark-petersen reviewed Jun 26, 2020

View reviewed changes

vanroekel mentioned this pull request Sep 8, 2020

Revert GM buoyancy gradient calculation to v1 form #687

Merged

matthewhoffman added bug Ocean labels Mar 17, 2021

xylar mentioned this pull request Aug 13, 2021

Mixed layer depth analysis member in global ocean test cases outputs every time step MPAS-Dev/compass#211

Closed

xylar mentioned this pull request Aug 13, 2021

Remove mixed layer depth AM output from global ocean steps MPAS-Dev/compass#212

Merged


		modifiedGroupName = trim(groupItr % memberName) // "_redi_term3"
		call mpas_pool_get_array(redi3Pool, trim(modifiedGroupName), redi_term3)

Bug fixes in Redi and new tapering options #445

Bug fixes in Redi and new tapering options #445

Conversation

vanroekel commented Feb 6, 2020

vanroekel commented Feb 6, 2020

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

vanroekel commented Feb 12, 2020 • edited Loading

vanroekel commented Feb 19, 2020

mark-petersen commented Feb 19, 2020

milenaveneziani commented Feb 19, 2020

vanroekel commented Feb 19, 2020

jonbob commented Feb 19, 2020

jonbob commented Feb 19, 2020

vanroekel commented Feb 20, 2020

mark-petersen commented Feb 20, 2020

Choose a reason for hiding this comment

vanroekel commented Feb 25, 2020

milenaveneziani commented Feb 25, 2020

vanroekel commented Feb 28, 2020

vanroekel commented Feb 28, 2020

milenaveneziani commented Feb 28, 2020

vanroekel commented Feb 28, 2020

vanroekel commented Mar 2, 2020

vanroekel commented Mar 23, 2020

mark-petersen commented Mar 25, 2020

mark-petersen commented Apr 15, 2020

vanroekel commented Apr 15, 2020

vanroekel commented Apr 16, 2020

mark-petersen commented Apr 16, 2020

mark-petersen left a comment

Choose a reason for hiding this comment

vanroekel commented Apr 16, 2020

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 Apr 16, 2020

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

vanroekel commented Feb 12, 2020 •

edited

Loading