@@ -30,14 +30,30 @@ real Fluid<Phys>::CheckDivB() {
3030 data->beg [IDIR], data->end [IDIR],
3131 KOKKOS_LAMBDA (int k, int j, int i, real &divBmax) {
3232 [[maybe_unused]] real dB1,dB2,dB3;
33+ [[maybe_unused]] real d1, d2, d3;
34+ [[maybe_unused]] real B1,B2,B3;
3335
3436 dB1=dB2=dB3=ZERO_F;
37+ d1=d2=d3=ZERO_F;
38+ B1=B2=B3=ZERO_F;
39+
3540
3641 D_EXPAND ( dB1=(Ax1 (k,j,i+1 )*Vs (BX1s,k,j,i+1 )-Ax1 (k,j,i)*Vs (BX1s,k,j,i)); ,
3742 dB2=(Ax2 (k,j+1 ,i)*Vs (BX2s,k,j+1 ,i)-Ax2 (k,j,i)*Vs (BX2s,k,j,i)); ,
3843 dB3=(Ax3 (k+1 ,j,i)*Vs (BX3s,k+1 ,j,i)-Ax3 (k,j,i)*Vs (BX3s,k,j,i)); )
3944
40- divBmax=FMAX (FABS (D_EXPAND (dB1, +dB2, +dB3))/dV (k,j,i),divBmax);
45+ D_EXPAND ( d1=0.5 *(Ax1 (k,j,i+1 ) + Ax1 (k,j,i)); ,
46+ d2=0.5 *(Ax2 (k,j+1 ,i) + Ax2 (k,j,i)); ,
47+ d3=0.5 *(Ax3 (k+1 ,j,i) + Ax3 (k,j,i)); )
48+
49+ D_EXPAND ( B1=0.5 *(Vs (BX1s,k,j,i+1 ) + Vs (BX1s,k,j,i)); ,
50+ B2=0.5 *(Vs (BX2s,k,j+1 ,i) + Vs (BX2s,k,j,i)); ,
51+ B3=0.5 *(Vs (BX3s,k+1 ,j,i) + Vs (BX3s,k,j,i)); )
52+
53+ real amplitude = 1e-40 ;
54+ amplitude += D_EXPAND ( std::fabs (B1)*d1, + std::fabs (B2)*d2, + std::fabs (B3)*d3 );
55+
56+ divBmax=FMAX (FABS (D_EXPAND (dB1, +dB2, +dB3))/amplitude,divBmax);
4157 },
4258 Kokkos::Max<real>(divB) // reduction
4359 );
0 commit comments