180 bool schemeSwitch =
false;
183 T omega =
parameters.template get<descriptors::OMEGA>();
184 auto v = cell.template getField<descriptors::VELOCITY>();
185 auto a = cell.template getField<descriptors::G>();
190 constexpr auto unknownIndices = util::subIndexOutgoing3DonEdges<DESCRIPTOR,Plane,Normal1,Normal2>();
191 assert(unknownIndices.size() == 2);
192 auto c0=descriptors::c<DESCRIPTOR>(unknownIndices[0]);
193 auto c1=descriptors::c<DESCRIPTOR>(unknownIndices[1]);
199 T w = descriptors::t<T,DESCRIPTOR>(1);
200 T gamma = -descriptors::invCs2<T,DESCRIPTOR>() * omega;
201 T b1 = a1-a2*gamma*(n[0]*v[0]+n[1]*v[1]+n[2]*v[2]);
205 T sum1 = cell.computeRho();
207 for (
unsigned iPop : unknownIndices) {
208 sum1 -= (cell[iPop]+w + cell[descriptors::opposite<DESCRIPTOR>(iPop)]+w);
209 sum2 += cell[descriptors::opposite<DESCRIPTOR>(iPop)] + w;
211 T beta = (a3-b1*sum1+2*b2*sum2)/(unknownIndices.size()*w*(b1+b2));
213 for (
unsigned iPop : unknownIndices) {
214 cell[iPop] = w * beta - (cell[descriptors::opposite<DESCRIPTOR>(iPop)] + w) - w;
218 T w = descriptors::t<T,DESCRIPTOR>(1);
219 T chi = (a2!=0.) ? -1/ (descriptors::invCs2<T,DESCRIPTOR>()*a2*omega) : (1/omega)/(1/omega-0.5);
221 T NdotV = n[0]*v[0] + n[1]*v[1] + n[2]*v[2];
222 T A = unknownIndices.
size()-1;
223 T alpha = (-chi*k+NdotV+2*w*(1-A)) / (chi*k-NdotV+2*w*(1+A));
227 for (
unsigned iPop : unknownIndices) {
228 B += 2*(cell[descriptors::opposite<DESCRIPTOR>(iPop)]+w);
232 for (
unsigned iPop : unknownIndices) {
233 B -= 2*(cell[descriptors::opposite<DESCRIPTOR>(iPop)]+w);
234 T beta = (2*w*chi*a3+2*w*B) / (chi*k-NdotV+2*w*(1+A));
235 cell[iPop] = alpha * (cell[descriptors::opposite<DESCRIPTOR>(iPop)] + w) + beta - w;
236 B += 2*(cell[descriptors::opposite<DESCRIPTOR>(iPop)]+w);