From 608297d1aa8835bf64a25f46cfbd6f6f4f2124df Mon Sep 17 00:00:00 2001 From: Alessandro De Rosis Date: Mon, 7 Mar 2022 12:35:59 +0000 Subject: [PATCH] written central moments in a smarter way --- .gitignore | 1 + src/D2Q9_CentralMoments_LIFE.m | 128 +++++++++++++++++++++++++++++++++ src/Grid.cpp | 117 ++++++++---------------------- 3 files changed, 159 insertions(+), 87 deletions(-) create mode 100644 src/D2Q9_CentralMoments_LIFE.m diff --git a/.gitignore b/.gitignore index dc7159b..22798a5 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ /examples/**/Results/ /obj/ /testing/**/ +.DS_Store diff --git a/src/D2Q9_CentralMoments_LIFE.m b/src/D2Q9_CentralMoments_LIFE.m new file mode 100644 index 0000000..e45a1f2 --- /dev/null +++ b/src/D2Q9_CentralMoments_LIFE.m @@ -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) diff --git a/src/Grid.cpp b/src/Grid.cpp index 25ba36f..a8b41f8 100644 --- a/src/Grid.cpp +++ b/src/Grid.cpp @@ -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]; @@ -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]; @@ -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 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++) {