aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/FieldSolver/WarpX_QED_Field_Pushers.cpp36
-rw-r--r--Source/FieldSolver/WarpX_QED_K.H33
2 files changed, 52 insertions, 17 deletions
diff --git a/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp b/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp
index 4892610c6..4aece5178 100644
--- a/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp
+++ b/Source/FieldSolver/WarpX_QED_Field_Pushers.cpp
@@ -1,4 +1,4 @@
-/* Copyright 2019-2020 Glenn Richardson, Maxence Thevenet
+/* Copyright 2019-2020 Glenn Richardson, Maxence Thevenet, Revathi Jambunathan, Axel Huebl
*
* This file is part of WarpX.
*
@@ -10,6 +10,7 @@
#include "BoundaryConditions/WarpX_PML_kernels.H"
#include "BoundaryConditions/PML_current.H"
#include "WarpX_FDTD.H"
+
#ifdef WARPX_USE_PY
# include "Python/WarpX_py.H"
#endif
@@ -63,7 +64,7 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
const Real dy = dx_vec[1];
const Real dz = dx_vec[2];
- MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
+ MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *Jx, *Jy, *Jz;
if (patch_type == PatchType::fine)
{
Ex = Efield_fp[lev][0].get();
@@ -72,6 +73,9 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
Bx = Bfield_fp[lev][0].get();
By = Bfield_fp[lev][1].get();
Bz = Bfield_fp[lev][2].get();
+ Jx = current_fp[lev][0].get();
+ Jy = current_fp[lev][1].get();
+ Jz = current_fp[lev][2].get();
}
else
{
@@ -81,6 +85,9 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
Bx = Bfield_cp[lev][0].get();
By = Bfield_cp[lev][1].get();
Bz = Bfield_cp[lev][2].get();
+ Jx = current_cp[lev][0].get();
+ Jy = current_cp[lev][1].get();
+ Jz = current_cp[lev][2].get();
}
amrex::Vector<amrex::Real>* cost = WarpX::getCosts(lev);
@@ -97,12 +104,19 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
}
Real wt = amrex::second();
- // Get boxes for E and B
- const Box& tbx = mfi.tilebox(Bx_nodal_flag);
+ // Get boxes for E, B, and J
+
+ const Box& tbx = mfi.tilebox(Bx->ixType().ixType());
+ const Box& tby = mfi.tilebox(By->ixType().ixType());
+ const Box& tbz = mfi.tilebox(Bz->ixType().ixType());
- const Box& tex = mfi.tilebox(Ex_nodal_flag);
- const Box& tey = mfi.tilebox(Ey_nodal_flag);
- const Box& tez = mfi.tilebox(Ez_nodal_flag);
+ const Box& tex = mfi.tilebox(Ex->ixType().ixType());
+ const Box& tey = mfi.tilebox(Ey->ixType().ixType());
+ const Box& tez = mfi.tilebox(Ez->ixType().ixType());
+
+ const Box& tjx = mfi.tilebox(Jx->ixType().ixType());
+ const Box& tjy = mfi.tilebox(Jy->ixType().ixType());
+ const Box& tjz = mfi.tilebox(Jz->ixType().ixType());
// Get field arrays
auto const& Bxfab = Bx->array(mfi);
@@ -111,6 +125,9 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
auto const& Exfab = Ex->array(mfi);
auto const& Eyfab = Ey->array(mfi);
auto const& Ezfab = Ez->array(mfi);
+ auto const& Jxfab = Jx->array(mfi);
+ auto const& Jyfab = Jy->array(mfi);
+ auto const& Jzfab = Jz->array(mfi);
// Define grown box with 1 ghost cell for finite difference stencil
const Box& gex = amrex::grow(tex,1);
@@ -148,14 +165,15 @@ WarpX::Hybrid_QED_Push (int lev, PatchType patch_type, Real a_dt)
// Make local copy of xi, to use on device.
const Real xi_c2 = WarpX::quantum_xi_c2;
+
// Apply QED correction to electric field, using temporary arrays.
amrex::ParallelFor(
tbx,
[=] AMREX_GPU_DEVICE (int j, int k, int l)
{
warpx_hybrid_QED_push(j,k,l, Exfab, Eyfab, Ezfab, Bxfab, Byfab,
- Bzfab, tmpEx, tmpEy, tmpEz, dx, dy, dz,
- a_dt, xi_c2);
+ Bzfab, tmpEx, tmpEy, tmpEz, Jxfab, Jyfab, Jzfab, dx, dy, dz,
+ a_dt, xi_c2);
}
);
diff --git a/Source/FieldSolver/WarpX_QED_K.H b/Source/FieldSolver/WarpX_QED_K.H
index d030280b0..e71f40cfc 100644
--- a/Source/FieldSolver/WarpX_QED_K.H
+++ b/Source/FieldSolver/WarpX_QED_K.H
@@ -74,8 +74,9 @@ void warpx_hybrid_QED_push (
amrex::Array4<amrex::Real> const& Bx, amrex::Array4<amrex::Real> const& By,
amrex::Array4<amrex::Real const> const& Bz, amrex::Array4<amrex::Real> const& tmpEx,
amrex::Array4<amrex::Real> const& tmpEy, amrex::Array4<amrex::Real> const& tmpEz,
- amrex::Real dx, amrex::Real dy, amrex::Real dz, amrex::Real dt,
- const amrex::Real xi_c2)
+ amrex::Array4<amrex::Real> const& Jx, amrex::Array4<amrex::Real> const& Jy,
+ amrex::Array4<amrex::Real> const& Jz, const amrex::Real dx, const amrex::Real dy,
+ const amrex::Real dz, const amrex::Real dt, const amrex::Real xi_c2)
{
using namespace amrex;
@@ -142,6 +143,9 @@ const amrex::Real dyi = 1./dy;
const amrex::Real bx = Bx(j,k,l);
const amrex::Real by = By(j,k,l);
const amrex::Real bz = Bz(j,k,l);
+ const amrex::Real mu0jx = PhysConst::mu0*Jx(j,k,l);
+ const amrex::Real mu0jy = PhysConst::mu0*Jy(j,k,l);
+ const amrex::Real mu0jz = PhysConst::mu0*Jz(j,k,l);
const amrex::Real ee = ex*ex + ey*ey + ez*ez;
const amrex::Real bb = bx*bx + by*by + bz*bz;
const amrex::Real eb = ex*bx + ey*by + ez*bz;
@@ -149,6 +153,8 @@ const amrex::Real dyi = 1./dy;
const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2];
const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2];
const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2];
+ const amrex::Real Emu0J = ex*mu0jx + ey*mu0jy + ez*mu0jz;
+ const amrex::Real Bmu0J = bx*mu0jx + by*mu0jy + bz*mu0jz;
const amrex::Real beta = 4._rt*xi_c2*( c2i*ee - bb ) + PhysConst::ep0;
@@ -159,9 +165,12 @@ const amrex::Real dyi = 1./dy;
};
const amrex::Real Omega[3] = {
- Alpha[0] + 2._rt*xi_c2*( 4._rt*ex*EVxB + 2._rt*VxB[0]*( ee - c2*bb ) + 7._rt*c2*bx*BVxB ),
- Alpha[1] + 2._rt*xi_c2*( 4._rt*ey*EVxB + 2._rt*VxB[1]*( ee - c2*bb ) + 7._rt*c2*by*BVxB ),
- Alpha[2] + 2._rt*xi_c2*( 4._rt*ez*EVxB + 2._rt*VxB[2]*( ee - c2*bb ) + 7._rt*c2*bz*BVxB )
+ Alpha[0] + 2._rt*xi_c2*( 4._rt*ex*(EVxB + Emu0J) + 2._rt*(VxB[0] + mu0jx)*( ee - c2*bb )
+ + 7._rt*c2*bx*(BVxB +Bmu0J) ),
+ Alpha[1] + 2._rt*xi_c2*( 4._rt*ey*(EVxB + Emu0J) + 2._rt*(VxB[1] + mu0jy)*( ee - c2*bb )
+ + 7._rt*c2*by*(BVxB + Bmu0J) ),
+ Alpha[2] + 2._rt*xi_c2*( 4._rt*ez*(EVxB + Emu0J) + 2._rt*(VxB[2] + mu0jz)*( ee - c2*bb )
+ + 7._rt*c2*bz*(BVxB + Bmu0J) )
};
// Calcualting matrix values for the QED correction algorithm
@@ -260,6 +269,9 @@ const amrex::Real dyi = 1./dy;
const amrex::Real bx = Bx(j,k,0);
const amrex::Real by = By(j,k,0);
const amrex::Real bz = Bz(j,k,0);
+ const amrex::Real mu0jx = PhysConst::mu0*Jx(j,k,0);
+ const amrex::Real mu0jy = PhysConst::mu0*Jy(j,k,0);
+ const amrex::Real mu0jz = PhysConst::mu0*Jz(j,k,0);
const amrex::Real ee = ex*ex + ey*ey + ez*ez;
const amrex::Real bb = bx*bx + by*by + bz*bz;
const amrex::Real eb = ex*bx + ey*by + ez*bz;
@@ -267,6 +279,8 @@ const amrex::Real dyi = 1./dy;
const amrex::Real BVxE = bx*VxE[0] + by*VxE[1] + bz*VxE[2];
const amrex::Real EVxB = ex*VxB[0] + ey*VxB[1] + ez*VxB[2];
const amrex::Real BVxB = bx*VxB[0] + by*VxB[1] + bz*VxB[2];
+ const amrex::Real Emu0J = ex*mu0jx + ey*mu0jy + ez*mu0jz;
+ const amrex::Real Bmu0J = bx*mu0jx + by*mu0jy + bz*mu0jz;
const amrex::Real beta = 4._rt*xi_c2*( c2i*ee - bb ) + PhysConst::ep0;
@@ -277,9 +291,12 @@ const amrex::Real dyi = 1./dy;
};
const amrex::Real Omega[3] = {
- Alpha[0] + 2._rt*xi_c2*( 4._rt*ex*EVxB + 2._rt*VxB[0]*( ee - c2*bb ) + 7._rt*c2*bx*BVxB ),
- Alpha[1] + 2._rt*xi_c2*( 4._rt*ey*EVxB + 2._rt*VxB[1]*( ee - c2*bb ) + 7._rt*c2*by*BVxB ),
- Alpha[2] + 2._rt*xi_c2*( 4._rt*ez*EVxB + 2._rt*VxB[2]*( ee - c2*bb ) + 7._rt*c2*bz*BVxB )
+ Alpha[0] + 2._rt*xi_c2*( 4._rt*ex*(EVxB + Emu0J) + 2._rt*(VxB[0] + mu0jx)*( ee - c2*bb )
+ + 7._rt*c2*bx*(BVxB + Bmu0J) ),
+ Alpha[1] + 2._rt*xi_c2*( 4._rt*ey*(EVxB + Emu0J) + 2._rt*(VxB[1] + mu0jy)*( ee - c2*bb )
+ + 7._rt*c2*by*(BVxB + Bmu0J) ),
+ Alpha[2] + 2._rt*xi_c2*( 4._rt*ez*(EVxB +Emu0J) + 2._rt*(VxB[2] + mu0jz)*( ee - c2*bb )
+ + 7._rt*c2*bz*(BVxB + Bmu0J) )
};
// Calcualting matrix values for the QED correction algorithm
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-12 13:22:25 +0000