Neutral Diffusion: Add option for discontinuous reconstruction of T/S#553
Conversation
Some initial code changes to update the neutral_diffusion algorithm to find neutral surface positions using the discontinuous T and S reconstructions. Changes: - Break out a subroutine to return both the left and right edges of PPM reconstructions within a cell - Update logic to neutral_diffusion_calc_coeffs so that either the continuous or discontinuous reconstructions can be used ToDo: - Come up with a search algorithm that works on each subinterface - Update interpolate_for_nondim_position so that it returns a more accurate position of the neutral surface within each layer
search for the discontinuous positions can begin now.
…diff_update_neutral_position
…discontinuous reconstruction
of find_neutral_surface_positions_discontinuous routine. Will need to think about how this works and possibly rewrite the algorithm to loop over layers instead of by neutral surface index. ToDo: -Think about logic in the discontinuous reconstruction neutral surface algorithm, might need to be kl<=2*nk. Also what it means if neutral surface connects at a discontinuity Changes: -Some defined dimensions are now defined based on the total number of neutral surfaces, set at runtime to be 2*nk+2 if using a continuous reconstruction or 4*nk if using discontinuous
find_neutral_surface_positions. Does not give expected results in tracer_mixing test case, more debugging needed.
still not as expected. Suspect a wrong value when setting position in either right or left search.
…, so need to debug. A problem for future Shao
After an examination, the odd behavior is due to the piecewise constant reconstruction of T and S in the top and bottom layers. ToDo: -Implement PPM boundary extrapolation either by replacing the current T/S reconstruction with the functions in ALE/PPM_functions.F90 or port over the code so that neutral diffusion code is self-contained. -Do something other than linear interpolation of the change in density within a layer to find the neutral surface position. neutral surface.
…ff_update_neutral_position
adcroft
left a comment
adcroft
left a comment
There was a problem hiding this comment.
This looks good.
- There's one little bug. Otherwise it looks correct.
- Unit tests are needed.
- We need to work through each of the thought experiments and build a unit test for each. Some might give the same answer as the continuous (like "left column everywhere denser than right"). Most I think we'll need to do the math on.
- Unit tests should precede debugging a 3d problem.
- Does it run in 3d?
(I'm clicking comment rather than approve since you said review-only).
| real, dimension(:,:,:), allocatable :: dRdS ! dRho/dS (kg/m3/ppt) at interfaces | ||
|
|
||
| if (CS%continuous_reconstruction) then | ||
| if (.not. ALLOCATED(dRdT)) ALLOCATE(dRdT(SZI_(G),SZJ_(G),SZK_(G)+1)) |
There was a problem hiding this comment.
Is the _init() routine called more than once? I'm wondering why the if (.not. ALLOCATED(...))?
There was a problem hiding this comment.
Whoops, this was an oversight before making the PR. I think when I wrote it I was trying to be super-safe that there was no memory leak.
| endif | ||
| call PPM_reconstruction(G%ke, h(i,j,:), S(i,j,:), S_i(i,j,:,:), ppoly_coefficients) | ||
| if (CS%boundary_extrapolation) then | ||
| call PPM_boundary_extrapolation(G%ke, h(i,j,:), S(i,j,:), S_i(i,j,:,:), ppoly_coefficients) |
There was a problem hiding this comment.
CS%boundary_extrapolation is false right now. Extrapolation obviously introduces new extrema in the buoyancy profiles but I think the stage where we limit the flux to be down gradient should stop tracers from forming new extrema. This assertion needs to be checked. Does extrapolation not work (wondering since you hard coded it)?
Talk to @Hallberg-NOAA about how the limiter should work and whether it will zero fluxes for extrapolated profiles.
There was a problem hiding this comment.
That's a good point that I hadn't thought of. My main reason for wanting to do some kind of extrapolation is that it seems somewhat unphysical that you would have a large jump in the sum of hEff in the case between two columns with identical columns versus the case where one is only very slightly lighter or denser than the other because you exclude the entirety of the top and bottom cells.
| do j = G%jsc,G%jec ; do I = G%isc-1,G%iec | ||
| if (G%mask2dCu(I,j)>0.) then | ||
| call neutral_surface_flux(nk, h(i,j,:), h(i+1,j,:), & | ||
| call neutral_surface_flux(nk, CS%nsurf, h(i,j,:), h(i+1,j,:), & |
There was a problem hiding this comment.
Why aren't you masking anymore?
There was a problem hiding this comment.
Mistake, I'll add it back in.
| hR = absolute_position(nk,Pr,KoR,PoR,k_surface) - absolute_position(nk,Pr,KoR,PoR,k_surface-1) | ||
| if ( hL + hR > 0.) then | ||
| hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) | ||
| hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean of layer thicknesses |
There was a problem hiding this comment.
I'm not sure why I didn't write this as 2. * hL * hR / ( ( hL + hR ) + h_neglect ) that would avoid the conditional. I'm not sure if doing so now would change answers.
|
|
||
| ! Identical columns | ||
| call find_neutral_surface_positions(3, & | ||
| call find_neutral_surface_positions_continuous(3, & |
There was a problem hiding this comment.
We'll need units tests for each or your _discontinuous functions.
There was a problem hiding this comment.
As you noted some time before, the discontinuous reconstructions don't always lend themselves nicely to whole number results. That being said, it's on my to-do list before submitting an actual PR.
|
|
||
| !> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns | ||
| !! of combined interfaces using intracell reconstructions of T/S | ||
| subroutine find_neutral_surface_positions_discontinuous(nk, Pres_l, Tl, Sl, dRdT_l, dRdS_l, & |
There was a problem hiding this comment.
It's a shame we can't reuse the same algorithm, or combine them. Debugging these sorts of algorithms is painful. Is this the same algorithm but with different data structures being passed in?
One test to try is to set the discontinuous edge values to the interface values of the "continuous" approach. It should give the same solution as the continuous but with 2*G%ke vanished sub-layers.
There was a problem hiding this comment.
I spent a bit of time trying to shoehorn your original algorithm into the discontinuous one, but to no avail. What it boiled down to was that the original algorithm relies on the fact that the layer index always corresponds to the top interface. With the discontinuous reconstructions a different mapping between interface and layer exists, plus you may potentially need to search across the discontinuity.
I did do something similar to the test you suggested when I was prototyping the logic for the new one for some really simple cases that produced qualitatively similar results (didn't check for bitwise reproduction). If they do, then I agree that the code could be collapsed so that there's not a huge bit of similar code.
| if (k<=G%ke) then | ||
| Pint(:,j,k+1) = Pint(:,j,k) + h(:,j,k) * GV%H_to_Pa ! Pressure at next interface, k+1 (Pa) | ||
| call calculate_density_derivs(T_i(:,j,k,2), S_i(:,j,k,2), Pint(:,j,k+1), & | ||
| dRdT_i(:,j,k,1), dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, EOS) |
There was a problem hiding this comment.
Bug: there's a 1 here that should be a 2.
There was a problem hiding this comment.
PS> For the record, the use of indices 1,2 for left,right in the reconstruction code is not my doing and I don't like it because you can't read the code. One day we'll clean up the reconstruction interfaces...
There was a problem hiding this comment.
Thanks for the catch. Also, I agree that it's a little bit difficult for 1,2 to map onto left right edge values. At the same time the solution that I can think of off the top I my head is to have separate variables for the two of them which can lead to a huge amount of input arguments (which was the original I had this algorithm coded).
|
Many thanks for the review. I'll close this PR for now and reopen when it's ready to be looked at again. |
…fusion is happening. Maybe in reconstructions?
implement a better way of determining the neutral surface position based on the polynomial reconstructions of T/S.
with the polynomial reconstructions
… from the linear and TEOS-10 equations of state
…ition. Need to test
combination of bisection and Newton's method with regula falsi or Brent's method
Brent's method to do neutral position finding. Need to set it up so that if either second derivatives of rho aren't available or if the Newton's step goes out of the interval [0,1], it switches to Brent's method.
derivatives. Also, wrote an extractor routine so that members of the EOS type are available outside of the MOM_EOS module
-The refinement of the neutral position now includes a different implementation of Brent's method. -Also, the refinement is now forced to be bracketed between the last position within the layer and the bottom interface. -Update the flux calculation so that the continuous case calculates the top and bottom bounding neutral surfaces using the parabolic reconstructions instead of linear interpolation. To Do: -Try to combine the discontinuous and continuous find_neutral_position routines
Maybe think about doing the same direction check regardless of reached_bottom.
…uous neutral position finding
-Rearrange and make optional arguments for find_neutral_surface_positions_discontinuous if refine_position is not requested -Simplify some of the logic for search_other_column
…ff_update_neutral_position # Conflicts: # src/equation_of_state/MOM_EOS_TEOS10.F90
seems to pass all tests. Performance in more extreme configurations of tracer_mixing work. All unit tests for discontinuous have been updated.
|
The new neutral diffusion code is ready to be pulled into the main code. For now, the original reconstruction code is retained and remains the default so answers should not change. The improvements made here (discontinuous reconstruction and iterating to find a true neutral position) increases the computational cost of tracer diffusion by a factor of 3-4 in OM4_05, about the same cost as the advection part of the routine. Unit tests for a wide variety of cases are included. |
ashao
changed the title
|
There are a couple of the density derivative routines that are being duplicated for scalar and array arguments, including some with messy calculations. I would prefer that the pattern of calculate_density be followed, in that the scalar versions are not coded separately but rather copy their arguments into 1-point arrays and then call the array version. Doing it that way instead would avoid having two separate (and potentially inconsistent) implementations of the algorithms. |
Based on the suggestion of @Hallberg-NOAA, the array and scalar versions of the EOS functions were wrapped into interfaces. For each of the available equations of state, the array and scalar versions were also put behind interfaces. For the Wright EOS, the first and second deriviative scalar variants of the routines call the array variant by promoting the scalar T, S, and P to 1-element arrays and demoting the 1-element output arrays to scalars. Answers not expected to change.
Previous commit failed Travis test because of compiler-dependent evaluation of logical expressions. (if (kl>1 .and. bot_connected(kl-1)) could yield to an out of bounds error if the compiler evaluates bother sides of the statement. The offending line was moved one level below an explicit check for (kl>1)
|
The lengthy list of commits in this PR appears to change both answers and (less concerningly) the parameter_doc files in several test cases. However, I do not see anything in the commit logs that indicate which commits actually change the answers. If there are answer-changing modifications, please break up this pull request into the larger part that does not change answers and a very limited set of code changes that do change answers. |
|
None of the commits should have changed answers and the cases checked by Travis passed. To help narrow down where I may have inadvertently introduced a change, could you send a subset of the test cases which failed their checks? |
|
Oh wait, I see changes that shouldn't have been there. Submitting yet another commit shortly |
Changes made in #a49373c are correct in that the continuous polynomials should be evaluated to calculate gradients. However, they were intended to be implemented in a different set of commits. Answers should now be unchanged.
3 participants
This PR adds the option to connect neutral surfaces using a discontinuous reconstruction of temperature and salinity (after drawing 100s of arrows and surfaces on paper, writing new logic, throwing it away, drawing more arrows, and finally coming up with a consistent algorithm). A notebook demonstrating the algorithm for some simple test cases can be found at: https://github.com/ashao/MOM6_analysis/blob/master/neutral_diffusion/neutral_diffusion_examples.ipynb. This is controlled by the runtime parameter NDIFF_CONTINUOUS which is by default True if using the original continuous parabolic reconstructions and false if doing a discontinuous PPM reconstruction.
The main deficiency of this algorithm as it stands is that T/S in the top and bottom boundary layers are using a piecewise continuous reconstruction. This leads to some odd behavior in the two bottom-most and topmost layers in the ocean_only/tracer_mixing experiment. If you have two columns with a density profile rho(z) and the other with rho(z)+epsilon, the entire top or bottom layer will be skipped whereas in actuality some part of the top and bottom layers should always overlap. However, in the interior the neutral surfaces seem to behave as expected.
Before finalizing this update to the neutral diffusion code, two more improvements could be made:
-Boundary extrapolation (to address the above deficiency) applied to both the T/S reconstructions but also the tracer being diffused.
-Something other than a linear interpolation could be used to find where the neutral surface connects within a layer.