From 107425cbb59bfa13b6e1de56881bd5f85dd8fd77 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 5 Jul 2019 16:47:14 -0700 Subject: Converted E push to c++ --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 53 ++++++++++++++++---------------- 1 file changed, 27 insertions(+), 26 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index c53e13f8f..9c5c3f834 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #ifdef WARPX_USE_PY #include #endif @@ -192,6 +193,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; + Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; if (patch_type == PatchType::fine) @@ -240,16 +242,17 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& jxfab = jx->array(mfi); + auto const& jyfab = jy->array(mfi); + auto const& jzfab = jz->array(mfi); + if (do_nodal) { - auto const& Exfab = Ex->array(mfi); - auto const& Eyfab = Ey->array(mfi); - auto const& Ezfab = Ez->array(mfi); - auto const& Bxfab = Bx->array(mfi); - auto const& Byfab = By->array(mfi); - auto const& Bzfab = Bz->array(mfi); - auto const& jxfab = jx->array(mfi); - auto const& jyfab = jy->array(mfi); - auto const& jzfab = jz->array(mfi); amrex::ParallelFor(tex, [=] AMREX_GPU_DEVICE (int j, int k, int l) { @@ -266,23 +269,21 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { - // Call picsar routine for each tile - warpx_push_evec( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - BL_TO_FORTRAN_3D((*jx)[mfi]), - BL_TO_FORTRAN_3D((*jy)[mfi]), - BL_TO_FORTRAN_3D((*jz)[mfi]), - &mu_c2_dt, - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, - &xmin, &dx[0]); + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); + }); } if (F) -- cgit v1.2.3 From 97165ec52b7f353c3392c66beb1e70bb3cfdfad4 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 5 Jul 2019 17:28:17 -0700 Subject: Converted Yee B push to c++ --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 42 +++++++++++----------- Source/FieldSolver/WarpX_FDTD.H | 60 ++++++++++++++++++++++++++++---- 2 files changed, 76 insertions(+), 26 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 9c5c3f834..6d5b38d04 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -42,6 +42,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; + Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -82,13 +83,13 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); if (do_nodal) { - auto const& Bxfab = Bx->array(mfi); - auto const& Byfab = By->array(mfi); - auto const& Bzfab = Bz->array(mfi); - auto const& Exfab = Ex->array(mfi); - auto const& Eyfab = Ey->array(mfi); - auto const& Ezfab = Ez->array(mfi); amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int j, int k, int l) { @@ -105,20 +106,21 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); } else { - // Call picsar routine for each tile - warpx_push_bvec( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - &dtsdx, &dtsdy, &dtsdz, - &xmin, &dx[0], - &WarpX::maxwell_fdtd_solver_id); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); + }); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); + }); } if (cost) { diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index c5fef2149..7092fea31 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -5,10 +5,58 @@ using namespace amrex; +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bx_yee(int i, int j, int k, Array4 const& Bx, + Array4 const& Ey, Array4 const& Ez, + Real dtsdy_c2, Real dtsdz_c2) +{ +#if (AMREX_SPACEDIM == 3) + Bx(i,j,k) = Bx(i,j,k) - dtsdy_c2 * (Ez(i,j+1,k ) - Ez(i,j,k)) + + dtsdz_c2 * (Ey(i,j ,k+1) - Ey(i,j,k)); +#else + // Note that the 2D Cartesian and RZ versions are the same + Bx(i,j,0) = Bx(i,j,0) + dtsdz_c2 * (Ey(i,j+1,0) - Ey(i,j,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_by_yee(int i, int j, int k, Array4 const& By, + Array4 const& Ex, Array4 const& Ez, + Real dtsdx_c2, Real dtsdz_c2) +{ +#if (AMREX_SPACEDIM == 3) + By(i,j,k) = By(i,j,k) + dtsdx_c2 * (Ez(i+1,j,k ) - Ez(i,j,k)) + - dtsdz_c2 * (Ex(i ,j,k+1) - Ex(i,j,k)); +#else + // Note that the 2D Cartesian and RZ versions are the same + By(i,j,0) = By(i,j,0) + dtsdx_c2 * (Ez(i+1,j ,0) - Ez(i,j,0)) + - dtsdz_c2 * (Ex(i ,j+1,0) - Ex(i,j,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bz_yee(int i, int j, int k, Array4 const& Bz, + Array4 const& Ex, Array4 const& Ey, + Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin) +{ +#if (AMREX_SPACEDIM == 3) + Bz(i,j,k) = Bz(i,j,k) - dtsdx_c2 * (Ey(i+1,j ,k) - Ey(i,j,k)) + + dtsdy_c2 * (Ex(i ,j+1,k) - Ex(i,j,k)); +#else +#ifndef WARPX_RZ + Bz(i,j,0) = Bz(i,j,0) - dtsdx_c2 * (Ey(i+1,j,0) - Ey(i,j,0)); +#else + Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); + Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); + Bz(i,j,0) = Bz(i,j,0) - dtsdx_c2 * (ru*Ey(i+1,j,0) - rd*Ey(i,j,0)); +#endif +#endif +} + AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ex_yee(int i, int j, int k, Array4 const& Ex, - Array4 const& By, Array4 const& Bz, Array4 const& Jx, - Real mu_c2_dt, Real dtsdy_c2, Real dtsdz_c2) + Array4 const& By, Array4 const& Bz, Array4 const& Jx, + Real mu_c2_dt, Real dtsdy_c2, Real dtsdz_c2) { #if (AMREX_SPACEDIM == 3) Ex(i,j,k) = Ex(i,j,k) + dtsdy_c2 * (Bz(i,j,k) - Bz(i,j-1,k )) @@ -23,8 +71,8 @@ void warpx_push_ex_yee(int i, int j, int k, Array4 const& Ex, AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ey_yee(int i, int j, int k, Array4 const& Ey, - Array4 const& Bx, Array4 const& Bz, Array4 const& Jy, - Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin) + Array4 const& Bx, Array4 const& Bz, Array4 const& Jy, + Real mu_c2_dt, Real dtsdx_c2, Real dtsdz_c2, Real rmin) { #if (AMREX_SPACEDIM == 3) Ey(i,j,k) = Ey(i,j,k) - dtsdx_c2 * (Bz(i,j,k) - Bz(i-1,j,k)) @@ -45,8 +93,8 @@ void warpx_push_ey_yee(int i, int j, int k, Array4 const& Ey, AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ez_yee(int i, int j, int k, Array4 const& Ez, - Array4 const& Bx, Array4 const& By, Array4 const& Jz, - Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin) + Array4 const& Bx, Array4 const& By, Array4 const& Jz, + Real mu_c2_dt, Real dtsdx_c2, Real dtsdy_c2, Real dxinv, Real rmin) { #if (AMREX_SPACEDIM == 3) Ez(i,j,k) = Ez(i,j,k) + dtsdx_c2 * (By(i,j,k) - By(i-1,j ,k)) -- cgit v1.2.3 From edb6e6a9fb08dcba169c395f5338deaf373356c7 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 5 Jul 2019 23:02:58 -0700 Subject: Converted CKC routines to c++ --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 34 ++++++- Source/FieldSolver/WarpX_FDTD.H | 157 +++++++++++++++++++++++++++++++ Source/FortranInterface/WarpX_picsar.F90 | 74 --------------- 3 files changed, 190 insertions(+), 75 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6d5b38d04..06b6b1b03 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -105,7 +105,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); - } else { + } else if (WarpX::maxwell_fdtd_solver_id == 0) { amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int j, int k, int l) { @@ -121,6 +121,38 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); }); + } else if (WarpX::maxwell_fdtd_solver_id == 1) { + Real betaxy, betaxz, betayx, betayz, betazx, betazy; + Real gammax, gammay, gammaz; + Real alphax, alphay, alphaz; + warpx_calculate_ckc_coefficients(dtsdx, dtsdy, dtsdz, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); } if (cost) { diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 45d78e44f..acd849f01 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -118,6 +118,163 @@ void warpx_push_ez_yee(int i, int j, int k, Array4 const& Ez, #endif } +static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, + Real &betaxy, Real &betaxz, Real &betayx, Real &betayz, Real &betazx, Real &betazy, + Real &gammax, Real &gammay, Real &gammaz, + Real &alphax, Real &alphay, Real &alphaz) +{ + // Cole-Karkkainen-Cowan push + // computes coefficients according to Cowan - PRST-AB 16, 041303 (2013) +#if (AMREX_SPACEDIM == 3) + Real delta = std::max( { dtsdx,dtsdy,dtsdz } ); + Real rx = (dtsdx/delta)*(dtsdx/delta); + Real ry = (dtsdy/delta)*(dtsdy/delta); + Real rz = (dtsdz/delta)*(dtsdz/delta); + Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry)); + betaxy = ry*beta; + betaxz = rz*beta; + betayx = rx*beta; + betayz = rz*beta; + betazx = rx*beta; + betazy = ry*beta; + gammax = ry*rz*(0.0625 - 0.125*ry*rz/(ry*rz + rz*rx + rx*ry)); + gammay = rx*rz*(0.0625 - 0.125*rx*rz/(ry*rz + rz*rx + rx*ry)); + gammaz = rx*ry*(0.0625 - 0.125*rx*ry/(ry*rz + rz*rx + rx*ry)); + alphax = 1. - 2.*betaxy - 2.*betaxz - 4.*gammax; + alphay = 1. - 2.*betayx - 2.*betayz - 4.*gammay; + alphaz = 1. - 2.*betazx - 2.*betazy - 4.*gammaz; + + betaxy = dtsdx*betaxy; + betaxz = dtsdx*betaxz; + betayx = dtsdy*betayx; + betayz = dtsdy*betayz; + betazx = dtsdz*betazx; + betazy = dtsdz*betazy; + alphax = dtsdx*alphax; + alphay = dtsdy*alphay; + alphaz = dtsdz*alphaz; + gammax = dtsdx*gammax; + gammay = dtsdy*gammay; + gammaz = dtsdz*gammaz; +#else + Real delta = std::max(dtsdx,dtsdz); + Real rx = (dtsdx/delta)*(dtsdx/delta); + Real rz = (dtsdz/delta)*(dtsdz/delta); + betaxz = 0.125*rz; + betazx = 0.125*rx; + alphax = 1. - 2.*betaxz; + alphaz = 1. - 2.*betazx; + + betaxz = dtsdx*betaxz; + betazx = dtsdz*betazx; + alphax = dtsdx*alphax; + alphaz = dtsdz*alphaz; +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bx_ckc(int j, int k, int l, Array4 const& Bx, + Array4 const& Ey, Array4 const& Ez, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + Bx(j,k,l) = Bx(j,k,l) - alphay * (Ez(j ,k+1,l ) - Ez(j, k ,l )) + - betayx * (Ez(j+1,k+1,l ) - Ez(j+1,k ,l ) + + Ez(j-1,k+1,l ) - Ez(j-1,k ,l )) + - betayz * (Ez(j ,k+1,l+1) - Ez(j ,k ,l+1) + + Ez(j ,k+1,l-1) - Ez(j ,k ,l-1)) + - gammay * (Ez(j+1,k+1,l+1) - Ez(j+1,k ,l+1) + + Ez(j-1,k+1,l+1) - Ez(j-1,k ,l+1) + + Ez(j+1,k+1,l-1) - Ez(j+1,k ,l-1) + + Ez(j-1,k+1,l-1) - Ez(j-1,k ,l-1)) + + alphaz * (Ey(j ,k ,l+1) - Ey(j, k, l )) + + betazx * (Ey(j+1,k ,l+1) - Ey(j+1,k ,l ) + + Ey(j-1,k ,l+1) - Ey(j-1,k ,l )) + + betazy * (Ey(j ,k+1,l+1) - Ey(j ,k+1,l ) + + Ey(j ,k-1,l+1) - Ey(j ,k-1,l )) + + gammaz * (Ey(j+1,k+1,l+1) - Ey(j+1,k+1,l ) + + Ey(j-1,k+1,l+1) - Ey(j-1,k+1,l ) + + Ey(j+1,k-1,l+1) - Ey(j+1,k-1,l ) + + Ey(j-1,k-1,l+1) - Ey(j-1,k-1,l )); +#else + Bx(j,k,0) = Bx(j,k,0) + alphaz * (Ey(j ,k+1,0) - Ey(j, k,0)) + + betazx * (Ey(j+1,k+1,0) - Ey(j+1,k,0) + + Ey(j-1,k+1,0) - Ey(j-1,k,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_by_ckc(int j, int k, int l, Array4 const& By, + Array4 const& Ex, Array4 const& Ez, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + By(j,k,l) = By(j,k,l) + alphax * (Ez(j+1,k ,l ) - Ez(j, k, l )) + + betaxy * (Ez(j+1,k+1,l ) - Ez(j ,k+1,l ) + + Ez(j+1,k-1,l ) - Ez(j ,k-1,l )) + + betaxz * (Ez(j+1,k ,l+1) - Ez(j ,k ,l+1) + + Ez(j+1,k ,l-1) - Ez(j ,k ,l-1)) + + gammax * (Ez(j+1,k+1,l+1) - Ez(j ,k+1,l+1) + + Ez(j+1,k-1,l+1) - Ez(j ,k-1,l+1) + + Ez(j+1,k+1,l-1) - Ez(j ,k+1,l-1) + + Ez(j+1,k-1,l-1) - Ez(j ,k-1,l-1)) + - alphaz * (Ex(j ,k ,l+1) - Ex(j ,k ,l )) + - betazx * (Ex(j+1,k ,l+1) - Ex(j+1,k ,l ) + + Ex(j-1,k ,l+1) - Ex(j-1,k ,l )) + - betazy * (Ex(j ,k+1,l+1) - Ex(j ,k+1,l ) + + Ex(j ,k-1,l+1) - Ex(j ,k-1,l )) + - gammaz * (Ex(j+1,k+1,l+1) - Ex(j+1,k+1,l ) + + Ex(j-1,k+1,l+1) - Ex(j-1,k+1,l ) + + Ex(j+1,k-1,l+1) - Ex(j+1,k-1,l ) + + Ex(j-1,k-1,l+1) - Ex(j-1,k-1,l )); +#else + By(j,k,0) = By(j,k,0) + alphax * (Ez(j+1,k ,0) - Ez(j,k ,0)) + + betaxz * (Ez(j+1,k+1,0) - Ez(j,k+1,0) + + Ez(j+1,k-1,0) - Ez(j,k-1,0)) + - alphaz * (Ex(j ,k+1,0) - Ex(j,k ,0)) + - betazx * (Ex(j+1,k+1,0) - Ex(j+1,k,0) + + Ex(j-1,k+1,0) - Ex(j-1,k,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_bz_ckc(int j, int k, int l, Array4 const& Bz, + Array4 const& Ex, Array4 const& Ey, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + Bz(j,k,l) = Bz(j,k,l) - alphax * (Ey(j+1,k ,l ) - Ey(j ,k ,l )) + - betaxy * (Ey(j+1,k+1,l ) - Ey(j ,k+1,l ) + + Ey(j+1,k-1,l ) - Ey(j ,k-1,l )) + - betaxz * (Ey(j+1,k ,l+1) - Ey(j ,k ,l+1) + + Ey(j+1,k ,l-1) - Ey(j ,k ,l-1)) + - gammax * (Ey(j+1,k+1,l+1) - Ey(j ,k+1,l+1) + + Ey(j+1,k-1,l+1) - Ey(j ,k-1,l+1) + + Ey(j+1,k+1,l-1) - Ey(j ,k+1,l-1) + + Ey(j+1,k-1,l-1) - Ey(j ,k-1,l-1)) + + alphay * (Ex(j ,k+1,l ) - Ex(j ,k ,l )) + + betayx * (Ex(j+1,k+1,l ) - Ex(j+1,k ,l ) + + Ex(j-1,k+1,l ) - Ex(j-1,k ,l )) + + betayz * (Ex(j ,k+1,l+1) - Ex(j ,k ,l+1) + + Ex(j ,k+1,l-1) - Ex(j ,k ,l-1)) + + gammay * (Ex(j+1,k+1,l+1) - Ex(j+1,k ,l+1) + + Ex(j-1,k+1,l+1) - Ex(j-1,k ,l+1) + + Ex(j+1,k+1,l-1) - Ex(j+1,k ,l-1) + + Ex(j-1,k+1,l-1) - Ex(j-1,k ,l-1)); + +#else + Bz(j,k,0) = Bz(j,k,0) - alphax * (Ey(j+1,k ,0) - Ey(j,k ,0)) + - betaxz * (Ey(j+1,k+1,0) - Ey(j,k+1,0) + + Ey(j+1,k-1,0) - Ey(j,k-1,0)); +#endif +} + AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_computedivb(int i, int j, int k, int dcomp, Array4 const& divB, Array4 const& Bx, Array4 const& By, Array4 const& Bz, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index fe2d51d30..f04b54acf 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -2,14 +2,11 @@ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb3d_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic -#define WRPX_PXR_PUSH_BVEC pxrpush_em3d_bvec -#define WRPX_PXR_PUSH_BVEC_CKC pxrpush_em3d_bvec_ckc #define WRPX_PXR_PUSH_EVEC_F pxrpush_em3d_evec_f #define WRPX_PXR_PUSH_EVEC_F_CKC pxrpush_em3d_evec_f_ckc #elif (AMREX_SPACEDIM == 2) -#define WRPX_PXR_PUSH_BVEC_CKC pxrpush_em2d_bvec_ckc #define WRPX_PXR_PUSH_EVEC_F pxrpush_em2d_evec_f #define WRPX_PXR_PUSH_EVEC_F_CKC pxrpush_em2d_evec_f_ckc @@ -19,13 +16,11 @@ #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz #define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho #define WRPX_PXR_RZ_VOLUME_SCALING_J apply_rz_volume_scaling_j -#define WRPX_PXR_PUSH_BVEC pxrpush_emrz_bvec #else #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2dxz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic_2d -#define WRPX_PXR_PUSH_BVEC pxrpush_em2d_bvec #endif @@ -526,75 +521,6 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n end subroutine warpx_particle_pusher_momenta - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the Maxwell solver (B field) - !> - !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds - !> where to apply the solver (PICSAR's valid cells) - !> @param[in] ex, ey, ez arrays for electric field - !> @param[inout] bx, by, bz arrays for magnetic field to update - !> @param[in] exlo, exhi, eylo, eyhi, ezlo, ezhi lower and higher bound for - !> the electric field arrays - !> @param[in] bxlo, bxhi, bylo, byhi, bzlo, bzhi lower and higher bound for - !> the magnetic field arrays - !> @param[in] dtsdx, dtsdy, dtsdz factors 0.5 * dt/(dx, dy, dz) - subroutine warpx_push_bvec( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - bx, bxlo, bxhi, & - by, bylo, byhi, & - bz, bzlo, bzhi, & - dtsdx, dtsdy, dtsdz, xmin, dx, & - maxwell_fdtd_solver_id) bind(C, name="warpx_push_bvec") - - integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), & - ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), & - exlo(BL_SPACEDIM), exhi(BL_SPACEDIM), eylo(BL_SPACEDIM), eyhi(BL_SPACEDIM), & - ezlo(BL_SPACEDIM), ezhi(BL_SPACEDIM), bxlo(BL_SPACEDIM), bxhi(BL_SPACEDIM), & - bylo(BL_SPACEDIM), byhi(BL_SPACEDIM), bzlo(BL_SPACEDIM), bzhi(BL_SPACEDIM), & - maxwell_fdtd_solver_id - - real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*) - - real(amrex_real), intent(IN):: bx(*), by(*), bz(*) - - real(amrex_real), intent(IN) :: dtsdx, dtsdy, dtsdz - - real(amrex_real), intent(IN) :: xmin, dx - - IF (maxwell_fdtd_solver_id .eq. 0) THEN - ! Yee FDTD solver - CALL WRPX_PXR_PUSH_BVEC( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - bx, bxlo, bxhi, & - by, bylo, byhi, & - bz, bzlo, bzhi, & - dtsdx,dtsdy,dtsdz & -#ifdef WARPX_RZ - ,xmin,dx & -#endif - ) - ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN - ! Cole-Karkkainen FDTD solver - CALL WRPX_PXR_PUSH_BVEC_CKC( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - bx, bxlo, bxhi, & - by, bylo, byhi, & - bz, bzlo, bzhi, & - dtsdx,dtsdy,dtsdz) - ENDIF - end subroutine warpx_push_bvec - ! _________________________________________________________________ !> !> @brief -- cgit v1.2.3 From 05bff4f6e3c1b2501ceae76bd15bfa2e48380146 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 5 Jul 2019 23:55:14 -0700 Subject: Converted evec_f routines to c++ --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 61 ++++++++++++++++---- Source/FieldSolver/WarpX_FDTD.H | 99 ++++++++++++++++++++++++++++++++ Source/FortranInterface/WarpX_picsar.F90 | 60 ------------------- 3 files changed, 150 insertions(+), 70 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 06b6b1b03..6bc674e06 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -322,16 +322,57 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) if (F) { - warpx_push_evec_f( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*F)[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, - &WarpX::maxwell_fdtd_solver_id); + auto const& Ffab = F->array(mfi); + if (WarpX::maxwell_fdtd_solver_id == 0) { + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_f_yee(j,k,l,Exfab,Ffab,dtsdx_c2); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_f_yee(j,k,l,Eyfab,Ffab,dtsdy_c2); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_f_yee(j,k,l,Ezfab,Ffab,dtsdz_c2); + }); + } + else if (WarpX::maxwell_fdtd_solver_id == 1) { + Real betaxy, betaxz, betayx, betayz, betazx, betazy; + Real gammax, gammay, gammaz; + Real alphax, alphay, alphaz; + warpx_calculate_ckc_coefficients(dtsdx_c2, dtsdy_c2, dtsdz_c2, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_f_ckc(j,k,l,Exfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_f_ckc(j,k,l,Eyfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_f_ckc(j,k,l,Ezfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + } } if (cost) { diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index acd849f01..dc1a9fae9 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -118,6 +118,37 @@ void warpx_push_ez_yee(int i, int j, int k, Array4 const& Ez, #endif } +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ex_f_yee(int j, int k, int l, Array4 const& Ex, + Array4 const& F, Real dtsdx_c2) +{ +#if (AMREX_SPACEDIM == 3) + Ex(j,k,l) = Ex(j,k,l) + dtsdx_c2 * (F(j+1,k,l) - F(j,k,l)); +#else + Ex(j,k,0) = Ex(j,k,0) + dtsdx_c2 * (F(j+1,k,0) - F(j ,k,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ey_f_yee(int j, int k, int l, Array4 const& Ey, + Array4 const& F, Real dtsdy_c2) +{ +#if (AMREX_SPACEDIM == 3) + Ey(j,k,l) = Ey(j,k,l) + dtsdy_c2 * (F(j,k+1,l) - F(j,k,l)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ez_f_yee(int j, int k, int l, Array4 const& Ez, + Array4 const& F, Real dtsdz_c2) +{ +#if (AMREX_SPACEDIM == 3) + Ez(j,k,l) = Ez(j,k,l) + dtsdz_c2 * (F(j,k,l+1) - F(j,k,l)); +#else + Ez(j,k,0) = Ez(j,k,0) + dtsdz_c2 * (F(j,k+1,0) - F(j,k,0)); +#endif +} + static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, Real &betaxy, Real &betaxz, Real &betayx, Real &betayz, Real &betazx, Real &betazy, Real &gammax, Real &gammay, Real &gammaz, @@ -275,6 +306,74 @@ void warpx_push_bz_ckc(int j, int k, int l, Array4 const& Bz, #endif } +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ex_f_ckc(int j, int k, int l, Array4 const& Ex, + Array4 const& F, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + Ex(j,k,l) = Ex(j,k,l) + alphax * (F(j+1,k ,l ) - F(j, k, l )) + + betaxy * (F(j+1,k+1,l ) - F(j,k+1,l ) + + F(j+1,k-1,l ) - F(j,k-1,l )) + + betaxz * (F(j+1,k ,l+1) - F(j,k ,l+1) + + F(j+1,k ,l-1) - F(j,k ,l-1)) + + gammax * (F(j+1,k+1,l+1) - F(j,k+1,l+1) + + F(j+1,k-1,l+1) - F(j,k-1,l+1) + + F(j+1,k+1,l-1) - F(j,k+1,l-1) + + F(j+1,k-1,l-1) - F(j,k-1,l-1)); +#else + Ex(j,k,0) = Ex(j,k,0) + alphax * (F(j+1,k ,0) - F(j,k ,0)) + + betaxz * (F(j+1,k+1,0) - F(j,k+1,0) + + F(j+1,k-1,0) - F(j,k-1,0)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ey_f_ckc(int j, int k, int l, Array4 const& Ey, + Array4 const& F, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + Ey(j,k,l) = Ey(j,k,l) + alphay * (F(j ,k+1,l ) - F(j ,k,l )) + + betayx * (F(j+1,k+1,l ) - F(j+1,k,l ) + + F(j-1,k+1,l ) - F(j-1,k,l )) + + betayz * (F(j ,k+1,l+1) - F(j ,k,l+1) + + F(j ,k+1,l-1) - F(j ,k,l-1)) + + gammay * (F(j+1,k+1,l+1) - F(j+1,k,l+1) + + F(j-1,k+1,l+1) - F(j-1,k,l+1) + + F(j+1,k+1,l-1) - F(j+1,k,l-1) + + F(j-1,k+1,l-1) - F(j-1,k,l-1)); +#endif +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void warpx_push_ez_f_ckc(int j, int k, int l, Array4 const& Ez, + Array4 const& F, + Real betaxy, Real betaxz, Real betayx, Real betayz, Real betazx, Real betazy, + Real gammax, Real gammay, Real gammaz, + Real alphax, Real alphay, Real alphaz) +{ +#if (AMREX_SPACEDIM == 3) + Ez(j,k,l) = Ez(j,k,l) + alphaz * (F(j ,k ,l+1) - F(j, k, l)) + + betazx * (F(j+1,k ,l+1) - F(j+1,k ,l) + + F(j-1,k ,l+1) - F(j-1,k ,l)) + + betazy * (F(j ,k+1,l+1) - F(j ,k+1,l) + + F(j ,k-1,l+1) - F(j ,k-1,l)) + + gammaz * (F(j+1,k+1,l+1) - F(j+1,k+1,l) + + F(j-1,k+1,l+1) - F(j-1,k+1,l) + + F(j+1,k-1,l+1) - F(j+1,k-1,l) + + F(j-1,k-1,l+1) - F(j-1,k-1,l)); +#else + Ez(j,k,0) = Ez(j,k,0) + alphaz * (F(j ,k+1,0) - F(j ,k,0)) + + betazx * (F(j+1,k+1,0) - F(j+1,k,0) + + F(j-1,k+1,0) - F(j-1,k,0)); +#endif +} + AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_computedivb(int i, int j, int k, int dcomp, Array4 const& divB, Array4 const& Bx, Array4 const& By, Array4 const& Bz, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index f04b54acf..dc47245dd 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -2,14 +2,9 @@ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb3d_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jxjyjz_generic -#define WRPX_PXR_PUSH_EVEC_F pxrpush_em3d_evec_f -#define WRPX_PXR_PUSH_EVEC_F_CKC pxrpush_em3d_evec_f_ckc #elif (AMREX_SPACEDIM == 2) -#define WRPX_PXR_PUSH_EVEC_F pxrpush_em2d_evec_f -#define WRPX_PXR_PUSH_EVEC_F_CKC pxrpush_em2d_evec_f_ckc - #ifdef WARPX_RZ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic @@ -521,59 +516,4 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n end subroutine warpx_particle_pusher_momenta - ! _________________________________________________________________ - !> - !> @brief - !> Subroutine for pushing E with the grad F terms - !> - !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds - !> where to apply the solver (PICSAR's valid cells) - !> @param[inout] ex, ey, ez arrays for electric field to update - !> @param[in] f array for the charge-correction term - !> @param[in] exlo, exhi, eylo, eyhi, ezlo, ezhi lower and higher bound for - !> the electric field arrays - !> @param[in] flo, fhi, lower and higher bound for the charge-correction term - !> @param[in] dtsdx_c2, dtsdy_c2, dtsdz_c2 factors dt/(dx, dy, dz)*c2 - subroutine warpx_push_evec_f( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - f, flo, fhi, & - dtsdx_c2, dtsdy_c2, dtsdz_c2, & - maxwell_fdtd_solver_id) bind(C, name="warpx_push_evec_f") - - integer(c_int), intent(in) :: xlo(BL_SPACEDIM), xhi(BL_SPACEDIM), & - ylo(BL_SPACEDIM), yhi(BL_SPACEDIM), zlo(BL_SPACEDIM), zhi(BL_SPACEDIM), & - exlo(BL_SPACEDIM), exhi(BL_SPACEDIM), eylo(BL_SPACEDIM), eyhi(BL_SPACEDIM), & - ezlo(BL_SPACEDIM), ezhi(BL_SPACEDIM), flo(BL_SPACEDIM), fhi(BL_SPACEDIM), & - maxwell_fdtd_solver_id - - real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*) - - real(amrex_real), intent(IN):: f - - real(amrex_real), intent(IN) :: dtsdx_c2, dtsdy_c2, dtsdz_c2 - - IF (maxwell_fdtd_solver_id .eq. 0) THEN - ! Yee FDTD solver - CALL WRPX_PXR_PUSH_EVEC_F( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - f, flo, fhi, & - dtsdx_c2, dtsdy_c2, dtsdz_c2) - ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN - ! Cole-Karkkainen FDTD solver - CALL WRPX_PXR_PUSH_EVEC_F_CKC( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - f, flo, fhi, & - dtsdx_c2, dtsdy_c2, dtsdz_c2) - ENDIF - end subroutine warpx_push_evec_f - end module warpx_to_pxr_module -- cgit v1.2.3 From 1434e6ef9cf3935ae31f94e09a71e6adf466d279 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Tue, 23 Jul 2019 09:07:24 -0700 Subject: For EB push conversion, added consts --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 12 ++++++------ Source/FieldSolver/WarpX_FDTD.H | 32 ++++++++++++++++---------------- 2 files changed, 22 insertions(+), 22 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6bc674e06..0887fa6a2 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -41,8 +41,8 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; - Real dxinv = 1./dx[0]; + const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; + const Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -224,10 +224,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; - Real dxinv = 1./dx[0]; + const Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; + const Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; if (patch_type == PatchType::fine) @@ -462,7 +462,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); const std::array dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index ee13875f6..e2d40f4c2 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -45,8 +45,8 @@ void warpx_push_bz_yee(int i, int j, int k, Array4 const& Bz, #elif defined WARPX_DIM_2D Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); #elif defined WARPX_DIM_RZ - Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); - Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); + const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); + const Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); Bz(i,j,0) += - dtsdx * (ru*Ey(i+1,j,0) - rd*Ey(i,j,0)); #endif } @@ -103,8 +103,8 @@ void warpx_push_ez_yee(int i, int j, int k, Array4 const& Ez, - mu_c2_dt * Jz(i,j,0); #elif defined WARPX_DIM_RZ if (i != 0 || rmin != 0.) { - Real ru = 1. + 0.5/(rmin*dxinv + i); - Real rd = 1. - 0.5/(rmin*dxinv + i); + const Real ru = 1. + 0.5/(rmin*dxinv + i); + const Real rd = 1. - 0.5/(rmin*dxinv + i); Ez(i,j,0) += + dtsdx_c2 * (ru*By(i,j,0) - rd*By(i-1,j,0)) - mu_c2_dt * Jz(i,j,0); } else { @@ -153,11 +153,11 @@ static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, // Cole-Karkkainen-Cowan push // computes coefficients according to Cowan - PRST-AB 16, 041303 (2013) #if defined WARPX_DIM_3D - Real delta = std::max( { dtsdx,dtsdy,dtsdz } ); - Real rx = (dtsdx/delta)*(dtsdx/delta); - Real ry = (dtsdy/delta)*(dtsdy/delta); - Real rz = (dtsdz/delta)*(dtsdz/delta); - Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry)); + const Real delta = std::max( { dtsdx,dtsdy,dtsdz } ); + const Real rx = (dtsdx/delta)*(dtsdx/delta); + const Real ry = (dtsdy/delta)*(dtsdy/delta); + const Real rz = (dtsdz/delta)*(dtsdz/delta); + const Real beta = 0.125*(1. - rx*ry*rz/(ry*rz + rz*rx + rx*ry)); betaxy = ry*beta; betaxz = rz*beta; betayx = rx*beta; @@ -184,9 +184,9 @@ static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, gammay *= dtsdy; gammaz *= dtsdz; #elif defined WARPX_DIM_2D - Real delta = std::max(dtsdx,dtsdz); - Real rx = (dtsdx/delta)*(dtsdx/delta); - Real rz = (dtsdz/delta)*(dtsdz/delta); + const Real delta = std::max(dtsdx,dtsdz); + const Real rx = (dtsdx/delta)*(dtsdx/delta); + const Real rz = (dtsdz/delta)*(dtsdz/delta); betaxz = 0.125*rz; betazx = 0.125*rx; alphax = 1. - 2.*betaxz; @@ -386,8 +386,8 @@ void warpx_computedivb(int i, int j, int k, int dcomp, Array4 const& divB, divB(i,j,0,dcomp) = (Bx(i+1,j ,0) - Bx(i,j,0))*dxinv + (Bz(i ,j+1,0) - Bz(i,j,0))*dzinv; #elif defined WARPX_DIM_RZ - Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); - Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); + const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); + const Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); divB(i,j,0,dcomp) = (ru*Bx(i+1,j,0) - rd*Bx(i,j,0))*dxinv + (Bz(i,j+1,0) - Bz(i,j,0))*dzinv; #endif @@ -417,8 +417,8 @@ void warpx_computedive(int i, int j, int k, int dcomp, Array4 const& divE, divE(i,j,0,dcomp) = 4.*Ex(i,j,0)*dxinv + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv; } else { - Real ru = 1. + 0.5/(rmin*dxinv + i); - Real rd = 1. - 0.5/(rmin*dxinv + i); + const Real ru = 1. + 0.5/(rmin*dxinv + i); + const Real rd = 1. - 0.5/(rmin*dxinv + i); divE(i,j,0,dcomp) = (ru*Ex(i,j,0) - rd*Ex(i-1,j,0))*dxinv + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv; } -- cgit v1.2.3