176 {
177
178 using V = typename CELLS::template value_t<names::A>::value_t;
179 using DESCRIPTOR = typename CELLS::template value_t<names::A>::descriptor_t;
180
181
182 auto& cellA = cells.template get<names::A>();
183
184
185 auto& cellB = cells.template get<names::B>();
186
187 V u[2] { };
188 V phi = cellA.computeRho();
189 V rho = cellB.computeRho();
190
191
192
193
194
195
196
197
198
199
200
201
202 auto cellA1 = cellA.neighbor(descriptors::c<DESCRIPTOR>(1));
203 auto cellA2 = cellA.neighbor(descriptors::c<DESCRIPTOR>(2));
204 auto cellA3 = cellA.neighbor(descriptors::c<DESCRIPTOR>(3));
205
206 auto cellA7 = cellA.neighbor(descriptors::c<DESCRIPTOR>(7));
207 auto cellA6 = cellA.neighbor(descriptors::c<DESCRIPTOR>(6));
208 auto cellA5 = cellA.neighbor(descriptors::c<DESCRIPTOR>(5));
209
210
211 V gradMuPhiX = 1./12. * (
212 - cellA1.template getField<descriptors::CHEM_POTENTIAL>()
213 - 4. * cellA2.template getField<descriptors::CHEM_POTENTIAL>()
214 - cellA3.template getField<descriptors::CHEM_POTENTIAL>()
215 + cellA7.template getField<descriptors::CHEM_POTENTIAL>()
216 + 4. * cellA6.template getField<descriptors::CHEM_POTENTIAL>()
217 + cellA5.template getField<descriptors::CHEM_POTENTIAL>());
218
219 auto cellA8 = cellA.neighbor(descriptors::c<DESCRIPTOR>(8));
220 auto cellA4 = cellA.neighbor(descriptors::c<DESCRIPTOR>(4));
221
222
223 V gradMuPhiY = 1./12. * (
224 + cellA1.template getField<descriptors::CHEM_POTENTIAL>()
225 + 4. * cellA8.template getField<descriptors::CHEM_POTENTIAL>()
226 + cellA7.template getField<descriptors::CHEM_POTENTIAL>()
227 - cellA3.template getField<descriptors::CHEM_POTENTIAL>()
228 - 4. * cellA4.template getField<descriptors::CHEM_POTENTIAL>()
229 - cellA5.template getField<descriptors::CHEM_POTENTIAL>());
230
231
232
233
234
235
236
237
238
239
240
241 auto cellB1 = cellB.neighbor(descriptors::c<DESCRIPTOR>(1));
242 auto cellB2 = cellB.neighbor(descriptors::c<DESCRIPTOR>(2));
243 auto cellB3 = cellB.neighbor(descriptors::c<DESCRIPTOR>(3));
244
245 auto cellB7 = cellB.neighbor(descriptors::c<DESCRIPTOR>(7));
246 auto cellB6 = cellB.neighbor(descriptors::c<DESCRIPTOR>(6));
247 auto cellB5 = cellB.neighbor(descriptors::c<DESCRIPTOR>(5));
248
249
250 V gradMuRhoX = 1./12. * (
251 - cellB1.template getField<descriptors::CHEM_POTENTIAL>()
252 - 4. * cellB2.template getField<descriptors::CHEM_POTENTIAL>()
253 - cellB3.template getField<descriptors::CHEM_POTENTIAL>()
254 + cellB7.template getField<descriptors::CHEM_POTENTIAL>()
255 + 4. * cellB6.template getField<descriptors::CHEM_POTENTIAL>()
256 + cellB5.template getField<descriptors::CHEM_POTENTIAL>());
257
258 auto cellB8 = cellB.neighbor(descriptors::c<DESCRIPTOR>(8));
259 auto cellB4 = cellB.neighbor(descriptors::c<DESCRIPTOR>(4));
260
261
262 V gradMuRhoY = 1./12. * (
263 + cellB1.template getField<descriptors::CHEM_POTENTIAL>()
264 + 4. * cellB8.template getField<descriptors::CHEM_POTENTIAL>()
265 + cellB7.template getField<descriptors::CHEM_POTENTIAL>()
266 - cellB3.template getField<descriptors::CHEM_POTENTIAL>()
267 - 4. * cellB4.template getField<descriptors::CHEM_POTENTIAL>()
268 - cellB5.template getField<descriptors::CHEM_POTENTIAL>());
269
270
271
272
273 if constexpr (CELLS::map_t::keys_t::template contains<names::C>()){
274 auto& cellC = cells.template get<names::C>();
275 V psi = cellC.computeRho();
276
277 auto cellC1 = cellC.neighbor(descriptors::c<DESCRIPTOR>(1));
278 auto cellC2 = cellC.neighbor(descriptors::c<DESCRIPTOR>(2));
279 auto cellC3 = cellC.neighbor(descriptors::c<DESCRIPTOR>(3));
280
281 auto cellC7 = cellC.neighbor(descriptors::c<DESCRIPTOR>(7));
282 auto cellC6 = cellC.neighbor(descriptors::c<DESCRIPTOR>(6));
283 auto cellC5 = cellC.neighbor(descriptors::c<DESCRIPTOR>(5));
284
285
286
287 V gradMuPsiX = 1./12. * (
288 - cellC1.template getField<descriptors::CHEM_POTENTIAL>()
289 - 4 * cellC2.template getField<descriptors::CHEM_POTENTIAL>()
290 - cellC3.template getField<descriptors::CHEM_POTENTIAL>()
291 + cellC7.template getField<descriptors::CHEM_POTENTIAL>()
292 + 4 * cellC6.template getField<descriptors::CHEM_POTENTIAL>()
293 + cellC5.template getField<descriptors::CHEM_POTENTIAL>());
294
295 auto cellC8 = cellC.neighbor(descriptors::c<DESCRIPTOR>(8));
296 auto cellC4 = cellC.neighbor(descriptors::c<DESCRIPTOR>(4));
297
298
299 V gradMuPsiY = 1./12. * (
300 + cellC1.template getField<descriptors::CHEM_POTENTIAL>()
301 + 4 * cellC8.template getField<descriptors::CHEM_POTENTIAL>()
302 + cellC7.template getField<descriptors::CHEM_POTENTIAL>()
303 - cellC3.template getField<descriptors::CHEM_POTENTIAL>()
304 - 4 * cellC4.template getField<descriptors::CHEM_POTENTIAL>()
305 - cellC5.template getField<descriptors::CHEM_POTENTIAL>());
306
307 cellB.template setField<descriptors::FORCE>({
308 - rho*gradMuRhoX - phi*gradMuPhiX - psi*gradMuPsiX,
309 - rho*gradMuRhoY - phi*gradMuPhiY - psi*gradMuPsiY
310 });
311
312 cellB.computeU(u);
313 cellA.template setField<descriptors::FORCE>(u);
314 cellC.template setField<descriptors::FORCE>(u);
315
316 } else {
317
318 cellB.template setField<descriptors::FORCE>({
319 - rho*gradMuRhoX - phi*gradMuPhiX,
320 - rho*gradMuRhoY - phi*gradMuPhiY
321 });
322
323 cellB.computeU(u);
324 cellA.template setField<descriptors::FORCE>(u);
325 }
326
327 }