www.pudn.com > collide_mrt.zip > collide_mrt.m, change:2017-04-12,size:2157b

```function f = collide_mrt(f, u, v, rho, om)
% D2Q9 collisions on 2-d matrix.
% Multiple relaxation time formulation.
% om is omega, the relaxation frequency, or 1/tau (inverse of relaxation time).

% Constants
M = [ones(1,9);...
-4, -ones(1,4), 2*ones(1,4);
4, -2*ones(1,4), ones(1,4);
0, 1, 0, -1, 0, 1, -1, -1, 1;
0, -2, 0, 2, 0, 1, -1, -1, 1;
0, 0, 1, 0, -1, 1, 1, -1, -1;
0, 0, -2, 0, 2, 1, 1, -1, -1;
0, 1, -1, 1, -1, zeros(1,4);
zeros(1,5), 1, -1, 1 -1];
Minv = 1/36*[4*ones(9,1), M(2,:)', M(3,:)', 6*M(4,:)', 3*M(5,:)',...
6*M(6,:)', 3*M(7,:)', 9*M(8,:)', 9*M(9,:)'];
% S_vec = [1, 1.4, 1.4, 1, 1.2, 1, 1.2, om, om]'; % Mohamad.
S_vec = [1, 1.2, 1, 1, 1.2, 1, 1.2, om, om]'; % 2015 Zhang et al.
MinvS = Minv*diag(S_vec);
[rows, cols] = size(rho);

% Different orientation of vector indices
reorder = [1,2,6,3,7,4,8,5,9];
M_ = M(:,reorder);
Minv_coeff = 36*inv(M_);
% Minv_ = Minv_coeff / 36;
% S_vec_ = S_vec(reorder);
% MinvS_ = Minv_*diag(S_vec_);
% MinvS_ = Minv_*diag(S_vec);
% MinvS_coeff = 36*(Minv_*diag(S_vec));
S_vec2 = S_vec;
S_vec2([8,9]) = 1;
mdiff = Minv*diag(S_vec2) - MinvS;
mdiff_coeff = mdiff*36;

% Determine m - meq.
deltam = zeros(rows, cols, 9);
% First fill with m.
for k = 1:9
for n = 1:9
deltam(:,:,k) =  deltam(:,:,k) + M(k,n)*f(:,:,n);
end
end
% Now let us subtract meq.
ru = rho.*u;
rv = rho.*v;
u2 = u.^2;
v2 = v.^2;
uu = u2+v2;
deltam(:,:,1) = deltam(:,:,1) - rho;
% deltam(:,:,2) = deltam(:,:,2) - rho.*( -2 + 3*rho.*uu ); % Mohamad.
% deltam(:,:,3) = deltam(:,:,3) - rho.*( 1 - 3*rho.*uu ); % Mohamad.
deltam(:,:,2) = deltam(:,:,2) - rho.*(-2 + 3*uu); % 2015 Zhang et al.
deltam(:,:,3) = deltam(:,:,3) - rho.*(1 - 3*uu); % 2015 Zhang et al.
deltam(:,:,4) = deltam(:,:,4) - ru;
deltam(:,:,5) = deltam(:,:,5) + ru;
deltam(:,:,6) = deltam(:,:,6) - rv;
deltam(:,:,7) = deltam(:,:,7) + rv;
deltam(:,:,8) = deltam(:,:,8) - rho.*(u2-v2);
deltam(:,:,9) = deltam(:,:,9) - rho.*u.*v;
% Update f.
for k = 1:9
for n = 1:9
f(:,:,k) =  f(:,:,k) ...
- MinvS(k,n) * deltam(:,:,n);
end
end```