New turbulence closure option AMD

[cbegeman]

This PR adds a new option for turbulence closure scheme, Anisotropic Minimum Dissipation (AMD).

• AMD is implemented following Abkar et al. 2016 doi:10.1103/PhysRevFluids.1.041701, which extends the scheme to include scalar diffusivities.
• In this scheme, scalar diffusivities differ depending on the scalar (and its gradients). Thus, there is a new output variable for the saline diffusivity and a new diffusivity input to diffusion_s (the subroutine that computes scalar diffusion tendencies).
• This PR also includes an unrelated commit that saves the current commit hash to the build directory

[cbegeman]

@vanroekel @xylar

So far I've run short test cases and compared diffusivity values with the Moeng-Wyngaard (MW) scheme (diffusivities near the boundary are orders of magnitude higher for AMD than MW but comparable elsewhere). I'll run longer test cases next week. You can wait to review until those results are available.

Timing is similar for both schemes; just about all the time we save with AMD by not solving a prognostic equation for SGS TKE is used computing the diffusivities.

[vanroekel]

@cbegeman please let me know when you'd like me to do a review here. Happy to dig in when you are ready. This will be a great feature for us to have.

[cbegeman]

@vanroekel, thanks. This is ready for your review.

[cbegeman]

@katsmith133 Are you interested in weighing in on this PR or taking on the role of reviewer?

[katsmith133]

@cbegeman Happy to try to help with this, though I'll admit I am not up to speed on the ins and outs of the AMD closure. I can do some reading though and see if I can then provide some help.

[cbegeman]

@katsmith133 Totally up to you whether you want to weigh in. No outside reading required.

[katsmith133]

@cbegeman I've never been a reviewer, so not sure exactly what that entails... But in terms of the km issue at hand, my first instinct would be the km_max option because, like you said, its more transparent and obvious when reached.

[cbegeman]

@katsmith133 Sure, I can speak a little bit to the norms of the PALM group. In this case, the reviewer would need to look at the source paper for the scheme to check my implementation, and I'm planning to post some figures from the test case for comparison. Depending on complexity of the PR, the reviewer has typically also checked out the branch and run a case themselves. Since @vanroekel has offered to review, you don't have to do all this. You or anyone in the group is welcome to post comments/recommendations/requests for clarification.

[katsmith133]

@cbegeman Ah, ok. That sounds like something I can easily do. I will take a closer look and put my 2cents in where I can then.

[vanroekel]

@cbegeman I'm working through this PR now, but had a deeper broader question to pose. As I track from Rozema et al to Vreugdenheil and Taylor I'm a bit concerned by having a ks and kh. The arguments Abkar et al are for a passive scalar which I'm not convinced applies to T&S directly. Physically I think we want to dissipate buoyancy variance not that in T&S so I "think" (not sure) we want kappa_sgs to be based on buoyancy not have one for T&S. I think of it this way, if you look at eq 1 and 6 of V&T 2018 -- the potential temperature appears here, but in equation 1 it is clear this is because it is the buoyancy term. So for T&S I would expect equation 1 to go from

g*alpha*Theta

to

g*b

Suggesting that b would appear in equation 6. I would suggest the same for kappa_sgs, but maybe this isn't correct. What do you think?

How different do you see kh and ks?

[vanroekel]

let me know if you'd like to discuss my question further offline via Webex of phone.

[cbegeman]

@vanroekel Equation 1 and 6 of V&T 2018 both pertain to the solution of km, and the term of which you speak does depend on buoyancy as I've implemented it.

Now, buoyancy does not appear at all in the solution of kh and ks. Could you clarify further if you think that buoyancy should appear there? A relevant equation is Equation 10 from Abkar and Moin (2017).

km and ks can be quite different depending on the degree of stratification and variance in each.

[vanroekel]

@cbegeman thinking further about my bigger comment and finishing a read through of the code, I think what you have is in the spirit of the AMD literature and we should leave as you have it. I would still like to see how different kh from ks is in your runs. Could you post a plot here?

[vanroekel]

@cbegeman sorry didn't see your comment. I had made my original question prior to finishing a read of your code. After reading I think you're right and we shouldn't change. Sorry for the premature comment!

[vanroekel]

Okay, I've finished going through everything. This was a herculean effort. Great work on this @cbegeman. Outside of a few cleanups and the issue with the tracer gradient, this looks great. I will try and do a test here as well, do you have a case I could test this with? I'll also try one of my past cases to verify nothing changes.

[cbegeman]

Okay, I've finished going through everything. This was a herculean effort. Great work on this @cbegeman. Outside of a few cleanups and the issue with the tracer gradient, this looks great. I will try and do a test here as well, do you have a case I could test this with? I'll also try one of my past cases to verify nothing changes.

I'll post a case shortly, but thought I should note that if you try it with a past case, be sure to compare after 5941195, for turbulence_closure = 'mw'. I've found it can make a difference.

[cbegeman]

@vanroekel Thanks for your patience. Here is the atmospheric test case: /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_atm_amd_stableflux_disturb-2 and an ocean case: /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_amd_mcphee_onesided

[cbegeman]

@katsmith133 In case it wasn't clear, this is ready for you to run whatever test case you desire. Note though that the combination of stratification and low tke can lead to small timesteps restricted by diffusion (high km). I've run into this problem in the rayleigh damping region, so I haven't been imposing (initial) stratification there.

[katsmith133]

@cbegeman testing it today! Will get back to you guys soon.

[cbegeman]

@cbegeman thinking further about my bigger comment and finishing a read through of the code, I think what you have is in the spirit of the AMD literature and we should leave as you have it. I would still like to see how different kh from ks is in your runs. Could you post a plot here?

Here are a few images of what I have so far. I'm still experimenting with the disturbance settings for melting cases to kick-start TBL formation a bit harder. As you can see, kh and ks are very similar.

[vanroekel]

thanks for the plots @cbegeman it is good to know Kh and Ks remain close.

[katsmith133]

@cbegeman I've run a couple of simple convection tests, comparing results between the master and amd_closure branches, with AMD turned off. With AMD turned off, they should produce the exact same results, except that there are noticeable differences in the momentum fields and minor differences in the temperature and salinity fields:

After talking with @vanroekel, I've pinned the difference down to the changes made in init_grid.f90 file. In the amd_closure branch these lines were changed such that there is no longer an if-statement: https://github.com/xylar/palm_les_lanl/blob/c433fb5152ee03bf4bd574f167d99a98ad621ab2/trunk/SOURCE/init_grid.f90#L1892-L1894 If I put the if-statement back in, then the results get a bit better (note: the errors are an order of magnitude smaller, so the color bar is also an order of magnitude smaller):

Then when I replace the init_grid.f90 file in the amd_closure branch completely with the version that is in master, the errors are completely gone:

I haven't pinned down what other lines in init_grid.f90 are causing the issue yet. Working on that now, but perhaps you know off-hand which ones they might be.

[katsmith133]

@cbegeman Just to clarify, when the legend says Master and AMD Closure, it is referring to the branch. AMD was not turned on in any of these results.

[cbegeman]

@katsmith133 Thanks for looking into this! I made a few changes to wall_flags_0, one of which is through the topo variable you pointed out. Let me remind myself which of these are necessary and which could be changed only for the amd case.

[cbegeman]

@katsmith133 Can you point me to the namelist file for this test case? Thanks!

[katsmith133]

@cbegeman It is just the test_oceanml_p3d file in the master repository. I added a gradient in salinity to it, but you can see the errors even without that.

[cbegeman]

@katsmith133, I reverted the changes in init_grid.f90, with some minor adjustments. The key issue was whether to solve for velocity and scalars at the bottom (k=nzb) grid point when there is no topography (constant_flux_layer = 'none'). I had set conditions such that PALM did solve at nzb, but I have now reverted those changes. I did run a short test to make sure there is a bit-for-bit match with master.

[katsmith133]

@cbegeman and @vanroekel I am getting the same matching with master as well with this fix. Other than that, I ran the case that Carolyn provided here: /lustre/scratch3/turquoise/cbegeman/palm/jobs/test_ocean_amd_mcphee_onesided. It runs and produces reasonable results as far as I can tell.

[cbegeman]

Thanks @katsmith133 for the testing! @vanroekel I think I addressed all your comments now with a few aesthetic changes, but let me know if I missed something.

[xylar]

Congratulations!!