30 Real bx, Real by, Real bz)
noexcept
36 Real hp, hm,
scale, out;
38 hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
39 if (levset(i+1,j,k) < Real(0.0)) {
40 tmp =
x(i+1,j,k) -
x(i,j,k);
42 tmp = (xeb(i+1,j,k) -
x(i,j,k)) / hp;
45 hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
46 if (levset(i-1,j,k) < Real(0.0)) {
47 tmp +=
x(i-1,j,k) -
x(i,j,k);
49 tmp += (xeb(i-1,j,k) -
x(i,j,k)) / hm;
52 out = tmp * bx * Real(2.0) / (hp+hm);
55 hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
56 if (levset(i,j+1,k) < Real(0.0)) {
57 tmp =
x(i,j+1,k) -
x(i,j,k);
59 tmp = (xeb(i,j+1,k) -
x(i,j,k)) / hp;
62 hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
63 if (levset(i,j-1,k) < Real(0.0)) {
64 tmp +=
x(i,j-1,k) -
x(i,j,k);
66 tmp += (xeb(i,j-1,k) -
x(i,j,k)) / hm;
69 out += tmp * by * Real(2.0) / (hp+hm);
72 hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
73 if (levset(i,j,k+1) < Real(0.0)) {
74 tmp =
x(i,j,k+1) -
x(i,j,k);
76 tmp = (xeb(i,j,k+1) -
x(i,j,k)) / hp;
79 hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
80 if (levset(i,j,k-1) < Real(0.0)) {
81 tmp +=
x(i,j,k-1) -
x(i,j,k);
83 tmp += (xeb(i,j,k-1) -
x(i,j,k)) / hm;
86 out += tmp * bz * Real(2.0) / (hp+hm);
141 int redblack)
noexcept
143 if ((i+j+k+redblack)%2 == 0) {
149 hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
150 if (levset(i+1,j,k) < Real(0.0)) {
154 tmp0 = Real(-1.0) / hp;
158 hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
159 if (levset(i-1,j,k) < Real(0.0)) {
163 tmp0 += Real(-1.0) / hm;
166 Real gamma = tmp0 * (bx * Real(2.0) / (hp+hm));
167 Real rho = tmp1 * (bx * Real(2.0) / (hp+hm));
170 hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
171 if (levset(i,j+1,k) < Real(0.0)) {
175 tmp0 = Real(-1.0) / hp;
179 hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
180 if (levset(i,j-1,k) < Real(0.0)) {
184 tmp0 += Real(-1.0) / hm;
187 gamma += tmp0 * (by * Real(2.0) / (hp+hm));
188 rho += tmp1 * (by * Real(2.0) / (hp+hm));
191 hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
192 if (levset(i,j,k+1) < Real(0.0)) {
196 tmp0 = Real(-1.0) / hp;
200 hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
201 if (levset(i,j,k-1) < Real(0.0)) {
205 tmp0 += Real(-1.0) / hm;
208 gamma += tmp0 * (bz * Real(2.0) / (hp+hm));
209 rho += tmp1 * (bz * Real(2.0) / (hp+hm));
212 Real Ax = rho + gamma*
x(i,j,k);
213 constexpr Real omega = Real(1.25);
214 x(i,j,k) += (rhs(i,j,k) - Ax*
scale) * (omega / (gamma*
scale));
338 Real bx, Real by, Real bz)
noexcept
341 y(i,j,k) = Real(0.0);
344 Real hp, hm,
scale, out;
346 sigma = (sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
347 sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
348 sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
349 sig(i ,j ,k ) * vfrc(i ,j ,k ))
350 / ( vfrc(i ,j-1,k-1) +
354 hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
355 if (levset(i+1,j,k) < Real(0.0)) {
356 tmp = sigma*(
x(i+1,j,k) -
x(i,j,k));
358 tmp = sigma*((xeb(i+1,j,k) -
x(i,j,k)) / hp);
361 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
362 sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
363 sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
364 sig(i-1,j ,k ) * vfrc(i-1,j ,k ))
365 / ( vfrc(i-1,j-1,k-1) +
369 hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
370 if (levset(i-1,j,k) < Real(0.0)) {
371 tmp += sigma*(
x(i-1,j,k) -
x(i,j,k));
373 tmp += sigma*((xeb(i-1,j,k) -
x(i,j,k)) / hm);
376 out = tmp * bx * Real(2.0) / (hp+hm);
379 sigma = (sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
380 sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
381 sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
382 sig(i ,j ,k ) * vfrc(i ,j ,k ))
383 / ( vfrc(i-1,j ,k-1) +
387 hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
388 if (levset(i,j+1,k) < Real(0.0)) {
389 tmp = sigma*(
x(i,j+1,k) -
x(i,j,k));
391 tmp = sigma*((xeb(i,j+1,k) -
x(i,j,k)) / hp);
394 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
395 sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
396 sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
397 sig(i ,j-1,k ) * vfrc(i ,j-1,k ))
398 / ( vfrc(i-1,j-1,k-1) +
402 hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
403 if (levset(i,j-1,k) < Real(0.0)) {
404 tmp += sigma*(
x(i,j-1,k) -
x(i,j,k));
406 tmp += sigma*((xeb(i,j-1,k) -
x(i,j,k)) / hm);
409 out += tmp * by * Real(2.0) / (hp+hm);
412 sigma = (sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
413 sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
414 sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
415 sig(i ,j ,k ) * vfrc(i ,j ,k ))
416 / ( vfrc(i-1,j-1,k ) +
420 hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
421 if (levset(i,j,k+1) < Real(0.0)) {
422 tmp = sigma*(
x(i,j,k+1) -
x(i,j,k));
424 tmp = sigma*((xeb(i,j,k+1) -
x(i,j,k)) / hp);
427 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
428 sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
429 sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
430 sig(i ,j ,k-1) * vfrc(i ,j ,k-1))
431 / ( vfrc(i-1,j-1,k-1) +
435 hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
436 if (levset(i,j,k-1) < Real(0.0)) {
437 tmp += sigma*(
x(i,j,k-1) -
x(i,j,k));
439 tmp += sigma*((xeb(i,j,k-1) -
x(i,j,k)) / hm);
442 out += tmp * bz * Real(2.0) / (hp+hm);
445 y(i,j,k) = out*
scale;
495 Real bx, Real by, Real bz,
int redblack)
noexcept
497 if ((i+j+k+redblack)%2 == 0) {
501 Real tmp0, tmp1, sigma;
504 sigma = (sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
505 sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
506 sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
507 sig(i ,j ,k ) * vfrc(i ,j ,k ))
508 / ( vfrc(i ,j-1,k-1) +
512 hp = (ecx(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecx(i,j,k));
513 if (levset(i+1,j,k) < Real(0.0)) {
514 tmp0 = sigma*Real(-1.0);
515 tmp1 = sigma*
x(i+1,j,k);
517 tmp0 = sigma*Real(-1.0) / hp;
521 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
522 sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
523 sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
524 sig(i-1,j ,k ) * vfrc(i-1,j ,k ))
525 / ( vfrc(i-1,j-1,k-1) +
529 hm = (ecx(i-1,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecx(i-1,j,k));
530 if (levset(i-1,j,k) < Real(0.0)) {
531 tmp0 += sigma*Real(-1.0);
532 tmp1 += sigma*
x(i-1,j,k);
534 tmp0 += sigma*Real(-1.0) / hm;
537 Real gamma = tmp0 * (bx * Real(2.0) / (hp+hm));
538 Real rho = tmp1 * (bx * Real(2.0) / (hp+hm));
541 sigma = (sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
542 sig(i ,j ,k-1) * vfrc(i ,j ,k-1) +
543 sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
544 sig(i ,j ,k ) * vfrc(i ,j ,k ))
545 / ( vfrc(i-1,j ,k-1) +
549 hp = (ecy(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.)*ecy(i,j,k));
550 if (levset(i,j+1,k) < Real(0.0)) {
551 tmp0 = sigma*Real(-1.0);
552 tmp1 = sigma*
x(i,j+1,k);
554 tmp0 = sigma*Real(-1.0) / hp;
558 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
559 sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
560 sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
561 sig(i ,j-1,k ) * vfrc(i ,j-1,k ))
562 / ( vfrc(i-1,j-1,k-1) +
566 hm = (ecy(i,j-1,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecy(i,j-1,k));
567 if (levset(i,j-1,k) < Real(0.0)) {
568 tmp0 += sigma*Real(-1.0);
569 tmp1 += sigma*
x(i,j-1,k);
571 tmp0 += sigma*Real(-1.0) / hm;
574 gamma += tmp0 * (by * Real(2.0) / (hp+hm));
575 rho += tmp1 * (by * Real(2.0) / (hp+hm));
578 sigma = (sig(i-1,j-1,k ) * vfrc(i-1,j-1,k ) +
579 sig(i ,j-1,k ) * vfrc(i ,j-1,k ) +
580 sig(i-1,j ,k ) * vfrc(i-1,j ,k ) +
581 sig(i ,j ,k ) * vfrc(i ,j ,k ))
582 / ( vfrc(i-1,j-1,k ) +
586 hp = (ecz(i,j,k) == Real(1.0)) ? Real(1.0) : (Real(1.0)+Real(2.0)*ecz(i,j,k));
587 if (levset(i,j,k+1) < Real(0.0)) {
588 tmp0 = sigma*Real(-1.0);
589 tmp1 = sigma*
x(i,j,k+1);
591 tmp0 = sigma*Real(-1.0) / hp;
595 sigma = (sig(i-1,j-1,k-1) * vfrc(i-1,j-1,k-1) +
596 sig(i ,j-1,k-1) * vfrc(i ,j-1,k-1) +
597 sig(i-1,j ,k-1) * vfrc(i-1,j ,k-1) +
598 sig(i ,j ,k-1) * vfrc(i ,j ,k-1))
599 / ( vfrc(i-1,j-1,k-1) +
603 hm = (ecz(i,j,k-1) == Real(1.0)) ? Real(1.0) : (Real(1.0)-Real(2.)*ecz(i,j,k-1));
604 if (levset(i,j,k-1) < Real(0.0)) {
605 tmp0 += sigma*Real(-1.0);
606 tmp1 += sigma*
x(i,j,k-1);
608 tmp0 += sigma*Real(-1.0) / hm;
611 gamma += tmp0 * (bz * Real(2.0) / (hp+hm));
612 rho += tmp1 * (bz * Real(2.0) / (hp+hm));
615 Real Ax = rho + gamma*
x(i,j,k);
616 constexpr Real omega = Real(1.25);
617 x(i,j,k) += (rhs(i,j,k) - Ax*
scale) * (omega / (gamma*
scale));