diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp')
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp | 258 |
1 files changed, 184 insertions, 74 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp index 0757d17cd..ea5aa7ec0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp @@ -18,6 +18,16 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace, const bool update_with_rho) // Initialize members of base class : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal), + // Initialize the centered finite-order modified k vectors: these are computed + // always with the assumption of centered grids (argument nodal = true), + // for both nodal and staggered simulations + modified_kx_vec_centered(spectral_kspace.getModifiedKComponent(dm,0,norder_x,true)), +#if (AMREX_SPACEDIM==3) + modified_ky_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_y,true)), + modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,2,norder_z,true)), +#else + modified_kz_vec_centered(spectral_kspace.getModifiedKComponent(dm,1,norder_z,true)), +#endif m_v_galilean(v_galilean), m_dt(dt), m_update_with_rho(update_with_rho) @@ -91,8 +101,8 @@ GalileanAlgorithm::pushSpectralFields (SpectralFieldData& f) const const Real ky = modified_ky_arr[j]; const Real kz = modified_kz_arr[k]; #else - constexpr Real ky = 0; - const Real kz = modified_kz_arr[j]; + constexpr Real ky = 0._rt; + const Real kz = modified_kz_arr[j]; #endif // Physical constant c**2 and imaginary unit constexpr Real c2 = PhysConst::c*PhysConst::c; @@ -168,11 +178,14 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp const Box& bx = ba[mfi]; // Extract pointers for the k vectors - const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); + const Real* kx = modified_kx_vec[mfi].dataPtr(); + const Real* kx_c = modified_kx_vec_centered[mfi].dataPtr(); #if (AMREX_SPACEDIM==3) - const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); + const Real* ky = modified_ky_vec[mfi].dataPtr(); + const Real* ky_c = modified_ky_vec_centered[mfi].dataPtr(); #endif - const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + const Real* kz = modified_kz_vec[mfi].dataPtr(); + const Real* kz_c = modified_kz_vec_centered[mfi].dataPtr(); // Extract arrays for the coefficients Array4<Real> C = C_coef[mfi].array(); @@ -194,13 +207,22 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Calculate norm of vector - const Real k_norm = std::sqrt( - std::pow(modified_kx[i], 2) + + const Real knorm = std::sqrt( + std::pow(kx[i], 2) + #if (AMREX_SPACEDIM==3) - std::pow(modified_ky[j], 2) + - std::pow(modified_kz[k], 2)); + std::pow(ky[j], 2) + + std::pow(kz[k], 2)); #else - std::pow(modified_kz[j], 2)); + std::pow(kz[j], 2)); +#endif + // Calculate norm of vector + const Real knorm_c = std::sqrt( + std::pow(kx_c[i], 2) + +#if (AMREX_SPACEDIM==3) + std::pow(ky_c[j], 2) + + std::pow(kz_c[k], 2)); +#else + std::pow(kz_c[j], 2)); #endif // Physical constants c, c**2, and epsilon_0, and imaginary unit constexpr Real c = PhysConst::c; @@ -213,142 +235,222 @@ void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& sp const Real dt3 = dt * dt2; Complex X2_old, X3_old; - // Calculate dot product of k vector with Galilean velocity - const Real kv = modified_kx[i]*vx + + // Calculate dot product of k vector with Galilean velocity: + // this has to be computed always with the centered finite-order modified k + // vectors, in order to work correctly for both nodal and staggered simulations + const Real kv = kx_c[i]*vx + #if (AMREX_SPACEDIM==3) - modified_ky[j]*vy + modified_kz[k]*vz; + ky_c[j]*vy + kz_c[k]*vz; #else - modified_kz[j]*vz; + kz_c[j]*vz; #endif // The coefficients in the following refer to the ones given in equations // (12a)-(12d) of (Lehe et al, PRE 94, 2016), used to update B and E // (equations (11a) and (11b) of the same reference, respectively) - if (k_norm != 0.) { + if (knorm != 0. && knorm_c != 0.) { // Auxiliary coefficients - const Real k2 = k_norm * k_norm; - const Real ck = c * k_norm; - const Real ckdt = ck * dt; - const Complex tmp1 = amrex::exp( I * ckdt); // limit of T2 for nu = 1 - const Complex tmp2 = amrex::exp(- I * ckdt); // limit of T2 for nu = -1 + const Real om = c * knorm; + const Real om2 = om * om; + const Real om3 = om * om2; + const Real om_c = c * knorm_c; + const Real om2_c = om_c * om_c; + const Real om3_c = om_c * om2_c; + const Complex tmp1 = amrex::exp( I * om * dt); + const Complex tmp2 = amrex::exp(- I * om * dt); // See equation (12a) - C (i,j,k) = std::cos(ckdt); - S_ck(i,j,k) = std::sin(ckdt) / ck; + C (i,j,k) = std::cos(om * dt); + S_ck(i,j,k) = std::sin(om * dt) / om; // See equation (12b) - const Real nu = kv / ck; - const Complex theta = amrex::exp( I * 0.5_rt * kv * dt); - const Complex theta_star = amrex::exp(- I * 0.5_rt * kv * dt); + const Real nu = kv / om_c; + const Complex theta = amrex::exp( I * nu * om_c * dt * 0.5_rt); + const Complex theta_star = amrex::exp(- I * nu * om_c * dt * 0.5_rt); // This is exp(i*(k \dot v_gal)*dt) T2(i,j,k) = theta * theta; - if ( (nu != 1.) && (nu != 0.) ) { + if ( (nu != om/om_c) && (nu != -om/om_c) && (nu != 0.) ) { // x1 is the coefficient chi_1 in equation (12c) - Complex x1 = 1._rt / (1._rt - nu*nu) - * (theta_star - C(i,j,k) * theta + I * kv * S_ck(i,j,k) * theta); + Complex x1 = om2_c / (om2 - nu * nu * om2_c) + * (theta_star - theta * C(i,j,k) + I * nu * om_c * theta * S_ck(i,j,k)); // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = theta * x1 / (ep0 * c2 * k2); + X1(i,j,k) = theta * x1 / (ep0 * om2_c); if (update_with_rho) { // X2 multiplies rho_new in the update equation for E // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2); - X3(i,j,k) = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2); + X2(i,j,k) = c2 * (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c) + / (theta_star - theta) / (ep0 * om2_c * om2); + X3(i,j,k) = c2 * (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c) + / (theta_star - theta) / (ep0 * om2_c * om2); } else { // X2_old is the coefficient chi_2 in equation (12d) // X3_old is the coefficient chi_3 in equation (12d) // X2 multiplies (k \dot E) in the update equation for E // X3 multiplies (k \dot J) in the update equation for E - X2_old = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta); - X3_old = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta); - X2(i,j,k) = T2(i,j,k) * (X2_old - X3_old) / k2; - X3(i,j,k) = I * X2_old * (T2(i,j,k) - 1._rt) / (ep0 * k2 * kv); + X2_old = (x1 * om2 - theta * (1._rt - C(i,j,k)) * om2_c) + / (theta_star - theta); + X3_old = (x1 * om2 - theta_star * (1._rt - C(i,j,k)) * om2_c) + / (theta_star - theta); + X2(i,j,k) = c2 * T2(i,j,k) * (X2_old - X3_old) + / (om2_c * om2); + X3(i,j,k) = I * c2 * X2_old * (T2(i,j,k) - 1._rt) + / (ep0 * nu * om3_c * om2); } // X4 multiplies J in the update equation for E - X4(i,j,k) = I * kv * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0; + X4(i,j,k) = I * nu * om_c * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0; } // Limits for nu = 0 if (nu == 0.) { // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * c2 * k2); + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2); if (update_with_rho) { // X2 multiplies rho_new in the update equation for E // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = (1._rt - S_ck(i,j,k) / dt) / (ep0 * k2); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * k2); + X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2); + X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2); } else { // X2 multiplies (k \dot E) in the update equation for E // X3 multiplies (k \dot J) in the update equation for E - X2(i,j,k) = (1._rt - C(i,j,k)) / k2; - X3(i,j,k) = (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * k2); + X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2; + X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2); } // Coefficient multiplying J in update equation for E X4(i,j,k) = - S_ck(i,j,k) / ep0; } - // Limits for nu = 1 - if (nu == 1.) { + // Limits for nu = omega/omega_c + if (nu == om/om_c) { // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2); + X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * om * dt) / (4._rt * ep0 * om2); if (update_with_rho) { // X2 multiplies rho_new in the update equation for E // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * ckdt) - / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); - X3(i,j,k) = (3._rt - 2._rt / tmp1 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * ckdt) - / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + X2(i,j,k) = c2 * (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * om * dt) + / (4._rt * ep0 * om2 * (tmp1 - 1._rt)); + X3(i,j,k) = c2 * (3._rt - 2._rt * tmp2 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * om * dt) + / (4._rt * ep0 * om2 * (tmp1 - 1._rt)); } else { // X2 multiplies (k \dot E) in the update equation for E // X3 multiplies (k \dot J) in the update equation for E - X2(i,j,k) = (1._rt - C(i,j,k)) * tmp1 / k2; - X3(i,j,k) = (2._rt * ckdt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I) - / (4._rt * ep0 * ck * k2); + X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp1 / om2; + X3(i,j,k) = c2 * (2._rt * om * dt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I) + / (4._rt * ep0 * om3); } // Coefficient multiplying J in update equation for E - X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * ckdt) / (4._rt * ep0 * ck); + X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * om * dt) / (4._rt * ep0 * om); } - // Limits for nu = -1 - if (nu == -1.) { + // Limits for nu = -omega/omega_c + if (nu == -om/om_c) { // X1 multiplies i*(k \times J) in the update equation for B - X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2); + X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * om * dt) / (4._rt * ep0 * om2); if (update_with_rho) { // X2 multiplies rho_new in the update equation for E // X3 multiplies rho_old in the update equation for E - X2(i,j,k) = (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * ckdt * tmp1) - / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); - X3(i,j,k) = (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * ckdt * tmp1) - / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + X2(i,j,k) = c2 * (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * om * dt * tmp1) + / (4._rt * ep0 * om2 * (tmp1 - 1._rt)); + X3(i,j,k) = c2 * (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * om * dt * tmp1) + / (4._rt * ep0 * om2 * (tmp1 - 1._rt)); } else { // X2 multiplies (k \dot E) in the update equation for E // X3 multiplies (k \dot J) in the update equation for E - X2(i,j,k) = (1._rt - C(i,j,k)) * tmp2 / k2; - X3(i,j,k) = (2._rt * ckdt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I) - / (4._rt * ep0 * ck * k2); + X2(i,j,k) = c2 * (1._rt - C(i,j,k)) * tmp2 / om2; + X3(i,j,k) = c2 * (2._rt * om * dt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I) + / (4._rt * ep0 * om3); } // Coefficient multiplying J in update equation for E - X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * ckdt) / (4._rt * ep0 * ck); + X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * om * dt) / (4._rt * ep0 * om); + } + } + + // Limits for omega_c = 0 only + else if (knorm != 0. && knorm_c == 0.) { + + const Real om = c * knorm; + const Real om2 = om * om; + + C (i,j,k) = std::cos(om * dt); + S_ck(i,j,k) = std::sin(om * dt) / om; + T2(i,j,k) = 1._rt; + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * om2); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (1._rt - S_ck(i,j,k) / dt) / (ep0 * om2); + X3(i,j,k) = c2 * (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * om2); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = c2 * (1._rt - C(i,j,k)) / om2; + X3(i,j,k) = c2 * (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * om2); } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = - S_ck(i,j,k) / ep0; + + } + + // Limits for omega = 0 only + else if (knorm == 0. && knorm_c != 0.) { + + const Real om_c = c * knorm_c; + const Real om2_c = om_c * om_c; + const Real om3_c = om_c * om2_c; + const Real nu = kv / om_c; + const Complex theta = amrex::exp(I * nu * om_c * dt * 0.5_rt); + + C(i,j,k) = 1._rt; + S_ck(i,j,k) = dt; + T2(i,j,k) = theta * theta; + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (-1._rt + T2(i,j,k) - I * nu * om_c * dt * T2(i,j,k)) + / (ep0 * nu * nu * om2_c); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k) + + 0.5_rt * nu * nu * om2_c * dt * dt * T2(i,j,k)) + / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt)); + X3(i,j,k) = c2 * (1._rt - T2(i,j,k) + I * nu * om_c * dt * T2(i,j,k) + + 0.5_rt * nu * nu * om2_c * dt * dt) + / (ep0 * nu * nu * om2_c * (T2(i,j,k) - 1._rt)); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = c2 * dt * dt * T2(i,j,k) * 0.5_rt; + X3(i,j,k) = c2 * (2._rt * I - 2._rt * nu * om_c * dt * T2(i,j,k) + + I * nu * nu * om2_c * dt * dt * T2(i,j,k)) + / (2._rt * ep0 * nu * nu * nu * om3_c); + } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = I * (T2(i,j,k) - 1._rt) / (ep0 * nu * om_c); } - // Limits for k = 0 - else { + // Limits for omega = 0 and omega_c = 0 + else if (knorm == 0. && knorm_c == 0.) { // Limits of cos(c*k*dt) and sin(c*k*dt)/(c*k) C(i,j,k) = 1._rt; @@ -404,11 +506,14 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data, amrex::Array4<Complex> fields = field_data.fields[mfi].array(); // Extract pointers for the k vectors - const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); + const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); + const amrex::Real* const modified_kx_arr_c = modified_kx_vec_centered[mfi].dataPtr(); #if (AMREX_SPACEDIM==3) - const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); + const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); + const amrex::Real* const modified_ky_arr_c = modified_ky_vec_centered[mfi].dataPtr(); #endif - const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + const amrex::Real* const modified_kz_arr_c = modified_kz_vec_centered[mfi].dataPtr(); // Local copy of member variables before GPU loop const amrex::Real dt = m_dt; @@ -429,13 +534,18 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data, const Complex rho_new = fields(i,j,k,Idx::rho_new); // k vector values, and coefficients - const amrex::Real kx = modified_kx_arr[i]; + const amrex::Real kx = modified_kx_arr[i]; + const amrex::Real kx_c = modified_kx_arr_c[i]; #if (AMREX_SPACEDIM==3) - const amrex::Real ky = modified_ky_arr[j]; - const amrex::Real kz = modified_kz_arr[k]; + const amrex::Real ky = modified_ky_arr[j]; + const amrex::Real kz = modified_kz_arr[k]; + const amrex::Real ky_c = modified_ky_arr_c[j]; + const amrex::Real kz_c = modified_kz_arr_c[k]; #else - constexpr amrex::Real ky = 0._rt; - const amrex::Real kz = modified_kz_arr[j]; + constexpr amrex::Real ky = 0._rt; + const amrex::Real kz = modified_kz_arr[j]; + constexpr amrex::Real ky_c = 0._rt; + const amrex::Real kz_c = modified_kz_arr_c[j]; #endif constexpr Complex I = Complex{0._rt,1._rt}; @@ -445,7 +555,7 @@ GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data, if (k_norm != 0._rt) { const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; - const amrex::Real k_dot_vg = kx * vgx + ky * vgy + kz * vgz; + const amrex::Real k_dot_vg = kx_c * vgx + ky_c * vgy + kz_c * vgz; if ( k_dot_vg != 0._rt ) { |