diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 28b5344347..2905bdb929 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -765,7 +765,7 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st m_multi_plasma.ExplicitDeposition(m_fields, m_3D_geom, lev); // Solves Bx, By using Sx, Sy and chi - ExplicitMGSolveBxBy(lev, WhichSlice::This); + // ExplicitMGSolveBxBy(lev, WhichSlice::This); } } else { // Solves Bx and By in the current slice and modifies the force terms of the plasma particles diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 46ec1b8106..36af79879a 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -964,17 +964,23 @@ Fields::SolvePoissonPsiExmByEypBxEzBz (amrex::Vector const& geo m_poisson_solver[lev]->SolvePoissonEquation(lhs_Psi); // Ez: right-hand side 1/(episilon0 *c0 )*(d_x(jx) + d_y(jy)) +// LinCombination(getStagingArea(lev), +// 1._rt/(phys_const.ep0*phys_const.c), +// derivative{getField(lev, WhichSlice::This, "jx"), geom[lev]}, +// 1._rt/(phys_const.ep0*phys_const.c), +// derivative{getField(lev, WhichSlice::This, "jy"), geom[lev]}); LinCombination(getStagingArea(lev), 1._rt/(phys_const.ep0*phys_const.c), - derivative{getField(lev, WhichSlice::This, "jx"), geom[lev]}, + getField(lev, WhichSlice::This, "jx"), 1._rt/(phys_const.ep0*phys_const.c), - derivative{getField(lev, WhichSlice::This, "jy"), geom[lev]}); + getField(lev, WhichSlice::This, "jy")); SetBoundaryCondition(geom, lev, WhichSlice::This, "Ez", getStagingArea(lev), m_poisson_solver[lev]->BoundaryOffset(), m_poisson_solver[lev]->BoundaryFactor()); m_poisson_solver[lev]->SolvePoissonEquation(lhs_Ez); +/* // Bz: right-hand side mu_0*(d_y(jx) - d_x(jy)) LinCombination(getStagingArea(lev), phys_const.mu0, @@ -986,6 +992,7 @@ Fields::SolvePoissonPsiExmByEypBxEzBz (amrex::Vector const& geo m_poisson_solver[lev]->BoundaryOffset(), m_poisson_solver[lev]->BoundaryFactor()); m_poisson_solver[lev]->SolvePoissonEquation(lhs_Bz); +*/ } EnforcePeriodic(false, {Comps[WhichSlice::This]["Psi"], diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 5a2cbb6b22..b1fce414ce 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/src/particles/deposition/PlasmaDepositCurrent.cpp @@ -217,10 +217,13 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, // y direction auto [shape_y, j] = shape_factor(ymid, iy); + auto [shape2_y, shape_dy, j2] = derivative_shape_factor<0, depos_order>(ymid, iy); + auto [shape2_x, shape_dx, i2] = derivative_shape_factor<0, depos_order>(xmid, ix); + const amrex::Real charge_density = q_invvol * shape_x * shape_y; // wqx, wqy wqz are particle current in each direction - const amrex::Real wqx = charge_density * clight * vx; - const amrex::Real wqy = charge_density * clight * vy; + const amrex::Real wqx = q_invvol * shape_dx * shape2_y * clight * vx * dx_inv; + const amrex::Real wqy = q_invvol * shape2_x * shape_dy * clight * vy * dy_inv; const amrex::Real wqz = charge_density * clight * (gamma_psi-1._rt); const amrex::Real wq = charge_density * gamma_psi; const amrex::Real wchi = charge_density * q_mu0_mass_ratio * psi_inv; @@ -228,8 +231,10 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, // Deposit current into arr if (depos_idx[0] != -1) { // deposit_jx_jy - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), wqx); - amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[1]), wqy); +// amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[0]), wqx); +// amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[1]), wqy); + amrex::Gpu::Atomic::Add(arr.ptr(i2, j2, depos_idx[0]), wqx); + amrex::Gpu::Atomic::Add(arr.ptr(i2, j2, depos_idx[1]), wqy); } if (depos_idx[2] != -1) { // deposit_jz amrex::Gpu::Atomic::Add(arr.ptr(i, j, depos_idx[2]), wqz); diff --git a/src/particles/particles_utils/FieldGather.H b/src/particles/particles_utils/FieldGather.H index 40d9c39a9a..cdcf8ade3c 100644 --- a/src/particles/particles_utils/FieldGather.H +++ b/src/particles/particles_utils/FieldGather.H @@ -85,7 +85,7 @@ void doGatherShapeN (const amrex::Real xp, } } - constexpr int derivative_type = 1; + constexpr int derivative_type = 0; HIPACE_LOOP_UNROLL for (int iy=0; iy<=depos_order_xy+derivative_type; iy++){