Results should not be sensitive to choice of budget decomposition

[BrodiePearson]

I think the results of the ADC closure are quite sensitive to how we decompose the tendency equations! (possibly due to rounding errors?)

I was playing around with the w2 budget, and found that simply combining two tendency terms led to significant differences in final simulated profiles. In this case I ran two simulations where I either left the return to isotropy and buoyancy terms separate:

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines 571 to 573 in 103e68e

	w2tend3(k,iCell) = tauvVel(k,iCell)*(u2(i1,k,iCell) + v2(i1,k,iCell))/3.0_RKIND
	w2tend4(k,iCell) = (2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_b)*Mc(k,iCell)* &
	(grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell))

or I combined them into one term (which theoretically should be an identical system of equations):

w2tend3(k,iCell) = tauvVel(k,iCell)*(u2(i1,k,iCell) + v2(i1,k,iCell))/3.0_RKIND + &
(2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_b)*Mc(k,iCell)* &
(grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell))
w2tend4(k,iCell) = 0.0_RKIND*(2.0_RKIND - 4.0_RKIND/3.0_RKIND*C_b)*Mc(k,iCell)* &
(grav*alphaT(k,iCell)*tumd(k,iCell) - grav*betaS(k,iCell)*sumd(k,iCell))

This figure shows the final profiles of w2 (left) after 24 hours of the conv2 setup using either split (standard) terms or the combined tendency terms. It also shows the sum w2tend3 + w2tend4 (middle) and the individual terms (right).

This apparent sensitivity of w2 profiles to equation formulation is not good, so it would be great if someone could do a similar test and see if these differences are real/reproducible.

The maximum difference of total tendency terms (middle panel) between the two simulations grows from 1e-23 on the first time step to 4e-8 after 24 hours. Perhaps this is a truncation error that grows - although it is surprising that it produces such a different w2 profile given that the model doesn't have the chaos of a truly turbulent flow.

[vanroekel]

@BrodiePearson this is very odd indeed. can you post your namelist and the configuration parameters for the init step? It seems highly unlikely this is truncation error. If you can share the setup I'll try run this.

[BrodiePearson]

@vanroekel I've attached the folder that I'm using. It should be identical to the one that you originally sent me, with modifications to work for the new set of constants. I hadn't created any other examples yet, so I would be interested to know if these differences are (a) an error on my part, or (b) only specific to this case.
cvmix_test_conv2_new_constants.zip

[katsmith133]

@BrodiePearson and @vanroekel
For me, combining those two terms does not seem to cause the same issue, even after 82 hours...

Note: The difference in the w2tend1 + w2tend_pressure term is to be expected as it involves w2ten3 and t2tend4.

...at least not for the default namelist parameters I am using. I am now going to test using the namelist parameters Brodie is using and let you know what I get...

[katsmith133]

Using Brodie's namelist, I get the following after 24 hours:

It does appear that combining the two terms does change the w2 and w3 profiles slightly with Brodie's namelist parameter values (note: I used my init and forcing files, so it is not expected that I would get the same differences as Brodie).

In case it helps up determine why this change is happening, here are the differences in the namelist values between the previous comment's run (where no change is seen) and this comment's run with Brodie's parameters (where changes are seen):
parameter = Brodie's namelist param values -> Kat's namelist param values
config_dt = 0:01:00 -> 1:00:00
config_adc_Cmom = 0.0 -> 4.0
config_adc_Ctherm = 0.0 -> 4.0
config_adc_alpha_0 = 10.8 -> 0.8
config_adc_bc_const = 1.8 -> 1.0
config_adc_bc_const_wp2 = 1.8 -> 0.0
config_adc_up2_vp2_factor = 4.0 -> 1.0
config_adc_use_splat_parameterization = .true. -> .false.

[vanroekel]

Thanks @katsmith133 a few questions

1. are you using the same code for both runs?
2. Do you have any ideas what parameter is the key here? none of them really make sense to me as tripping this issue. The most likely to me would be Cmom, Ctherm or config_dt. I'll try test as well.

[katsmith133]

@vanroekel

1. I assume we are in as far as I am using the most recent version of the code on Github and made only that small same change as Brodie showed before. But I am not using the same executable, no. I can try that.
2. Agreed. I do not know which one it is, but would assume at first it is one of those you mentioned. I can change each one in turn and see which one(s) it is.

[vanroekel]

In my testing I see the same as @BrodiePearson even with config_dt = 00:00:01 so that is not the issue. If I reenable Cmom and Ctherm the diffs are much smaller (1E-17), but the model crashes 6 hours in. It makes no sense to me why disabling Cmom/Ctherm exposes this issue. I'll try look through the code some more.

[BrodiePearson]

@vanroekel @katsmith133 Thanks for both taking a look. I followed Kat's lead and found that her setup gave me very small differences in the final w2 profiles (~1e-21; we might expect this from small truncation differences between combining the terms or calculating them separately right?). It looks to me like config_adc_bc_const_wp2 is probably the main issue (amplifies this error?). I can run with all my parameters that Kat listed above, except with config_adc_bc_const_wp2 = 0.0 and the difference is larger (1e-9) but arguably negligible (<0.1%). But switching to config_adc_bc_const_wp2 = 1.8 increases the difference to ~1e-5 or ~100%. Currently stumped why that is...

On Lukes point, certain combinations of these parameters seem to cause crashes (I think Kat ran into this issue previously), I tried different combinations and found that I could disable Cmom/Ctherm as long as I switched splatting on. In this case, with all other parameters using Kats settings, the differences are ~1e-9

[vanroekel]

Hey all, just an update. If I do config_adc_bc_const_wp2 = 0.0 and config_use_splat_parameterization = .false.. the two options listed here are BFB (diff is hard zero). This is starting to have the feel of a memory index issue to me. I'll keep looking.

[BrodiePearson]

Hmmm. I just looked at the w2tend6 (splat) term and noticed there is a typo in the allocation code:

MPAS-Model/src/core_ocean/shared/mpas_ocn_turbulence.F

Line 380 in 103e68e

w3tend6 = w2tend6Tmp

and it is not deallocated:

MPAS-Model/src/core_ocean/shared/mpas_ocn_vmix_adc.F

Line 532 in 103e68e

w2tend5, wttend1, wttend2, wttend3, &

Could that be the problem?

[vanroekel]

Hah, just found that too! I'm testing if this is the problem. The deallocate won't matter, but the other might. I'm a bit doubtful though given it only influences splat, but I'll let you know.

[vanroekel]

unfortunately that wasn't it. Did not change my results. Important to fix though. Thanks for pointing that out.