Repository navigation
add KPP mixing - #512
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.
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.
|
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 |
6ccf433 to
de1c24d
Compare
|
@vanroekel, can you check this? |
|
Claude says that your results are consistent with running clang-format 22.x but that's not what we use in the omega dev environment, we're pinned to 18.1.x: |
|
The next CI issue is related to the docs. |
|
@xylar thanks for looking into this. I definitely thought was the version. I had regenerated with 18.1.8 and it still showed no issue, but just did it again and I see the fail in lint. Should I fix that here even though it's not touched in this PR? or should I make a separate PR to fix the linting issue? |
|
I'm happy to do either. |
|
@vanroekel, if there's a lint problem here and it isn't from a file you charged, a rebase must fix it. We don't have lint problems on develop. I'll give it another look over the weekend. |
5b7db4f to
287a9e4
Compare
|
@vanroekel, I see what you mean. You changed the time stepping file but not that part. How strange! We should fix this elsewhere first. |
|
Okay very weird, but after a rebase, those lines do show up as changed by claude in one of the KPP commits. I don't understand why or how, but I went ahead and fixed the linting and pushed |
mwarusz
left a comment
mwarusz
left a comment
There was a problem hiding this comment.
This is a big PR and I am still digesting it, but I have a few initial comments.
I also have been thinking how to make the time stepping sequence a bit better for parametrizations that are called once per step, either at the start or end of the step . My rough idea is to refactor doStep to do the following:
-
doStepStart- common for all time steppersprescribeState- compute all auxiliary vars
- compute KPP and maybe submesoscale eddies if we decide to compute it only once, like MPAS-O does
-
doStepMain- depends on time stepper- mostly the same as before, except, in the first stage, it can assume that all auxiliary variables have been computed, and hence it can use the
computeXTendenciesOnlyfunctions
- mostly the same as before, except, in the first stage, it can assume that all auxiliary variables have been computed, and hence it can use the
-
doStepFinish- common for all time steppers- implicit vertical mixing
@hyungyukang
What do you think of this plan ?
| // Called once per time step by the active time stepper. | ||
| void computeKPPFields(const OceanState *State, | ||
| const Array3DReal &TracerArray, int ThickTimeLevel, | ||
| int VelTimeLevel); |
There was a problem hiding this comment.
Same here, why isn't this just a method of the KPP class.
There was a problem hiding this comment.
okay I think i got claude to finally move it. It's now in KPPMix.h as updateKPPFields
|
I had Claude push a fix, so that's why. |
|
@mwarusz and @katsmith133 this should be ready for review again. I think all comments are addressed. I'm going to start testing again on the various machines. |
|
@katsmith133 - I still have one comment left from you I wasn't sure what you were asking about in it. |
mwarusz
left a comment
mwarusz
left a comment
There was a problem hiding this comment.
I have a few more comments. I also reopened two of my previous comments, since they don't seem to be resolved. Moreover, the current version of this branch doesn't compile for me. @vanroekel are you sure you pushed all of your fixes ?
One general issue I noticed is that Claude (I assume) is adding a lot of checks for conditions that can never happen on a valid Omega mesh. I highlighted some of them, but there are probably more. Most of the time it is hard to figure out which checks and index clamps are actually necessary.
| // Edge weights are the MPAS kite areas (0.25*dc*dv) normalized by | ||
| // the cell area. |
There was a problem hiding this comment.
I think this comment is misleading. 0.25*dc*dv is the area of triangle formed by cell center and edge vertices. "kite areas" makes me think of KiteAreasOnVertex, which the code is not using.
| const I4 NEdges = NEdgesOnCell(ICell); | ||
| const I4 NEdgesEff = Kokkos::min(NEdges, I4(MaxEdgesBound)); |
There was a problem hiding this comment.
If MaxEdgesBound is a bound on the number of edges then NEdgesEff == NEdges always, and we don't need both variables.
There was a problem hiding this comment.
I agree, removed
| if (NEdgesEff > 0) { | ||
| Real SumW = 0.0_Real; | ||
| for (I4 J = 0; J < NEdgesEff; ++J) { | ||
| if (EdgeValid[J]) { | ||
| SumW += EdgeWeights[J]; | ||
| } | ||
| } | ||
| const Real InvSumW = 1.0_Real / SumW; | ||
| for (I4 J = 0; J < NEdgesEff; ++J) { | ||
| if (EdgeValid[J]) { | ||
| EdgeWeights[J] *= InvSumW; | ||
| } | ||
| } | ||
| } |
There was a problem hiding this comment.
I don't think we want to do this. I believe that 0.25_Real * DcEdge(IEdge) * DvEdge(IEdge) * InvAreaCell is already the correct normalized weight for edge to cell interpolation.
| if (NEdgesEff > 0) { | |
| Real SumW = 0.0_Real; | |
| for (I4 J = 0; J < NEdgesEff; ++J) { | |
| if (EdgeValid[J]) { | |
| SumW += EdgeWeights[J]; | |
| } | |
| } | |
| const Real InvSumW = 1.0_Real / SumW; | |
| for (I4 J = 0; J < NEdgesEff; ++J) { | |
| if (EdgeValid[J]) { | |
| EdgeWeights[J] *= InvSumW; | |
| } | |
| } | |
| } |
There was a problem hiding this comment.
yes this is odd. looks like it is computing the sum of the weights which is weird given the invAreaCell, this is removed
| /// Compile-time upper bound on MaxEdges, for sizing on-stack (per-thread) | ||
| /// arrays in device kernels. Meshes with more edges per cell are rejected. | ||
| static constexpr I4 MaxEdgesBound = 10; |
There was a problem hiding this comment.
Despite what the comment says, I don't see a check for this in HorzMesh. There is one in tracer advection, but only if this tendency term is enabled. Let's move this check into HorzMesh.
There was a problem hiding this comment.
I believe this should be moved now
| const Real InvAreaCell = | ||
| 1.0_Real / Kokkos::fmax(AreaCell(ICell), KPP::Tiny); |
There was a problem hiding this comment.
I don't think any valid mesh has cell areas smaller than 1e-20.
| const Real InvAreaCell = | |
| 1.0_Real / Kokkos::fmax(AreaCell(ICell), KPP::Tiny); | |
| const Real InvAreaCell = | |
| 1.0_Real / AreaCell(ICell); |
There was a problem hiding this comment.
yes this is crazy, this guard is removed.
|
|
||
| // Resolved shear |V_r - V(d)|^2, averaged over the cell edges | ||
| Real DeltaVSq = 0.0_Real; | ||
| if (NEdges > 0) { |
There was a problem hiding this comment.
A valid mesh cannot have cells with zero edges.
| if (NEdges > 0) { |
There was a problem hiding this comment.
similarly removed
| if (INeighbor == NCellsAll) { | ||
| continue; | ||
| } | ||
|
|
||
| const I4 KMinNbr = MinLayerCell(INeighbor); | ||
| if (KMinNbr < 0 || KMinNbr >= NVertLayers) { | ||
| continue; | ||
| } |
There was a problem hiding this comment.
I believe the first check is enough to handle horizontal boundaries.
| if (INeighbor == NCellsAll) { | |
| continue; | |
| } | |
| const I4 KMinNbr = MinLayerCell(INeighbor); | |
| if (KMinNbr < 0 || KMinNbr >= NVertLayers) { | |
| continue; | |
| } | |
| if (INeighbor == NCellsAll) { | |
| continue; | |
| } |
There was a problem hiding this comment.
yes, not needed for a 2-D field smoothing, removed
| if (KMin < 0 || KMax < KMin || KMin >= NVertLayers) { | ||
| return; | ||
| } |
There was a problem hiding this comment.
Don't think this can ever happen.
| if (KMin < 0 || KMax < KMin || KMin >= NVertLayers) { | |
| return; | |
| } |
There was a problem hiding this comment.
Yes, agreed. Vertical checks also not necessary for 2D field routines.
| const bool IsValidCell = | ||
| KMin >= 0 && KMin < NVertLayers && KMax >= KMin; | ||
| if (IsValidCell && K >= KMin && K <= KMax + 1) { |
There was a problem hiding this comment.
Same here.
| const bool IsValidCell = | |
| KMin >= 0 && KMin < NVertLayers && KMax >= KMin; | |
| if (IsValidCell && K >= KMin && K <= KMax + 1) { | |
| if (K >= KMin && K <= KMax + 1) { |
@mwarusz , Sorry for the delayed response. I think that's a good idea. It would make it clearer where the physical parameterizations are applied in the time-stepping sequence, and I think it would also make the overall code structure cleaner. This is not directly related to this PR, but during the refactoring, could we also address the time-lagged auxiliary variable issue? We could move the auxiliary variable computations to
This is just a rough idea for now. We can discuss it further in another channel later. |
|
Sorry @mwarusz something must have gone wrong in the push. I had seen computeKPPFields moved previously (and it compiled on pm-cpu/gpu). I'll look into it. yes Co-pilot is very keen in protecting against things that never happen. I removed a number of them but see more to remove. |
|
I just merged #517. There are a number of conflicts between that branch and this one. I had previously tested KPP+SE Time Stepping using the merged branch remotes/hyungyukang/omega/se-time-stepper-kpp, (made by @hyungyukang). When you rebase this KPP branch, that will be a useful resource when you resolve conflicts. |
|
Okay I think I addressed the comments. I also worked with copilot to remove unneeded guards, such as those guaranteed by a valid mesh, no vertical bounds checking in 2D loops, or cases that don't happen (e.g. AreaCell < 1e-20). I'll look through the code more tomorrow to check for more places to streamline The code builds okay for me on pm-cpu and pm-gpu |
|
With the current head 9c25f68, I get two cTest failures on perlmutter cpus with gnu. The driver test detects state validation failure after the first time step and I get this for the KPP OSBL test Here is my command sequence:
|
|
Perlmutter GPU cTests failed identically to above, with same error messages. |
Previously MatchBoth silently passed back to SimpleShapes. It is now implemented
51614ff to
8c56ca4
Compare
|
Local merge passes CTests and Polaris
|
|
Thanks @sbrus89! |
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
New tests: