#pragma once

#include "concepts.h"

struct SmagorinskyBgkCollideO {

template <concepts::Arithmetic V, concepts::Arithmetic T>
static void apply(V f_curr[19], V f_next[19], T tau, T smagorinsky) {
	V rho;
	V u[3];
	V x0 = f_curr[10] + f_curr[13] + f_curr[17];
	V x1 = f_curr[14] + f_curr[5] + f_curr[6];
	V x2 = f_curr[1] + f_curr[2] + f_curr[4];
	V x3 = x0 + x1 + x2 + f_curr[0] + f_curr[11] + f_curr[12] + f_curr[15] + f_curr[16] + f_curr[18] + f_curr[3] + f_curr[7] + f_curr[8] + f_curr[9];
	V x4 = -f_curr[11] + f_curr[7];
	V x5 = -f_curr[15] + f_curr[3];
	V x6 = T{1}/x3;
	V x7 = f_curr[0] - f_curr[18];
	V x8 = T{2.00000000000000}*tau + T{-1.00000000000000};
	V x9 = T{72.0000000000000}*f_curr[5];
	V x26 = T{72.0000000000000}*f_curr[13];
	V x33 = T{72.0000000000000}*f_curr[11];
	V x44 = T{72.0000000000000}*f_curr[1];
	V x50 = T{72.0000000000000}*f_curr[17];
	V x56 = T{72.0000000000000}*f_curr[15];
	V x64 = T{72.0000000000000}*f_curr[4];
	V x68 = T{72.0000000000000}*f_curr[14];
	V x72 = T{72.0000000000000}*f_curr[18];
	rho = x3;
	V x78 = T{2.00000000000000}*rho;
	u[0] = x6*(x0 + x4 + x5 - f_curr[1] - f_curr[5] - f_curr[8]);
	V x12 = -u[0];
	V x14 = T{6.00000000000000}*u[0];
	V x17 = u[0]*u[0];
	V x18 = T{3.00000000000000}*x17;
	V x27 = -x14;
	V x38 = T{2.00000000000000} - x18;
	V x79 = T{6.00000000000000}*x17;
	u[1] = x6*(x1 + x4 + x7 - f_curr[12] - f_curr[13] - f_curr[4]);
	V x10 = T{6.00000000000000}*u[1];
	V x11 = -x10;
	V x13 = x12 + u[1];
	V x19 = u[1]*u[1];
	V x20 = T{3.00000000000000}*x19;
	V x21 = x18 + x20 + T{-2.00000000000000};
	V x28 = -u[1];
	V x29 = x28 + u[0];
	V x34 = u[0] + u[1];
	V x35 = T{9.00000000000000}*(x34*x34);
	V x37 = -x20;
	V x39 = x37 + x38;
	V x85 = T{6.00000000000000}*x19;
	u[2] = x6*(x2 + x5 + x7 - f_curr[14] - f_curr[16] - f_curr[17]);
	V x15 = u[2]*u[2];
	V x16 = T{3.00000000000000}*x15;
	V x22 = x16 + x21;
	V x23 = x14 + x22;
	V x24 = rho*(x11 + x23 - T{9.00000000000000}*x13*x13);
	V x25 = x24 + x9;
	V x30 = x10 + x22;
	V x31 = rho*(x27 + x30 - T{9.00000000000000}*x29*x29);
	V x32 = x26 + x31;
	V x36 = rho*(x14 + x30 - x35);
	V x40 = -x16;
	V x41 = x10 + x40;
	V x42 = rho*(x14 + x35 + x39 + x41) - T{72.0000000000000}*f_curr[7];
	V x43 = -x33 - x36 + x42;
	V x45 = T{6.00000000000000}*u[2];
	V x46 = -x45;
	V x47 = x12 + u[2];
	V x48 = rho*(x23 + x46 - T{9.00000000000000}*x47*x47);
	V x49 = x44 + x48;
	V x51 = -u[2];
	V x52 = x51 + u[0];
	V x53 = x22 + x45;
	V x54 = rho*(x27 + x53 - T{9.00000000000000}*x52*x52);
	V x55 = x50 + x54;
	V x57 = u[0] + u[2];
	V x58 = T{9.00000000000000}*(x57*x57);
	V x59 = rho*(x14 + x53 - x58);
	V x60 = x39 + x45;
	V x61 = x14 + x40;
	V x62 = rho*(x58 + x60 + x61) - T{72.0000000000000}*f_curr[3];
	V x63 = -x56 - x59 + x62;
	V x65 = x28 + u[2];
	V x66 = rho*(x30 + x46 - T{9.00000000000000}*x65*x65);
	V x67 = x64 + x66;
	V x69 = x51 + u[1];
	V x70 = rho*(x11 + x53 - T{9.00000000000000}*x69*x69);
	V x71 = x68 + x70;
	V x73 = u[1] + u[2];
	V x74 = T{9.00000000000000}*(x73*x73);
	V x75 = rho*(x30 + x45 - x74);
	V x76 = rho*(x41 + x60 + x74) - T{72.0000000000000}*f_curr[0];
	V x77 = -x72 - x75 + x76;
	V x80 = x16 + T{-2.00000000000000};
	V x81 = x14 + x20 - x79 + x80;
	V x82 = x37 + x61 + x79 + T{2.00000000000000};
	V x83 = -x44 - x48 - x50 - x54 + x63;
	V x84 = -x24 - x26 - x31 + x43 - x9;
	V x86 = x10 + x18 + x80 - x85;
	V x87 = x38 + x41 + x85;
	V x88 = -x64 - x66 - x68 - x70 + x77;
	V x89 = T{6.00000000000000}*x15;
	V x90 = x21 + x45 - x89;
	V x91 = x60 + x89;
	V x92 = T{1}/(T{20.0000000000000}*tau + T{2.00000000000000} + T{5.04537849152229}*simd::sqrt<T>((smagorinsky*smagorinsky)*simd::sqrt<T>((x25 + x32 + x43)*(x25 + x32 + x43) + (x49 + x55 + x63)*(x49 + x55 + x63) + (x67 + x71 + x77)*(x67 + x71 + x77) + T{0.500000000000000}*((-x78*x81 + x78*x82 + x83 + x84 - T{72}*f_curr[10] - T{72}*f_curr[8])*(-x78*x81 + x78*x82 + x83 + x84 - T{72}*f_curr[10] - T{72}*f_curr[8])) + T{0.500000000000000}*((-x78*x86 + x78*x87 + x84 + x88 - T{72}*f_curr[12] - T{72}*f_curr[6])*(-x78*x86 + x78*x87 + x84 + x88 - T{72}*f_curr[12] - T{72}*f_curr[6])) + T{0.500000000000000}*((-x78*x90 + x78*x91 + x83 + x88 - T{72}*f_curr[16] - T{72}*f_curr[2])*(-x78*x90 + x78*x91 + x83 + x88 - T{72}*f_curr[16] - T{72}*f_curr[2]))) + T{0.157134840263677}*(x8*x8)));
	f_next[0] = T{0.333333333333333}*x76*x92 + f_curr[0];
	f_next[1] = -T{0.333333333333333}*x49*x92 + f_curr[1];
	f_next[2] = T{0.666666666666667}*x92*(rho*x91 - T{36.0000000000000}*f_curr[2]) + f_curr[2];
	f_next[3] = T{0.333333333333333}*x62*x92 + f_curr[3];
	f_next[4] = -T{0.333333333333333}*x67*x92 + f_curr[4];
	f_next[5] = -T{0.333333333333333}*x25*x92 + f_curr[5];
	f_next[6] = T{0.666666666666667}*x92*(rho*x87 - T{36.0000000000000}*f_curr[6]) + f_curr[6];
	f_next[7] = T{0.333333333333333}*x42*x92 + f_curr[7];
	f_next[8] = -T{0.666666666666667}*x92*(rho*x81 + T{36.0000000000000}*f_curr[8]) + f_curr[8];
	f_next[9] = -T{4.00000000000000}*x92*(rho*x22 + T{6.00000000000000}*f_curr[9]) + f_curr[9];
	f_next[10] = T{0.666666666666667}*x92*(rho*x82 - T{36.0000000000000}*f_curr[10]) + f_curr[10];
	f_next[11] = -T{0.333333333333333}*x92*(x33 + x36) + f_curr[11];
	f_next[12] = -T{0.666666666666667}*x92*(rho*x86 + T{36.0000000000000}*f_curr[12]) + f_curr[12];
	f_next[13] = -T{0.333333333333333}*x32*x92 + f_curr[13];
	f_next[14] = -T{0.333333333333333}*x71*x92 + f_curr[14];
	f_next[15] = -T{0.333333333333333}*x92*(x56 + x59) + f_curr[15];
	f_next[16] = -T{0.666666666666667}*x92*(rho*x90 + T{36.0000000000000}*f_curr[16]) + f_curr[16];
	f_next[17] = -T{0.333333333333333}*x55*x92 + f_curr[17];
	f_next[18] = -T{0.333333333333333}*x92*(x72 + x75) + f_curr[18];
}

};