[WIP] add KPP mixing - #512
Conversation
|
This PR is rebased on #461 so it has extra changes. Anything related to forcing can be ignored here |
There was a problem hiding this comment.
Nice job on this @vanroekel! This must have been a lot of work.
This was just my first pass as this was a lot of code to go through, but its a place to start. I think in addition to the specific comments I made below, some commenting throughout the KPP guts would be helpful for me to review more in depth. There are a lot of variable names that are ambiguous, so I am not sure what is being calculated. I haven't done the job of going back and forth between the CVMix code and this code to make sure its all of the guts are exactly correct. I think the single column testing might actually do a better job of catching those or letting us know there is a problem with them in the first place. More so I looked at how everything (functions, variables, etc) was being used and passed around from a higher level.
Also, a lot of the variable names don't conform to the Omega standard. I'm guessing because they were direct variable name ports from CVMix?
Lastly... I think @mwarusz input on how things are called from a computational perspective will be super helpful here as it was a lot for me to go through and I am sure I did not catch those things.
| Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); | ||
|
|
||
| Array1DReal SurfacePressure("KPP-SurfacePressure", NCellsAll); | ||
| deepCopy(SurfacePressure, 1.0e5_Real); |
There was a problem hiding this comment.
Does this override the Fill Values?
| const Real ct_freezing = | ||
| Eos::calcCtFreezing(LocEosChoice, surface_salinity, | ||
| PressureMid(ICell, KSurf) * Pa2Db, 0.0_Real); | ||
| const Real heat_flux = |
There was a problem hiding this comment.
We will probably want to update these once the final equation is decided upon in #461.
Also, all of this is already computed in TendencyTerms.h by #461... can we use the temp flux and salt flux from there instead to save computation and also to lessen the risk of these calculations not matching if say the thermo forcing assumptions are updated in the future, then they don't have to be updated in multiple places? Maybe we can save surface temp flux and salt flux as variables calculated in TendencyTerms.h that can be passed in to here?
There was a problem hiding this comment.
Yes, @alicebarthel and I discussed the same things. I have a potential change to share code now that I'm posting soon here
There was a problem hiding this comment.
What was the reason for moving this to the .cpp side? And was the equivalent code in VertMix.h removed?
There was a problem hiding this comment.
This was a mistake in a copilot PR, I'm working to fix this now
There was a problem hiding this comment.
Actually what is 'this'? I see some computation here but I think that was added in the tendency hookup not in this PR. I just added to it.
| for (int KVec = 0; KVec < KLen; ++KVec) { | ||
| const I4 K = KStart + KVec; | ||
| const bool ApplyConv = | ||
| (!LocKPPEnabled) || |
There was a problem hiding this comment.
Maybe I'm not understanding something, but why is this logic applied to just convective adjustment and not the shear mixing part too?
There was a problem hiding this comment.
This is an assumption made in KPP that convective diffusivity should not be counted in the BL since it roughly resembles the non local flux, so there would be double counting. I should add this to the docs.
StopOBLSearch only rescaled the bulk Richardson crossing threshold, it never terminated the OBL search. At its default of 1.0 the threshold equals CriticalBulkRichardsonNumber exactly, so the option duplicated an existing knob. Decoupling the two is also physically inconsistent: RiCritical normalizes the unresolved shear Vt^2 in Large et al. (1994) Eq. 23, which is derived assuming the same critical Richardson number used for the crossing test. Use CriticalBulkRichardsonNumber to adjust the OBL depth criterion instead. This is answer preserving since every configuration used the default value.
The bulk Richardson search reset the surface-layer reference averages to the top of the column at every trial depth and then replayed the accumulation downward, making the search O(K^2) per column, and O(K^2*nEdges) for the per-edge velocity averages. The averaging window spans the top epsilon*d, which grows monotonically with the trial depth, so the window end only ever moves downward. Hoist the running sums and their level pointers out of the trial-depth loop so the accumulation becomes a two-pointer scan that is amortized O(1) per level. This is bit for bit identical since the accumulation order from the surface is unchanged, it is simply no longer repeated.
The bulk Richardson search traversed every model layer even after the crossing was found. Only the crossing level and the two above it feed the quadratic refinement of the boundary layer depth, so the remaining levels contributed nothing to the diagnosed depth or to any mixing coefficient. Break out of the loop once the crossing is latched. The four Ri diagnostic fields written inside the loop, BulkRichardsonNumber, BulkRichardsonShear, UnresolvedShear and BuoyancyJump, are consequently zero below the boundary layer base; they are zero filled before the kernel, so the values remain well defined. Setting DebugDiagnostics suppresses the early exit and restores the full water column profiles for those fields. BoundaryLayerDepth and all mixing coefficients are unchanged.
|
This KPP branch works great on the SOMA test case! SOMA after 4 months, 32km resolution. RK4 with 2 min timestep, viscosity del4 of 2e11. No explicit tracer diffusion. The instability creeps northward from the southern boundary where cooler water is being advected upwards. KPP saves the day! Using this branch with linear EOS on both (but JM/TEOS10 looks the same) I used default KPP flags in both models. Thank you @vanroekel for this addition! Thanks to @hyungyukang for his consultation to get this working. |
|
You should remove these. They are printing with every RHS compute |
Factor repeated KPP device math into KPPConstants.h helpers, decouple the non-local flux shape from the coefficient matching option, and correct several array extents that made KPP abort on multi-rank runs. Shared device math (answer-preserving): - Add kppTurbScales, kppMatchShape, kppNonLocalCs, kppClampOBLDepth and kppOBLIndex to KPPConstants.h. - Replace the six duplicated blocks in KPPMix.cpp: two turbulent-scale computations, two match-shape computations, the inline C_s constant, and two OBL clamp plus index-search loops. Collapse the repeated UseInteriorMix && MatchBoth predicate into LocUseMatchedShapes. Non-local shape (changes answers under MatchBoth only): - Always build the non-local flux from the unmatched scalar shape rather than reusing the diffusivity shape. The matched shape is non-zero at the boundary layer base by construction, so gamma was finite there and then dropped discontinuously to zero immediately below. CVMix and MPAS-Ocean likewise treat matching as a diffusivity choice, separate from the non-local shape. - Update testMatchBothInteriorCoefficients, which previously asserted both the matched value inside the boundary layer and zero at its base. Those two expectations cannot both hold. - Update design documentation, which had described the shared shape as intentional. Array extents (answer-preserving): - Size the KPP scratch arrays in Tendencies::computeKPPFields and the boundary layer index copy in VertMix::computeVertMix with NCellsSize rather than NCellsAll, matching the KPPMix members they are copied to and from. Kokkos rejects mismatched deep_copy extents, so KPP aborted on any run where NCellsAll and NCellsSize differ. Single-rank tests could not detect this. - Size the KPPMixTest scratch and expected buffers the same way. Add unit tests for the new helpers. Full Omega suite passes 55/55 on pm-cpu with gnu.
5 participants
This PR adds the mixing scheme of Large et al (1994), the K-Profile Parameterization (KPP) to omega. It includes the local and non-local components and theory wave langmuir mixing. The code
The bulk of the code is written through co-pilot with physics guidance and corrections by me.
Checklist
Documentation:
Linting
Building
Testing
aurora, oneapi-ifx, mpich
chrysalis, oneapi-ifx, openmpi
frontier, craygnu, mpich
frontier, craygnu-mphipcc, mpich
pm-cpu, gnu, mpich
pm-gpu, gnugpu, mpich
Provide relevant details in a comment to the PR titled
Testingwith the following:have been run on and indicate that are all passing.
has passed, using the Polaris
e3sm_submodules/Omegabaseline-pfor both the baseline (Polarise3sm_submodules/Omega) and the PR buildPerformance related PRs: Please include a relevant PACE experiment link documenting performance before and after.
New tests: