Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

written central moments in a smarter way #4

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

written central moments in a smarter way #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,3 +4,4 @@
/examples/**/Results/
/obj/
/testing/**/
.DS_Store
128 changes: 128 additions & 0 deletions src/D2Q9_CentralMoments_LIFE.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
clear all
clc

%% Initialize some symbolic variables
syms U V R omega f0 f1 f2 f3 f4 f5 f6 f7 f8 Fx Fy real
%% Define lattice directions, weights and other useful quantities of the D2Q9 model
cx = [0 1 -1 0 0 1 -1 1 -1];
cy = [0 0 0 1 -1 1 -1 -1 1];
w = [4./9., 1./9., 1./9., 1./9., 1./9., 1./36., 1./36., 1./36., 1./36.];
cs = 1/sqrt(3);
cs2 = cs^2;
cs3 = cs^3;
cs4 = cs^4;
cs6 = cs^6;
cs8 = cs^8;

f = [f0 f1 f2 f3 f4 f5 f6 f7 f8]'; % hydrodynamic populations
feq = sym(zeros(9,1));

Force = sym(zeros(9,1)); % generic force vector
T = sym(zeros(9,9)); %transformation matrix for central moments
M = zeros(9,9); %transformation matrix for raw moments
Lambda = diag([1, 1, 1, 1, omega, omega, 1, 1, 1]);
Id = eye(9,9); % identity matrix
for i=1:9
% build the complete equilibria
first_order = 1/cs2*(U*cx(i)+V*cy(i));
second_order = 1/(2*cs4)*((cx(i)*cx(i)-1/3)*U^2+...
(cy(i)*cy(i)-1/3)*V^2+...
2*cx(i)*cy(i)*U*V);
third_order = 1/(2*cs6)*((cx(i)^2-1/3)*cy(i)*U*U*V+(cy(i)^2-1/3)*cx(i)*U*V*V);
fourth_order = 1/(4*cs8)*((cx(i)^2-1/3)*(cy(i)^2-1/3)*U*U*V*V);
feq(i) = w(i)*R*(1+first_order+second_order+third_order+fourth_order);

% build the complete forcing terms
hat_cx = cx(i)/cs;
hat_cy = cy(i)/cs;
hat_ = [hat_cx, hat_cy];
for a=1:2
H1(a) = hat_(a);
for b=1:2
H2(a,b) = hat_(a)*hat_(b)-Id(a,b);
for c=1:2
hat_I = hat_(a)*Id(b,c)+hat_(b)*Id(a,c)+hat_(c)*Id(a,b);
H3(a,b,c) = hat_(a)*hat_(b)*hat_(c)-hat_I;
for d=1:2
hat_II = hat_(a)*hat_(b)*Id(c,d)+hat_(a)*hat_(c)*Id(b,d)+...
hat_(a)*hat_(d)*Id(b,c)+hat_(b)*hat_(c)*Id(a,d)+...
hat_(b)*hat_(d)*Id(a,c)+hat_(c)*hat_(d)*Id(a,b);
II = Id(a,b)*Id(c,d)+Id(a,c)*Id(b,d)+Id(a,d)*Id(b,c);
H4(a,b,c,d) = hat_(a)*hat_(b)*hat_(c)*hat_(d)-hat_II+II;
end
end
end
end
first_order = 1/cs*(Fx*H1(1)+Fy*H1(2));
second_order = 1/(2*cs2)*( (Fx*U+U*Fx)*H2(1,1) +...
(Fy*V+V*Fy)*H2(2,2) +...
(Fx*V+Fy*U)*H2(1,2) +...
(Fy*U+Fx*V)*H2(2,1) );
third_order = 1/(6*cs3)*( (Fx*V*V+U*Fy*V+U*V*Fy)*H3(1,2,2) +...
(Fy*U*V+V*Fx*V+V*U*Fy)*H3(2,1,2) +...
(Fy*V*U+V*Fy*U+V*V*Fx)*H3(2,2,1) +...
(Fx*U*V+U*Fx*V+U*U*Fy)*H3(1,1,2) +...
(Fy*U*U+V*Fx*U+V*U*Fx)*H3(2,1,1) +...
(Fx*V*U+U*Fy*U+U*V*Fx)*H3(1,2,1) );
fourth_order = 1/(24*cs4)*( (Fx*U*V*V+U*Fx*V*V+U*U*Fy*V+U*U*V*Fy)*H4(1,1,2,2)+...
(Fx*V*U*V+U*Fy*U*V+U*V*Fx*V+U*V*U*Fy)*H4(1,2,1,2)+...
(Fx*V*V*U+U*Fy*V*U+U*V*Fy*U+U*V*V*Fx)*H4(1,2,2,1)+...
(Fy*U*U*V+V*Fx*U*V+V*U*Fx*V+V*U*U*Fy)*H4(2,1,1,2)+...
(Fy*U*V*U+V*Fx*V*U+V*U*Fy*U+V*U*V*Fx)*H4(2,1,2,1)+...
(Fy*V*U*U+V*Fy*U*U+V*V*Fx*U+V*V*U*Fx)*H4(2,2,1,1) );
Force(i) = w(i)*(first_order + second_order + third_order + fourth_order)/R;

% build the transformation matrix T
CX = cx(i)-U;
CY = cy(i)-V;
T(1,i) = 1;
T(2,i) = CX;
T(3,i) = CY;
T(4,i) = CX*CX+CY*CY;
T(5,i) = CX*CX-CY*CY;
T(6,i) = CX*CY;
T(7,i) = CX*CX*CY;
T(8,i) = CX*CY*CY;
T(9,i) = CX*CX*CY*CY;

% build the tranformation matrix M
CX = cx(i);
CY = cy(i);
M(1,i) = 1;
M(2,i) = CX;
M(3,i) = CY;
M(4,i) = CX*CX+CY*CY;
M(5,i) = CX*CX-CY*CY;
M(6,i) = CX*CY;
M(7,i) = CX*CX*CY;
M(8,i) = CX*CY*CY;
M(9,i) = CX*CX*CY*CY;
end
T = simplify(T);
N = simplify(T*M^(-1)); %shift matrix
%N = Id;
%T = M;
%% HYDRO
syms k0_pre k1_pre k2_pre k3_pre k4_pre k5_pre k6_pre k7_pre k8_pre real
syms k1_star k2_star k3_star k4_star k5_star k6_star k7_star k8_star real
K_pre = simplify(T*f); %pre-collision central moments
K_eq = simplify(T*feq); % equilibrium central moments
K_force = simplify(T*Force); % forcing term central moments
K_pre(5) = k4_pre;
K_pre(6) = k5_pre;
K_star = simplify((Id-Lambda)*K_pre + Lambda*K_eq + (Id-Lambda/2)*K_force ) %post-collision central moments

%post-collision populations
K_sym = [R k1_star k2_star k3_star k4_star k5_star k6_star k7_star k8_star];
for i=1:9
if(K_star(i)~=sym(0))
K_star(i) = K_sym(i);
end
end
f_post_collision_onestep = collect(simplify(T \ K_star), K_star);

% two-steps approach
raw_moments = simplify(N\K_star)
syms r0 r1 r2 r3 r4 r5 r6 r7 r8 real
r = [r0 r1 r2 r3 r4 r5 r6 r7 r8]'; %symbolic raw moments
f_post_collision_twosteps = collect(simplify(M\r),K_star)
117 changes: 30 additions & 87 deletions src/Grid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ inline void GridClass::streamCollide(int i, int j, int id) {
double k5Pre = 0.0;

// Loop through velocities on src_id
for (int v = 0; v < nVels; v++) {
/*for (int v = 0; v < nVels; v++) {

// Shift lattice velocities
double cx = c[v * dims + eX] - u_n[id * dims + eX];
Expand All @@ -119,7 +119,16 @@ inline void GridClass::streamCollide(int i, int j, int id) {
// Transform f to central moments (k)
k4Pre += f_n[id * nVels + v] * (SQ(cx) - SQ(cy));
k5Pre += f_n[id * nVels + v] * cx * cy;
}
}*/
// Get velocity
double ux = u_n[id * dims + eX];
double uy = u_n[id * dims + eY];
double ux2 = ux*ux;
double uy2 = uy*uy;

// Required pre-collision central moments
double k4Pre = f_n[id * nVels + 1] + f_n[id * nVels + 2] - f_n[id * nVels + 3] - f_n[id * nVels + 4] + rho_n[id]*(uy*uy-ux*ux);
double k5Pre = f_n[id * nVels + 5] + f_n[id * nVels + 6] - f_n[id * nVels + 7] - f_n[id * nVels + 8] + rho_n[id]*ux*uy;

// Get post-collision central moments
double k0 = rho_n[id];
Expand All @@ -132,95 +141,29 @@ inline void GridClass::streamCollide(int i, int j, int id) {
double k7 = 0.5 * (force_xy[id * dims + eX] + force_ibm[id * dims + eX]) * SQ(c_s);
double k8 = rho_n[id] * QU(c_s);

// Get velocity
double ux = u_n[id * dims + eX];
double uy = u_n[id * dims + eY];
double r0 = k0;
double r1 = k1 + rho_n[id]*ux;
double r2 = k2 + rho_n[id]*uy;
double r3 = k3 + 2.0*ux*k1 + 2.0*uy*k2 + rho_n[id]*(ux2 + uy2);
double r4 = k4 + 2.0*ux*k1 - 2.0*uy*k2 + rho_n[id]*(ux2 - uy2);
double r5 = k5 + ux*k2 + uy*k1 + rho_n[id]*ux*uy;
double r6 = k6 + 2.0*ux*k5 + 0.5*uy*(k3+k4) + ux2*k2 + 2.0*ux*uy*k1 + rho_n[id]*ux2*uy;
double r7 = k7 + 0.5*ux*(k3-k4) + 2.0*uy*k5 + uy2*k1 + 2.0*ux*uy*k2 + rho_n[id]*ux*uy2;
double r8 = k8 + 2.0*ux*k7 + 2.0*uy*k6 + 0.5*k3*(ux2 + uy2) - 0.5*k4*(ux2 - uy2) + rho_n[id]*ux2*uy2 + 4.0*ux*uy*k5 + 2.0*ux*uy2*k1 + 2.0*ux2*uy*k2;

// Declare array
array<double, nVels> fStar;

// Set array (transform from central moments to f)
fStar[0] = (SQ(ux) * SQ(uy) - SQ(ux) - SQ(uy) + 1.0) * k0
+ (2.0 * ux * SQ(uy) - 2.0 * ux) * k1
+ (2.0 * uy * SQ(ux) - 2.0 * uy) * k2
+ (0.5 * SQ(ux) + 0.5 * SQ(uy) - 1.0) * k3
+ (0.5 * SQ(uy) - 0.5 * SQ(ux)) * k4
+ 4.0 * ux * uy * k5
+ 2.0 * uy * k6
+ 2.0 * ux * k7
+ k8;
fStar[1] = (0.5 * SQ(ux) - 0.5 * (SQ(ux) * SQ(uy)) - 0.5 * (ux * SQ(uy)) + 0.5 * ux) * k0
+ (ux - ux * SQ(uy) - 0.5 * SQ(uy) + 0.5) * k1
+ (-uy * SQ(ux) - uy * ux) * k2
+ (-0.25 * SQ(ux) - 0.25 * ux - 0.25 * SQ(uy) + 0.25) * k3
+ (0.25 * SQ(ux) + 0.25 * ux - 0.25 * SQ(uy) + 0.25) * k4
+ (-uy - 2.0 * ux * uy) * k5
+ (-uy) * k6
+ (-ux - 0.5) * k7
- 0.5 * k8;
fStar[2] = (-0.5 * (SQ(ux) * SQ(uy)) + 0.5 * SQ(ux) + 0.5 * (ux * SQ(uy)) - 0.5 * ux) * k0
+ (ux - ux * SQ(uy) + 0.5 * SQ(uy) - 0.5) * k1
+ (-uy * SQ(ux) + uy * ux) * k2
+ (-0.25 * SQ(ux) + 0.25 * ux - 0.25 * SQ(uy) + 0.25) * k3
+ (0.25 * SQ(ux) - 0.25 * ux - 0.25 * SQ(uy) + 0.25) * k4
+ (uy - 2.0 * ux * uy) * k5
+ (-uy) * k6
+ (0.5 - ux) * k7
- 0.5 * k8;
fStar[3] = (0.5 * SQ(uy) - 0.5 * (SQ(ux) * uy) - 0.5 * (SQ(ux) * SQ(uy)) + 0.5 * uy) * k0
+ (-ux * SQ(uy) - ux * uy) * k1
+ (uy - SQ(ux) * uy - 0.5 * SQ(ux) + 0.5) * k2
+ (-0.25 * SQ(ux) - 0.25 * SQ(uy) - 0.25 * uy + 0.25) * k3
+ (0.25 * SQ(ux) - 0.25 * SQ(uy) - 0.25 * uy - 0.25) * k4
+ (-ux - 2.0 * ux * uy) * k5
+ (-uy - 0.5) * k6
+ (-ux) * k7
- 0.5 * k8;
fStar[4] = (-0.5 * (SQ(ux) * SQ(uy)) + 0.5 * (SQ(ux) * uy) + 0.5 * SQ(uy) - 0.5 * uy) * k0
+ (-ux * SQ(uy) + ux * uy) * k1
+ (uy - SQ(ux) * uy + 0.5 * SQ(ux) - 0.5) * k2
+ (-0.25 * SQ(ux) - 0.25 * SQ(uy) + 0.25 * uy + 0.25) * k3
+ (0.25 * SQ(ux) - 0.25 * SQ(uy) + 0.25 * uy - 0.25) * k4
+ (ux - 2.0 * ux * uy) * k5
+ (0.5 - uy) * k6
+ (-ux) * k7
- 0.5 * k8;
fStar[5] = (0.25 * (SQ(ux) * SQ(uy)) + 0.25 * (SQ(ux) * uy) + 0.25 * (ux * SQ(uy)) + 0.25 * (ux * uy)) * k0
+ (0.25 * uy + 0.5 * (ux * uy) + 0.5 * (ux * SQ(uy)) + 0.25 * SQ(uy)) * k1
+ (0.25 * ux + 0.5 * (ux * uy) + 0.5 * (SQ(ux) * uy) + 0.25 * SQ(ux)) * k2
+ (0.125 * SQ(ux) + 0.125 * ux + 0.125 * SQ(uy) + 0.125 * uy) * k3
+ (-0.125 * SQ(ux) - 0.125 * ux + 0.125 * SQ(uy) + 0.125 * uy) * k4
+ (0.5 * ux + 0.5 * uy + ux * uy + 0.25) * k5
+ (0.5 * uy + 0.25) * k6
+ (0.5 * ux + 0.25) * k7
+ 0.25 * k8;
fStar[6] = (0.25 * (SQ(ux) * SQ(uy)) - 0.25 * (SQ(ux) * uy) - 0.25 * (ux * SQ(uy)) + 0.25 * (ux * uy)) * k0
+ (0.25 * uy - 0.5 * (ux * uy) + 0.5 * (ux * SQ(uy)) - 0.25 * SQ(uy)) * k1
+ (0.25 * ux - 0.5 * (ux * uy) + 0.5 * (SQ(ux) * uy) - 0.25 * SQ(ux)) * k2
+ (0.125 * SQ(ux) - 0.125 * ux + 0.125 * SQ(uy) - 0.125 * uy) * k3
+ (-0.125 * SQ(ux) + 0.125 * ux + 0.125 * SQ(uy) - 0.125 * uy) * k4
+ (ux * uy - 0.5 * uy - 0.5 * ux + 0.25) * k5
+ (0.5 * uy - 0.25) * k6
+ (0.5 * ux - 0.25) * k7
+ 0.25 * k8;
fStar[7] = (0.25 * (SQ(ux) * SQ(uy)) - 0.25 * (SQ(ux) * uy) + 0.25 * (ux * SQ(uy)) - 0.25 * (ux * uy)) * k0
+ (0.5 * (ux * SQ(uy)) - 0.5 * (ux * uy) - 0.25 * uy + 0.25 * SQ(uy)) * k1
+ (0.5 * (ux * uy) - 0.25 * ux + 0.5 * (SQ(ux) * uy) - 0.25 * SQ(ux)) * k2
+ (0.125 * SQ(ux) + 0.125 * ux + 0.125 * SQ(uy) - 0.125 * uy) * k3
+ (-0.125 * SQ(ux) - 0.125 * ux + 0.125 * SQ(uy) - 0.125 * uy) * k4
+ (0.5 * uy - 0.5 * ux + ux * uy - 0.25) * k5
+ (0.5 * uy - 0.25) * k6
+ (0.5 * ux + 0.25) * k7
+ 0.25 * k8;
fStar[8] = (0.25 * (SQ(ux) * SQ(uy)) + 0.25 * (SQ(ux) * uy) - 0.25 * (ux * SQ(uy)) - 0.25 * (ux * uy)) * k0
+ (0.5 * (ux * uy) - 0.25 * uy + 0.5 * (ux * SQ(uy)) - 0.25 * SQ(uy)) * k1
+ (0.5 * (SQ(ux) * uy) - 0.5 * (ux * uy) - 0.25 * ux + 0.25 * SQ(ux)) * k2
+ (0.125 * SQ(ux) - 0.125 * ux + 0.125 * SQ(uy) + 0.125 * uy) * k3
+ (-0.125 * SQ(ux) + 0.125 * ux + 0.125 * SQ(uy) + 0.125 * uy) * k4
+ (0.5 * ux - 0.5 * uy + ux * uy - 0.25) * k5
+ (0.5 * uy + 0.25) * k6
+ (0.5 * ux - 0.25) * k7
+ 0.25 * k8;
// Set array (transform from raw moments to f)
fStar[0] = r0 - r3 + r8;
fStar[1] = 0.5*(r1-r7-r8) + 0.25*(r3 + r4);
fStar[2] = 0.5*(-r1+r7-r8) + 0.25*(r3 + r4);
fStar[3] = 0.5*(r2-r6-r8) + 0.25*(r3 - r4);
fStar[4] = 0.5*(-r2+r6-r8) + 0.25*(r3 - r4);
fStar[5] = 0.25*(r5 + r6 + r7 + r8);
fStar[6] = 0.25*(r5 - r6 - r7 + r8);
fStar[7] = 0.25*(-r5 - r6 + r7 + r8);
fStar[8] = 0.25*(-r5 + r6 - r7 + r8);

// Now set to f
for (int v = 0; v < nVels; v++) {
Expand Down