106{
107 using V = typename CELL::value_t;
109
110 auto prevCell = cell.template getField<PREV_CELL>();
111
112 for (unsigned i=0; i < missing.size(); ++i) {
113 cell[missing[i]] = prevCell[i];
114 }
115
116 V rho0, u0[3];
117 V rho1, u1[3];
118 V rho2, u2[3];
119
120 cell.computeRhoU(rho0, u0);
121
122 static_assert(direction == 0 || direction == 1 || direction ==2,
123 "Direction must be one of 0, 1 or 2");
124 if constexpr (direction == 0) {
125 cell.neighbor({-orientation ,0,0}).computeRhoU(rho1, u1);
126 cell.neighbor({-orientation*2,0,0}).computeRhoU(rho2, u2);
127 }
128 else if constexpr (direction == 1) {
129 cell.neighbor({0,-orientation ,0}).computeRhoU(rho1, u1);
130 cell.neighbor({0,-orientation*2,0}).computeRhoU(rho2, u2);
131 }
132 else if constexpr (direction == 2) {
133 cell.neighbor({0,0,-orientation }).computeRhoU(rho1, u1);
134 cell.neighbor({0,0,-orientation*2}).computeRhoU(rho2, u2);
135 }
136
137 V uDelta[3];
138 V uAverage = rho0*u0[direction];
139
140
141
142
143
144
145
146 uDelta[0] = -uAverage*0.5*(3*rho0*u0[0]-4*rho1*u1[0]+rho2*u2[0]);
147 uDelta[1] = -uAverage*0.5*(3*rho0*u0[1]-4*rho1*u1[1]+rho2*u2[1]);
148 uDelta[2] = -uAverage*0.5*(3*rho0*u0[2]-4*rho1*u1[2]+rho2*u2[2]);
149
150 for (unsigned i=0; i < missing.size(); ++i) {
151 auto iPop = missing[i];
152 prevCell[i] = cell[iPop]
153 + descriptors::invCs2<V,DESCRIPTOR>()*descriptors::t<V,DESCRIPTOR>(iPop)
154 * ( uDelta[0]*descriptors::c<DESCRIPTOR>(iPop,0)
155 + uDelta[1]*descriptors::c<DESCRIPTOR>(iPop,1)
156 + uDelta[2]*descriptors::c<DESCRIPTOR>(iPop,2));
157 }
158
159 cell.template setField<PREV_CELL>(prevCell);
160}
constexpr auto populationsContributingToVelocity() any_platform
Return array of population indices where c[iVel] == value.