18
18
// 4. Fourier-space projection methods for the computation of the
19
19
// curl and divergence of the velocity field
20
20
//
21
- // Author: Julian Adamek (Université de Genève & Observatoire de Paris & Queen Mary University of London)
21
+ // Author: Julian Adamek (Université de Genève & Observatoire de Paris & Queen Mary University of London & Universität Zürich )
22
22
//
23
- // Last modified: April 2019
23
+ // Last modified: August 2022
24
24
//
25
25
// ////////////////////////
26
26
@@ -580,7 +580,8 @@ Real update_q(double dtau, double dx, part_simple * part, double * ref_dist, par
580
580
Real pgradB[3 ]={0 ,0 ,0 };
581
581
Real v2 = (*part).vel [0 ] * (*part).vel [0 ] + (*part).vel [1 ] * (*part).vel [1 ] + (*part).vel [2 ] * (*part).vel [2 ];
582
582
Real e2 = v2 + params[0 ] * params[0 ];
583
-
583
+
584
+ #if GRADIENT_ORDER == 1
584
585
gradphi[0 ] = (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (phi (xphi+0 ) - phi (xphi));
585
586
gradphi[1 ] = (1 .-ref_dist[0 ]) * (1 .-ref_dist[2 ]) * (phi (xphi+1 ) - phi (xphi));
586
587
gradphi[2 ] = (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * (phi (xphi+2 ) - phi (xphi));
@@ -593,6 +594,21 @@ Real update_q(double dtau, double dx, part_simple * part, double * ref_dist, par
593
594
gradphi[0 ] += ref_dist[1 ] * ref_dist[2 ] * (phi (xphi+2 +1 +0 ) - phi (xphi+2 +1 ));
594
595
gradphi[1 ] += ref_dist[0 ] * ref_dist[2 ] * (phi (xphi+2 +1 +0 ) - phi (xphi+2 +0 ));
595
596
gradphi[2 ] += ref_dist[0 ] * ref_dist[1 ] * (phi (xphi+2 +1 +0 ) - phi (xphi+1 +0 ));
597
+ #elif GRADIENT_ORDER == 2
598
+ for (int i=0 ; i<3 ; i++)
599
+ {
600
+ gradphi[i] = 0.5 * (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (phi (xphi+i) - phi (xphi-i));
601
+ gradphi[i] += 0.5 * ref_dist[0 ] * (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (phi (xphi+i+0 ) - phi (xphi-i+0 ));
602
+ gradphi[i] += 0.5 * (1 .-ref_dist[0 ]) * ref_dist[1 ] * (1 .-ref_dist[2 ]) * (phi (xphi+i+1 ) - phi (xphi-i+1 ));
603
+ gradphi[i] += 0.5 * ref_dist[0 ] * ref_dist[1 ] * (1 .-ref_dist[2 ]) * (phi (xphi+i+1 +0 ) - phi (xphi-i+1 +0 ));
604
+ gradphi[i] += 0.5 * (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * ref_dist[2 ] * (phi (xphi+2 +i) - phi (xphi+2 -i));
605
+ gradphi[i] += 0.5 * ref_dist[0 ] * (1 .-ref_dist[1 ]) * ref_dist[2 ] * (phi (xphi+2 +i+0 ) - phi (xphi+2 -i+0 ));
606
+ gradphi[i] += 0.5 * (1 .-ref_dist[0 ]) * ref_dist[1 ] * ref_dist[2 ] * (phi (xphi+2 +i+1 ) - phi (xphi+2 -i+1 ));
607
+ gradphi[i] += 0.5 * ref_dist[0 ] * ref_dist[1 ] * ref_dist[2 ] * (phi (xphi+2 +i+1 +0 ) - phi (xphi+2 -i+1 +0 ));
608
+ }
609
+ #else
610
+ #error GRADIENT_ORDER must be set to 1 or 2
611
+ #endif
596
612
597
613
gradphi[0 ] *= (v2 + e2 ) / e2 ;
598
614
gradphi[1 ] *= (v2 + e2 ) / e2 ;
@@ -698,7 +714,8 @@ Real update_q_Newton(double dtau, double dx, part_simple * part, double * ref_di
698
714
#define xchi (sites[1 ])
699
715
700
716
Real gradpsi[3 ]={0 ,0 ,0 };
701
-
717
+
718
+ #if GRADIENT_ORDER == 1
702
719
gradpsi[0 ] = (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (psi (xpsi+0 ) - psi (xpsi));
703
720
gradpsi[1 ] = (1 .-ref_dist[0 ]) * (1 .-ref_dist[2 ]) * (psi (xpsi+1 ) - psi (xpsi));
704
721
gradpsi[2 ] = (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * (psi (xpsi+2 ) - psi (xpsi));
@@ -711,6 +728,21 @@ Real update_q_Newton(double dtau, double dx, part_simple * part, double * ref_di
711
728
gradpsi[0 ] += ref_dist[1 ] * ref_dist[2 ] * (psi (xpsi+2 +1 +0 ) - psi (xpsi+2 +1 ));
712
729
gradpsi[1 ] += ref_dist[0 ] * ref_dist[2 ] * (psi (xpsi+2 +1 +0 ) - psi (xpsi+2 +0 ));
713
730
gradpsi[2 ] += ref_dist[0 ] * ref_dist[1 ] * (psi (xpsi+2 +1 +0 ) - psi (xpsi+1 +0 ));
731
+ #elif GRADIENT_ORDER == 2
732
+ for (int i=0 ; i<3 ; i++)
733
+ {
734
+ gradpsi[i] = 0.5 * (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (psi (xpsi+i) - psi (xpsi-i));
735
+ gradpsi[i] += 0.5 * ref_dist[0 ] * (1 .-ref_dist[1 ]) * (1 .-ref_dist[2 ]) * (psi (xpsi+i+0 ) - psi (xpsi-i+0 ));
736
+ gradpsi[i] += 0.5 * (1 .-ref_dist[0 ]) * ref_dist[1 ] * (1 .-ref_dist[2 ]) * (psi (xpsi+i+1 ) - psi (xpsi-i+1 ));
737
+ gradpsi[i] += 0.5 * ref_dist[0 ] * ref_dist[1 ] * (1 .-ref_dist[2 ]) * (psi (xpsi+i+1 +0 ) - psi (xpsi-i+1 +0 ));
738
+ gradpsi[i] += 0.5 * (1 .-ref_dist[0 ]) * (1 .-ref_dist[1 ]) * ref_dist[2 ] * (psi (xpsi+2 +i) - psi (xpsi+2 -i));
739
+ gradpsi[i] += 0.5 * ref_dist[0 ] * (1 .-ref_dist[1 ]) * ref_dist[2 ] * (psi (xpsi+2 +i+0 ) - psi (xpsi+2 -i+0 ));
740
+ gradpsi[i] += 0.5 * (1 .-ref_dist[0 ]) * ref_dist[1 ] * ref_dist[2 ] * (psi (xpsi+2 +i+1 ) - psi (xpsi+2 -i+1 ));
741
+ gradpsi[i] += 0.5 * ref_dist[0 ] * ref_dist[1 ] * ref_dist[2 ] * (psi (xpsi+2 +i+1 +0 ) - psi (xpsi+2 -i+1 +0 ));
742
+ }
743
+ #else
744
+ #error GRADIENT_ORDER must be set to 1 or 2
745
+ #endif
714
746
715
747
if (nfield>=2 && fields[1 ] != NULL )
716
748
{
0 commit comments