From 98e808ff1579603b1467debddd61c9b31468c61b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 14 Feb 2022 17:15:39 -0800 Subject: Allow flux injection in the out-of-plane direction for RZ/2D geometry (#2788) * Implement injection orthogal to plane * Generalize momentum distribution for flux injection * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Revert "[pre-commit.ci] auto fixes from pre-commit.com hooks" This reverts commit b0cd1891771a4c49c14abb7cb9df7374cee4458c. * Revert "Generalize momentum distribution for flux injection" This reverts commit 0a22b1d8fa68a3a5705d8f4824f757b6dee497f0. * Rotate momentum initialization * Correct flux number when the direction is normal to plane * Update distribution of particles within a cell * Clean-up injection code * Add more documentation * Add more comments * Handle 1D case * Only do the rotation for Gaussian flux profile * Fix compilation error * Correct compilation for GPU * Start adding automated test * Correct sign of velocity * Update to add continuous injection * Finalize test * Correct processing of flux_normal_axis * Add checksum * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix bug * Update script * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update checksum Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- Source/Particles/PhysicalParticleContainer.cpp | 73 +++++++++++++++++++++----- 1 file changed, 61 insertions(+), 12 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0e81d4000..27d94a5f1 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1235,9 +1235,23 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) const auto problo = geom.ProbLoArray(); Real scale_fac = 0._rt; - if (dx[plasma_injector->flux_normal_axis]*num_ppc_real != 0._rt) { - scale_fac = AMREX_D_TERM(dx[0], *dx[1], *dx[2])/dx[plasma_injector->flux_normal_axis]/num_ppc_real; - } + // Scale particle weight by the area of the emitting surface, within one cell +#if defined(WARPX_DIM_3D) + scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector->flux_normal_axis]/num_ppc_real; +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + scale_fac = dx[0]*dx[1]/num_ppc_real; + // When emission is in the r direction, the emitting surface is a cylinder. + // The factor 2*pi*r is added later below. + if (plasma_injector->flux_normal_axis == 0) scale_fac /= dx[0]; + // When emission is in the z direction, the emitting surface is an annulus + // The factor 2*pi*r is added later below. + if (plasma_injector->flux_normal_axis == 2) scale_fac /= dx[1]; + // When emission is in the theta direction (flux_normal_axis == 1), + // the emitting surface is a rectangle, within the plane of the simulation +#elif defined(WARPX_DIM_1D_Z) + scale_fac = dx[0]/num_ppc_real; + if (plasma_injector->flux_normal_axis == 2) scale_fac /= dx[0]; +#endif amrex::LayoutData* cost = WarpX::getCosts(0); @@ -1301,7 +1315,16 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) bool no_overlap = false; for (int dir=0; dirflux_normal_axis) { +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + if (2*dir == plasma_injector->flux_normal_axis) { + // The above formula captures the following cases: + // - flux_normal_axis=0 (emission along x/r) and dir=0 + // - flux_normal_axis=2 (emission along z) and dir=1 +#elif defined(WARPX_DIM_1D_Z) + if ( (dir==0) && (plasma_injector->flux_normal_axis==2) ) { +#endif if (plasma_injector->flux_direction > 0) { if (plasma_injector->surface_flux_pos < tile_realbox.lo(dir) || plasma_injector->surface_flux_pos >= tile_realbox.hi(dir)) { @@ -1459,6 +1482,9 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; +#ifdef WARPX_DIM_RZ + int const loc_flux_normal_axis = plasma_injector->flux_normal_axis; +#endif // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). @@ -1528,8 +1554,23 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) } else { theta = 2._rt*MathConst::pi*r.y; } - pos.x = xb*std::cos(theta); - pos.y = xb*std::sin(theta); + Real const cos_theta = std::cos(theta); + Real const sin_theta = std::sin(theta); + // Rotate the position + pos.x = xb*cos_theta; + pos.y = xb*sin_theta; + if ( (inj_mom->type == InjectorMomentum::Type::gaussianflux) && + (loc_flux_normal_axis != 2) ) { + // Rotate the momentum + // This because, when the flux direction is e.g. "r" + // the `inj_mom` objects generates a v*Gaussian distribution + // along the Cartesian "x" directionm by default. This + // needs to be rotated along "r". + Real ur = u.x; + Real ut = u.y; + u.x = cos_theta*ur - sin_theta*ut; + u.y = sin_theta*ur + cos_theta*ut; + } #endif // Lab-frame simulation @@ -1569,13 +1610,21 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) // Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); Real weight = dens * scale_fac * dt; #ifdef WARPX_DIM_RZ - if (radially_weighted) { - weight *= 2._rt*MathConst::pi*xb; - } else { - // This is not correct since it might shift the particle - // out of the local grid - pos.x = std::sqrt(xb*rmax); - weight *= dx[0]; + // The particle weight is proportional to the user-specified + // flux (denoted as `dens` here) and the emission surface within + // one cell (captured partially by `scale_fac`). + // For cylindrical emission (flux_normal_axis==0 + // or flux_normal_axis==2), the emission surface depends on + // the radius ; thus, the calculation is finalized here + if (loc_flux_normal_axis != 1) { + if (radially_weighted) { + weight *= 2._rt*MathConst::pi*xb; + } else { + // This is not correct since it might shift the particle + // out of the local grid + pos.x = std::sqrt(xb*rmax); + weight *= dx[0]; + } } #endif pa[PIdx::w ][ip] = weight; -- cgit v1.2.3