Bug in LHS evaluation of tensors with no bases (for a ball problem, maybe others)

[evanhanders]

Hi folks,

I've been working on getting the fully compressible equations running and conserving energy. In those equation, there is a viscous term like div(e - (1/3) * div_u * I), where I is the identity matrix. I've found that I get incorrect results when I define the identity matrix like so:

eye = dist.TensorField(coords, name='eye')

But I get correct results when I define it like this:

eye = dist.TensorField(coords, name='eye', bases=basis.radial_basis)

If I evaluate operations involving these two definitions of eye, then I find that they give the same result when they're included on the RHS of equations but not when they're included on the LHS, which is where I want them for my diffusive term!

I've thrown together a MWE by modifying the ball internally heated example; see here: Imatrix_internally_heated_convection.py.txt. I've done the following:

1. Set n_phi = 1 for azimuthal symmetry (runs faster). Also lowered Ra and nr, ntheta so that it runs fast on one core.
2. Defined eye using both of the definitions above. I'm calling the first definition "smol" and the second "ncc".
3. Added test fields with equations of the form "dt(test_field) = div(eye*T), which is kinda like the term in the fully compressible equations.
4. measured the value of the test field for different formulations (smol vs ncc, RHS vs LHS evaluation).

I'm taking the ncc case where the div() term is on the LHS to be the "truth" and comparing the other cases to it. I get something like this (which the attached script will produce if you run it):

...so we can see that the smol identity matrix when evaluated on the LHS gives the wrong answer. A simple fix for now seems to be to use the ncc version, but I wanted to make this known!

[kburns]

Thanks for the report. Here's a minimum working example showing the failure:
implicit_identity_mwe.py.txt

[evanhanders]

Thanks, @kburns ! Here's an alternate MWE in test form inside of test_spherical_ncc.py that Daniel and I came up with: test_spherical_ncc.py.txt

(see "test_identity_multiply")

[kburns]

Looks like the issue is that the intertwiners are not applied to the components of constant fields (a choice we actively made so that spherical BC fields were in spin space instead of regularity space). Not sure there's a really quick fix for this unfortunately.

[kburns]

Maybe the best we can do right now is raise an error and force you to use NCCs with bases in curvilinear problems.

[evanhanders]

As a temporary fix, I think raising an error and saying to make it an NCC works! That's a trivial difference to make to my script as a user, especially if it makes the solution correct.