From a0870a3063e9e655e281cc31e2d1b6580294696e Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Thu, 4 Mar 2021 13:00:13 -0800 Subject: Implement averaged algo on staggered grids & merge spectral classes (#1544) * Refactor and clean up some spectral classes * Abort when current correction or Vay deposition are not implemented * Implement general equations for averaged Galilean * Allocate averaged MultiFabs also when aux_is_nodal=1 and do_nodal=0 * Allocate +ngextra guard cells also for averaged MultiFabs * Make alias MultiFabs for averaged aux data * With averaging, interpolate from avg_fp (not fp) to aux * Fix some limits of the coefficients * Fix bug causing NaNs in spectral coefficients * Add 2D CI test with same analysis as nodal test * Add 3D CI test with same analysis as nodal test * Add limit that was not covered (knorm=0 && knorm_c!=0 && nu=0) * Allocate T2_coef only if Galilean algorithm is used * Allocate X4_coef only if Galilean algorithm is used * Remove extra ghost cell from 'avg_fp' MultiFabs --- .../SpectralAlgorithms/AvgGalileanAlgorithm.cpp | 376 --------------------- 1 file changed, 376 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp deleted file mode 100644 index ea381235e..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp +++ /dev/null @@ -1,376 +0,0 @@ -#include "FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H" -#include "Utils/WarpXConst.H" -#include - -using namespace amrex; - -/* \brief Initialize coefficients for the update equation */ -AvgGalileanAlgorithm::AvgGalileanAlgorithm(const SpectralKSpace& spectral_kspace, - const DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const bool nodal, - const amrex::Array& v_galilean, - const Real dt) - // Initialize members of base class - : SpectralBaseAlgorithm( spectral_kspace, dm, - norder_x, norder_y, norder_z, nodal ) -{ - const BoxArray& ba = spectral_kspace.spectralspace_ba; - - // Allocate the arrays of coefficients - C_coef = SpectralRealCoefficients(ba, dm, 1, 0); - S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); - - C1_coef = SpectralRealCoefficients(ba, dm, 1, 0); - S1_coef = SpectralRealCoefficients(ba, dm, 1, 0); - C3_coef = SpectralRealCoefficients(ba, dm, 1, 0); - S3_coef = SpectralRealCoefficients(ba, dm, 1, 0); - - Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Psi3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - - X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - - A1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - A2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - - Rhoold_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Rhonew_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - Jcoef_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - - InitializeSpectralCoefficients(spectral_kspace, dm, v_galilean, dt); - -} - -void AvgGalileanAlgorithm::InitializeSpectralCoefficients( - const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const Array& v_galilean, - const amrex::Real dt) -{ - const BoxArray& ba = spectral_kspace.spectralspace_ba; - // Fill them with the right values: - // Loop over boxes and allocate the corresponding coefficients - // for each box owned by the local MPI proc - for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ - - const Box& bx = ba[mfi]; - - // Extract pointers for the k vectors - const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); -#if (AMREX_SPACEDIM==3) - const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); - // Extract arrays for the coefficients - Array4 C = C_coef[mfi].array(); - Array4 S_ck = S_ck_coef[mfi].array(); - Array4 C1 = C1_coef[mfi].array(); - Array4 S1 = S1_coef[mfi].array(); - Array4 C3 = C3_coef[mfi].array(); - Array4 S3 = S3_coef[mfi].array(); - - Array4 Psi1 = Psi1_coef[mfi].array(); - Array4 Psi2 = Psi2_coef[mfi].array(); - Array4 Psi3 = Psi3_coef[mfi].array(); - Array4 X1 = X1_coef[mfi].array(); - Array4 X2 = X2_coef[mfi].array(); - Array4 X3 = X3_coef[mfi].array(); - Array4 X4 = X4_coef[mfi].array(); - Array4 Theta2 = Theta2_coef[mfi].array(); - Array4 A1 = A1_coef[mfi].array(); - Array4 A2 = A2_coef[mfi].array(); - - Array4 CRhoold = Rhoold_coef[mfi].array(); - Array4 CRhonew = Rhonew_coef[mfi].array(); - Array4 Jcoef = Jcoef_coef[mfi].array(); - // Extract reals (for portability on GPU) - Real vx = v_galilean[0]; -#if (AMREX_SPACEDIM==3) - Real vy = v_galilean[1]; -#endif - Real vz = v_galilean[2]; - - // Loop over indices within one box - 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) + -#if (AMREX_SPACEDIM==3) - std::pow(modified_ky[j], 2) + - std::pow(modified_kz[k], 2)); -#else - std::pow(modified_kz[j], 2)); -#endif - - // Calculate coefficients - constexpr Real c = PhysConst::c; - constexpr Real c2 = PhysConst::c*PhysConst::c; - constexpr Real ep0 = PhysConst::ep0; - const Complex I{0.,1.}; - if (k_norm != 0){ - - C(i,j,k) = std::cos(c*k_norm*dt); - S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); - - C1(i,j,k) = std::cos(0.5_rt*c*k_norm*dt); - S1(i,j,k) = std::sin(0.5_rt*c*k_norm*dt); - C3(i,j,k) = std::cos(1.5_rt*c*k_norm*dt); - S3(i,j,k) = std::sin(1.5_rt*c*k_norm*dt); - - // Calculate dot product with galilean velocity - const Real kv = modified_kx[i]*vx + -#if (AMREX_SPACEDIM==3) - modified_ky[j]*vy + - modified_kz[k]*vz; -#else - modified_kz[j]*vz; -#endif - - const Real nu = kv/(k_norm*c); - const Complex theta = amrex::exp( 0.5_rt*I*kv*dt ); - const Complex theta_star = amrex::exp( -0.5_rt*I*kv*dt ); - const Complex e_theta = amrex::exp( I*c*k_norm*dt ); - - Theta2(i,j,k) = theta*theta; - - if ( (nu != 1.) && (nu != 0) ) { - - // Note: the coefficients X1, X2, X3 do not correspond - // exactly to the original Galilean paper, but the - // update equation have been modified accordingly so that - // the expressions/ below (with the update equations) - // are mathematically equivalent to those of the paper. - Complex x1 = 1._rt/(1._rt-nu*nu) * - (theta_star - C(i,j,k)*theta + I*kv*S_ck(i,j,k)*theta); - - Complex C_rho = I* c2 /( (1._rt-theta*theta) * ep0); - - Psi1(i,j,k) = theta * ((S1(i,j,k) + I*nu*C1(i,j,k)) - - Theta2(i,j,k) * (S3(i,j,k) + I*nu*C3(i,j,k))) /(c*k_norm*dt * (nu*nu - 1._rt)); - Psi2(i,j,k) = theta * ((C1(i,j,k) - I*nu*S1(i,j,k)) - - Theta2(i,j,k) * (C3(i,j,k) - I*nu*S3(i,j,k))) /(c2*k_norm*k_norm*dt * (nu*nu - 1._rt)); - Psi3(i,j,k) = I * theta * (1._rt - theta*theta) /(c*k_norm*dt*nu); - - A1(i,j,k) = (Psi1(i,j,k) - 1._rt + I * kv*Psi2(i,j,k) )/ (c2* k_norm*k_norm * (nu*nu - 1._rt)); - A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / (c2*k_norm*k_norm); - - CRhoold(i,j,k) = C_rho * (theta*theta * A1(i,j,k) - A2(i,j,k)); - CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k)); - Jcoef(i,j,k) = (I*kv*A1(i,j,k) + Psi2(i,j,k))/ep0; - // x1, above, is identical to the original paper - X1(i,j,k) = theta*x1/(ep0*c*c*k_norm*k_norm); - // The difference betwen X2 and X3 below, and those - // from the original paper is the factor ep0*k_norm*k_norm - X2(i,j,k) = (x1 - theta*(1._rt - C(i,j,k))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X3(i,j,k) = (x1 - theta_star*(1._rt - C(i,j,k))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X4(i,j,k) = I*kv*X1(i,j,k) - theta*theta*S_ck(i,j,k)/ep0; - } - if ( nu == 0) { - X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0*c*c*k_norm*k_norm); - X2(i,j,k) = (1._rt - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm); - X4(i,j,k) = -S_ck(i,j,k)/ep0; - - Psi1(i,j,k) = (-S1(i,j,k) + S3(i,j,k)) / (c*k_norm*dt); - Psi2(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*k_norm*k_norm*dt); - Psi3(i,j,k) = 1._rt; - A1(i,j,k) = (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt); - A2(i,j,k) = (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt); - CRhoold(i,j,k) = 2._rt * I * S1(i,j,k) * ( dt*C(i,j,k) - S_ck(i,j,k)) - / (c*k_norm*k_norm*k_norm*dt*dt*ep0); - CRhonew(i,j,k) = - I * (c2* k_norm*k_norm * dt*dt - C1(i,j,k) + C3(i,j,k)) - / (c2 * k_norm*k_norm*k_norm*k_norm * ep0 * dt*dt); - Jcoef(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*ep0*k_norm*k_norm*dt); - } - if ( nu == 1.) { - X1(i,j,k) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm); - X2(i,j,k) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt- e_theta)); - X3(i,j,k) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm); - X4(i,j,k) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*c*k_norm); - } - - } else { // Handle k_norm = 0, by using the analytical limit - C(i,j,k) = 1._rt; - S_ck(i,j,k) = dt; - C1(i,j,k) = 1._rt; - S1(i,j,k) = 0._rt; - C3(i,j,k) = 1._rt; - S3(i,j,k) = 0._rt; - - X1(i,j,k) = dt*dt/(2._rt * ep0); - X2(i,j,k) = c2*dt*dt/(6._rt * ep0); - X3(i,j,k) = - c2*dt*dt/(3._rt * ep0); - X4(i,j,k) = -dt/ep0; - Theta2(i,j,k) = 1._rt; - - Psi1(i,j,k) = 1._rt; - Psi2(i,j,k) = -dt; - Psi3(i,j,k) = 1._rt; - A1(i,j,k) = 13._rt * dt*dt /24._rt; - A2(i,j,k) = 13._rt * dt*dt /24._rt; - CRhoold(i,j,k) = -I*c2 * dt*dt / (3._rt * ep0); - CRhonew(i,j,k) = -5._rt*I*c2 * dt*dt / (24._rt * ep0); - Jcoef(i,j,k) = -dt/ep0; - } - - }); - } -} - -/* Advance the E and B field in spectral space (stored in `f`) - * over one time step */ -void -AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ - - // Loop over boxes - for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ - - const Box& bx = f.fields[mfi].box(); - - // Extract arrays for the fields to be updated - Array4 fields = f.fields[mfi].array(); - // Extract arrays for the coefficients - Array4 C_arr = C_coef[mfi].array(); - Array4 S_ck_arr = S_ck_coef[mfi].array(); - Array4 X1_arr = X1_coef[mfi].array(); - Array4 X2_arr = X2_coef[mfi].array(); - Array4 X3_arr = X3_coef[mfi].array(); - Array4 X4_arr = X4_coef[mfi].array(); - Array4 Theta2_arr = Theta2_coef[mfi].array(); - Array4 Psi1_arr = Psi1_coef[mfi].array(); - Array4 Psi2_arr = Psi2_coef[mfi].array(); - - Array4 A1_arr = A1_coef[mfi].array(); - Array4 Rhonew_arr = Rhonew_coef[mfi].array(); - Array4 Rhoold_arr = Rhoold_coef[mfi].array(); - Array4 Jcoef_arr =Jcoef_coef[mfi].array(); - // Extract pointers for the k vectors - const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); -#if (AMREX_SPACEDIM==3) - const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); - - // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Record old values of the fields to be updated - using Idx = SpectralAvgFieldIndex; - - const Complex Ex_old = fields(i,j,k,Idx::Ex); - const Complex Ey_old = fields(i,j,k,Idx::Ey); - const Complex Ez_old = fields(i,j,k,Idx::Ez); - const Complex Bx_old = fields(i,j,k,Idx::Bx); - const Complex By_old = fields(i,j,k,Idx::By); - const Complex Bz_old = fields(i,j,k,Idx::Bz); - - // Shortcut for the values of J and rho - const Complex Jx = fields(i,j,k,Idx::Jx); - const Complex Jy = fields(i,j,k,Idx::Jy); - const Complex Jz = fields(i,j,k,Idx::Jz); - const Complex rho_old = fields(i,j,k,Idx::rho_old); - const Complex rho_new = fields(i,j,k,Idx::rho_new); - - // k vector values, and coefficients - const Real kx = modified_kx_arr[i]; - -#if (AMREX_SPACEDIM==3) - 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]; -#endif - constexpr Real c2 = PhysConst::c*PhysConst::c; - constexpr Real inv_ep0 = 1._rt/PhysConst::ep0; - constexpr Complex I = Complex{0,1}; - - const Real C = C_arr(i,j,k); - const Real S_ck = S_ck_arr(i,j,k); - - const Complex X1 = X1_arr(i,j,k); - const Complex X2 = X2_arr(i,j,k); - const Complex X3 = X3_arr(i,j,k); - const Complex X4 = X4_arr(i,j,k); - const Complex T2 = Theta2_arr(i,j,k); - - const Complex Psi1 = Psi1_arr(i,j,k); - const Complex Psi2 = Psi2_arr(i,j,k); - const Complex A1 = A1_arr(i,j,k); - const Complex CRhoold= Rhoold_arr(i,j,k); - const Complex CRhonew= Rhonew_arr(i,j,k); - const Complex Jcoef = Jcoef_arr(i,j,k); - - //Update E (see the original Galilean article) - fields(i,j,k,Idx::Ex) = T2*C*Ex_old - + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old) - + X4*Jx - I*(X2*rho_new - T2*X3*rho_old)*kx; - fields(i,j,k,Idx::Ey) = T2*C*Ey_old - + T2*S_ck*c2*I*(kz*Bx_old - kx*Bz_old) - + X4*Jy - I*(X2*rho_new - T2*X3*rho_old)*ky; - fields(i,j,k,Idx::Ez) = T2*C*Ez_old - + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old) - + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz; - - // Update B (see the original Galilean article) - // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where - // x1 has the same definition as in the original paper - fields(i,j,k,Idx::Bx) = T2*C*Bx_old - - T2*S_ck*I*(ky*Ez_old - kz*Ey_old) - + X1*I*(ky*Jz - kz*Jy); - fields(i,j,k,Idx::By) = T2*C*By_old - - T2*S_ck*I*(kz*Ex_old - kx*Ez_old) - + X1*I*(kz*Jx - kx*Jz); - fields(i,j,k,Idx::Bz) = T2*C*Bz_old - - T2*S_ck*I*(kx*Ey_old - ky*Ex_old) - + X1*I*(kx*Jy - ky*Jx); - - //Update the averaged E,B fields in time on the interval [(n-1/2)dx, (n+1/2)dx] - fields(i,j,k,Idx::Ex_avg) = Psi1*Ex_old - - Psi2*c2*I*(ky*Bz_old - kz*By_old) - + Jcoef*Jx + ( CRhonew * rho_new + CRhoold*rho_old )*kx; - fields(i,j,k,Idx::Ey_avg) = Psi1*Ey_old - - Psi2*c2*I*(kz*Bx_old - kx*Bz_old) - + Jcoef*Jy +( CRhonew * rho_new + CRhoold*rho_old )*ky; - fields(i,j,k,Idx::Ez_avg) = Psi1*Ez_old - - Psi2*c2*I*(kx*By_old - ky*Bx_old) - + Jcoef*Jz + ( CRhonew * rho_new + CRhoold*rho_old )*kz; - - fields(i,j,k,Idx::Bx_avg) = Psi1*Bx_old - + I*Psi2*(ky*Ez_old - kz*Ey_old) - + A1*I*(ky*Jz - kz*Jy)*inv_ep0; - fields(i,j,k,Idx::By_avg) = Psi1*By_old - + I*Psi2*(kz*Ex_old - kx*Ez_old) - + A1*I*(kz*Jx - kx*Jz)*inv_ep0; - fields(i,j,k,Idx::Bz_avg) = Psi1*Bz_old - + I*Psi2*(kx*Ey_old - ky*Ex_old) - + A1*I*(kx*Jy - ky*Jx)*inv_ep0; - }); - } -} - -void -AvgGalileanAlgorithm::CurrentCorrection (SpectralFieldData& /*field_data*/, - std::array,3>& /*current*/, - const std::unique_ptr& /*rho*/) -{ - amrex::Abort("Current correction not implemented for averaged Galilean PSATD"); -} - -void -AvgGalileanAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/, - std::array,3>& /*current*/) -{ - amrex::Abort("Vay deposition not implemented for averaged Galilean PSATD"); -} -- cgit v1.2.3