Diffusion with TRT collision operator

Hey everybody!
After having programmed a D2Q9 and D3Q19 model with BGK, TRT collision operators and the corresponding Smagorinsky model with C++ and OpenMP (I validated the code with several benchmark scenarios), I wanted to add a basic advection-diffusion equation that accounts for the passive scalar transport of a mass fraction Y.
Everything worked fine with the BGK operator but I ran into issues with the two relaxation time TRT collision operator. Weirdly enough my simulated physics change as I adjust the second relaxation time. If I enforce omega- = omega+ the result is correct but as soon as I deviate from this, for example by enforcing the magical parameter, the model behaves in a way I can’t yet understand: If I set omega- higher than omega+ the diffusion happens slower and if I set it lower the diffusion happens faster than it should. The contribution of the anti-symmetric part becomes comparably large but I can’t see my mistake.

I then took out the advection part, now the code is basically diffusion only but still exhibits the same behaviour. Here the relevant code segment:

const double tau_d = (6.0*diffusivity + 1.0) / 2.0;
const double magic_param = 1.0 / 6.0;

const double omegap = 1/tau_d;
const double omegam = (1.0/omegap – 1.0/2.0) / (magic_param + 1.0/2.0*( 1.0/omegap – 1.0/2.0));

for(unsigned int y = 0; y < NY; ++y)
{
for(unsigned int x = 0; x < NX; ++x)
{
/// temporary populations (load populations from main memory)
double gt0 = g0[D2Q9_ScalarIndex(x,y)];
double gt1 = g1[D2Q9_FieldIndex(x,y,1)];
double gt2 = g1[D2Q9_FieldIndex(x,y,2)];
double gt3 = g1[D2Q9_FieldIndex(x,y,3)];
double gt4 = g1[D2Q9_FieldIndex(x,y,4)];
double gt5 = g1[D2Q9_FieldIndex(x,y,5)];
double gt6 = g1[D2Q9_FieldIndex(x,y,6)];
double gt7 = g1[D2Q9_FieldIndex(x,y,7)];
double gt8 = g1[D2Q9_FieldIndex(x,y,8)];

/// microscopic to macroscopic: mass fractions
double Ymt = gt0 + gt1 + gt2 + gt3 + gt4 + gt5 + gt6 + gt7 + gt8;
Ym[D2Q9_ScalarIndex(x,y)] = Ymt;

/// equilibrium distribution: only diffusion (no velocity!)
double geq1 = Ymt*wc;
double geq2 = Ymt*wc;
double geq3 = Ymt*wc;
double geq4 = Ymt*wc;
double geq5 = Ymt*wd;
double geq6 = Ymt*wd;
double geq7 = Ymt*wd;
double geq8 = Ymt*wd;

/// symmetric and anti-symmetric populations
double g13p = (gt1 + gt3 – (geq1 + geq3))/2;
double g13m = (gt1 – gt3 – (geq1 – geq3))/2;
double g24p = (gt2 + gt4 – (geq2 + geq4))/2;
double g24m = (gt2 – gt4 – (geq2 – geq4))/2;
double g57p = (gt5 + gt7 – (geq5 + geq7))/2;
double g57m = (gt5 – gt7 – (geq5 – geq7))/2;
double g68p = (gt6 + gt8 – (geq6 + geq8))/2;
double g68m = (gt6 – gt8 – (geq6 – geq8))/2;

/// neighbouring cells coordinates (assume periodic boundaries
unsigned int xp = (x + 1) % NX;
unsigned int yp = (y + 1) % NY;
unsigned int xm = (NX + x – 1) % NX;
unsigned int ym = (NY + y – 1) % NY;

/// collision & streaming
g0[D2Q9_ScalarIndex(x,y)] = gt0 + omegap*(Ymt*w0 – gt0);
g2[D2Q9_FieldIndex(xp, y, 1)] = gt1 – omegap*g13p – omegam*g13m;
g2[D2Q9_FieldIndex(xm, y, 3)] = gt3 – omegap*g13p + omegam*g13m;
g2[D2Q9_FieldIndex(x, yp, 2)] = gt2 – omegap*g24p – omegam*g24m;
g2[D2Q9_FieldIndex(x, ym, 4)] = gt4 – omegap*g24p + omegam*g24m;
g2[D2Q9_FieldIndex(xp, yp, 5)] = gt5 – omegap*g57p – omegam*g57m;
g2[D2Q9_FieldIndex(xm, ym, 7)] = gt7 – omegap*g57p + omegam*g57m;
g2[D2Q9_FieldIndex(xm, yp, 6)] = gt6 – omegap*g68p – omegam*g68m;
g2[D2Q9_FieldIndex(xp, ym, 8)] = gt8 – omegap*g68p + omegam*g68m;
}
}

I am thankful for every hint!
Regards