Climatemachine.jl: Taylor-Green Vortex oscillations

[slifer50]

Description

Running the Taylor Green vortex case yields unphysical behavior in terms of the dissipation rate of kinetic energy. There are oscillations where we should be seeing a smoothe curve. The two plots below present the dissipation rate, the first is what I have been getting with CLIMA where oscillations are visible, the second was obtained with NUMA. The data for both of these plots were processed in the exact same way using the exact same scripts. @lcw @jkozdon @fxgiraldo

Steps to reproduce

Steps to reproduce the behavior:

The case can be found here https://github.com/CliMA/ClimateMachine.jl/blob/yt/tg_runs/experiments/AtmosLES/taylor-green.jl. If anyone would like to reproduce the results I can send you the post-processing scripts.

Expected behavior

We should be getting a smooth curve similar to the one produced by the data from NUMA.

For CLIMA developers

• [ ] The bug can be reproduced
• [ ] Labels have been added and developers assigned
• [ ] The bug can be fixed (if not, please explain why not in a comment below)
• [ ] There is a timeline to fix the bug
• [ ] The bug has been fixed (Please link the PR)

[tapios]

Why are the NUMA dissipation rates 2 orders of magnitude smaller? Are the domain sizes the same?

[smarras79]

@slifer50 following up on what Tapio asked: in the plot that came through with the issue notification email you plot the same for CLIMA but the the values are very different from the ones plotted above also for CLIMA. Are you sure you are plotting the same quantities?

[slifer50]

@tapios @smarras79 Sorry for not clarifying that, everything is the same. There is a division by 100 that is not present for the CLIMA plot shown above that is done for NUMA plot for the sake of non-dimensionality. I'll update the the post with the corrected CLIMA plot.

[slifer50]

@tapios @smarras79 The plots should both be scaled the same now. Only Vreman is shown since the results are the same for smagorinsky.

[smarras79]

Given that the results above for numa are CG, the difference can only be in a possible error in the imposition of continuity across elements in our DG implementation. Since this is run with explicit time integration at Courant < 1, the LES should be properly resolved so that time should not be be the cause of error.

DG should be, if anything more dissipative than CG; certainly not less. So, I really think that the issue is in a possible bug in the flux calculation, flux usage, or continuity enforcing across elements.

Another possible problem could be simply due to how the data are written to VTK, but I don't think that's the issue.

cc: @lcw @jkozdon @fxgiraldo @thomasgibson does this diagnosis seem plausible?

[jkozdon]

Why don't you run NUMA with DG then it's an apple-to-apples comparison?

[slifer50]

I've started one and will post results here as soon as I have them.

[slifer50]

@jkozdon @smarras79 @tapios @fxgiraldo @thomasgibson @lcw Here is the same plot using NUMA with DG : No oscillations, we get a smooth decay:

[tapios]

This is NUMA DG with Smagorinsky and same Smagorinsky coefficient as in you CliMA simulation? (So a direct one-on-one comparison?)

[tapios]

And are you using the Lilly stratification correction in CliMA, or is it straight Smagorinsky?

[mwarusz]

@fxgiraldo Does NUMA use LDG to implement higher-order operators ? From a quick look at the NUMA code I saw references to symmetric interior penalty methods. Wondering if that could explain the difference.

[slifer50]

@tapios NUMA is using Smagorinsky with the same coefficient. The Lilly stratification is active in CLIMA.

[smarras79]

I think @mwarusz has the good intuition here; Frank can confirm, but I think Clima and NUMA discretize the diffusion operators with two different methods.

[tapios]

@tapios NUMA is using Smagorinsky with the same coefficient. The Lilly stratification is active in CLIMA.

I suggest you turn off the Lilly factor for a direct comparison. I could imagine relaxation oscillations arising from the interaction of the Ri correction with weak stratification that periodically builds up. The long timescales of the oscillations point in that direction.

[slifer50]

@tapios I'll start a run without the lilly correction and post the results here

[smarras79]

But if gravity is zero there shouldn't be any correction activating though, am I right? this test has zero gravity.

[tapios]

But if gravity is zero there shouldn't be any correction activating though, am I right? this test has zero gravity.

Good point.

[slifer50]

@smarras79 Yes you're right there is no active gravity for this test

[fxgiraldo]

@fxgiraldo Does NUMA use LDG to implement higher-order operators ? From a quick look at the NUMA code I saw references to symmetric interior penalty methods. Wondering if that could explain the difference.

@mwarusz, yes we use SIPDG in NUMA but would be surprised that the difference in how we construct the second order operators is the culprit. Wish we hadn't removed LDG from NUMA in order to compare. Let me see if I can carve out some time to hack it into NUMA for a more apples to apples comparison.

[slifer50]

@mwarusz @jkozdon @fxgiraldo @tapios I've done a test by calculating the Kinetic energy on the DG grid during the simulation in order to verify that this is not a problem with any of the outputting/post-processing. Still seeing the same results as the plot above. The only things I can think of that might be causing this would be have to be the numerics. Any insight in how to deal with this or what I can try to test?

[jkozdon]

The only things I can think of that might be causing this would be have to be the numerics. Any insight in how to deal with this or what I can try to test?

@slifer50 Thanks for taking a careful look at things.

Which NUMA formulation are you using? Are you using the same equation set and DG formulation. (On a box with a straight mesh, I think that NUMA's weak form of DG should be close, though ClimateMachine uses LDG and NUMA does not. Not sure which equation set you should use. @fxgiraldo can probably comment here.) Also, are you using the same time step and time integrator in both codes?

If you think that it's something in the numerics and you think that the codes (ClimataMachine and NUMA) are using identical numerics, then there is actually a difference.

What you could try to do is unit compare each of the codes with known values. You can driver the both codes numerical flux routines with identical values and the SGS models.

I don't really know any of the details about the turbulence models in NUMA or ClimateMachine, but doing direct comparison with identical inputs would reveal the differences.

Another option would be to do a single tendency evaluation of each code and incrementally turn physics back on (though the differences in treatment of 2nd order terms would cause problems here). If you ran for a short time without turbulence with both codes, they should produce the same results with the same time step and time integrator.

You could also run both codes with a constant viscosity which would reveal if the problem is in the turbulence or not.

[slifer50]

Hey @jkozdon Thanks for the input. I actually just ran CLIMA with constant viscosity and got the same behavior as with using the Smagorinsky or Vreman models so this is definitely not a turbulence closure issue, the result I've posted above for NUMA dg was using Smagorinsky which was also used for CLIMA. This does not seem to be a time stepping issue either since I've tried running CLIMA with a few of the different explicit time integrators with smaller and smaller time steps and have observed the same behavior. I'm not quite sure how to go about driving the numerical fluxes for either code. There should be some differences between NUMA and CLIMA, I'm pretty sure CLIMA and NUMA are using two different equation sets. However, the paper I've posted here uses the same equation set as CLIMA with DG and does not seem to exhibit this issue
Murman_S_AIAA_2014_2784.pdf

[smarras79]

@jkozdon @fxgiraldo @mwarusz do you think that the way the 2nd order operator is calculated in NUMA is different than CLIMA? Either that or the calculation of the viscous fluxes. Just ideas here.

[jkozdon]

I'm not quite sure how to go about driving the numerical fluxes for either code.

You'd have to pull the two out of each codes and write a test harness to drive them. This is probably easier with ClimateMachine than NUMA.

There should be some differences between NUMA and CLIMA, I'm pretty sure CLIMA and NUMA are using two different equation sets.

Would be good to run NUMA with the same equation set and same DG formulation as close as we can.

I would suggest running NUMA with the same equation set and run with the weak form (if you are not already).

(Perhaps you could post your NUMA input file?)

However, the paper I've posted here uses the same equation set as CLIMA with DG and does not seem to exhibit this issue Murman_S_AIAA_2014_2784.pdf

This paper also uses a number of other tricks such as over integration. It's also unclear to me whether they are running on conservative or non-conservative form of the equation.

You could run ClimateMachine in a similar state by running with a high-polynomial order and then applying a cutoff filter after every RK stage. In the paper they run 8th order, but are over integrating to 16th order. So you'd run ClimateMachine with N = 16 with an 8th order cutoff filter.

The easiest way to apply the cutoff filter after every stage would be to apply them during the update aux. This will work on a single MPI rank, but with multiple MPI ranks you will need to do it after the update (or before the communication); update aux happen after communication which is why you need to hack things.

@jkozdon @fxgiraldo @mwarusz do you think that the way the 2nd order operator is calculated in NUMA is different than CLIMA? Either that or the calculation of the viscous fluxes. Just ideas here.

Maybe, but seems unlikely to me. But I am often wrong. Someone could try either adding adding symmetric interior penalty to ClimateMachine or adding LDG to NUMA.

I am wondering if it's the form of the equation. Are you using strong form in NUMA? This is a non-conservative form (only conserves mass). Can you try the weak form in NUMA? or strong_no_product?

Also, have you vized the solution at all to see how different they are visually to see if there are any obvious differences?

[slifer50]

Visually both NUMA and CLIMA present very similar flows with no obvious differences. I am not super familiar with NUMA so I'd prefer if @smarras79 or @fxgiraldo answer the question of which equation form is being used in NUMA.
I will try what you suggested in terms of using ClimateMachine with N=16 and an 8th order cutoff filter.

[fxgiraldo]

I need to merge Simone's/Yassine's branch into the main branch and can try running the simulation myself to see if I can see what could account for these differences. Will get to that soon.

[jkozdon]

I am not super familiar with NUMA so I'd prefer if @smarras79 or @fxgiraldo answer the question of which equation form is being used in NUMA.

If it's the vanilla input file for case 9002 it's using a non-conservative formulation with a non-conservative discretization.

I would suggest using
eqn_set='set4c'
and
form_method = 'weak'

I think this should run NUMA in a form similar to what we do in ClimateMachine. Not 100% sure, but from looking at the source I believe this will be should be weak form DG with the same equations. (Big thing I am interested in is whether the product rule tricks NUMA does in the strong form are impacting the results.)

[slifer50]

@jkozdon, I will run NUMA with your suggested input and post the result here

[slifer50]

@fxgiraldo Is there a reference case in NUMA I can use to check that I am setting up the initial conditions correctly for eqn_set='set4c'?

[smarras79]

@fxgiraldo Is there a reference case in NUMA I can use to check that I am setting up the initial conditions correctly for eqn_set='set4c'?

@slifer50 if all is setup properly, you should be able to simply replace whatever eqn_set you have, with eqn_set='set4c' and form_method = 'weak'.
this should work out of the box.
Try and let us know.

[slifer50]

@smarras79, I've tried running it out of the box and it breaks. Pretty sure, the initial conditions need to be changed since the energy equation is different.

[smarras79]

@smarras79, I've tried running it out of the box and it breaks. Pretty sure, the initial conditions need to be changed since the energy equation is different.

I'll let @fxgiraldo confirm, but numa is setup so that the initial condition is always defined in the same way and, based on the eqns value, numa runs through the necessary transformation.

[bischtob]

@slifer50, what is the status of this? Should this still be pinned?