diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
6 files changed, 282 insertions, 7 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index f5bba9dea..493e91dc9 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -31,6 +31,7 @@ #include <AMReX_LayoutData.H> #include <AMReX_MFIter.H> #include <AMReX_MultiFab.H> +#include <AMReX_iMultiFab.H> #include <AMReX_REAL.H> #include <AMReX_Utility.H> @@ -49,8 +50,17 @@ void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl, + std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell, + std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing, int lev, amrex::Real const dt ) { +#ifndef AMREX_USE_EB + amrex::ignore_unused(area_mod, ECTRhofield, Venl, flag_info_cell, borrowing); +#endif + // Select algorithm (The choice of algorithm is a runtime option, // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ @@ -69,7 +79,12 @@ void FiniteDifferenceSolver::EvolveB ( } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, Gfield, face_areas, lev, dt ); +#ifdef AMREX_USE_EB + } else if (m_fdtd_algo == MaxwellSolverAlgo::ECT) { + EvolveBCartesianECT(Bfield, face_areas, area_mod, ECTRhofield, Venl, flag_info_cell, + borrowing, lev, dt); +#endif #endif } else { amrex::Abort("EvolveB: Unknown algorithm"); @@ -197,6 +212,224 @@ void FiniteDifferenceSolver::EvolveBCartesian ( } } +void FiniteDifferenceSolver::EvolveRhoCartesianECT ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, const int lev ) { +#ifdef AMREX_USE_EB + amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev); + + // Loop through the grids, and over the tiles within each grid +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*ECTRhofield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) { + amrex::Gpu::synchronize(); + } + Real wt = amrex::second(); + + // Extract field data for this grid/tile + Array4<Real> const &Ex = Efield[0]->array(mfi); + Array4<Real> const &Ey = Efield[1]->array(mfi); + Array4<Real> const &Ez = Efield[2]->array(mfi); + Array4<Real> const &Rhox = ECTRhofield[0]->array(mfi); + Array4<Real> const &Rhoy = ECTRhofield[1]->array(mfi); + Array4<Real> const &Rhoz = ECTRhofield[2]->array(mfi); + amrex::Array4<amrex::Real> const &lx = edge_lengths[0]->array(mfi); + amrex::Array4<amrex::Real> const &ly = edge_lengths[1]->array(mfi); + amrex::Array4<amrex::Real> const &lz = edge_lengths[2]->array(mfi); + amrex::Array4<amrex::Real> const &Sx = face_areas[0]->array(mfi); + amrex::Array4<amrex::Real> const &Sy = face_areas[1]->array(mfi); + amrex::Array4<amrex::Real> const &Sz = face_areas[2]->array(mfi); + + // Extract tileboxes for which to loop + Box const &trhox = mfi.tilebox(ECTRhofield[0]->ixType().toIntVect()); + Box const &trhoy = mfi.tilebox(ECTRhofield[1]->ixType().toIntVect()); + Box const &trhoz = mfi.tilebox(ECTRhofield[2]->ixType().toIntVect()); + + amrex::ParallelFor(trhox, trhoy, trhoz, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + if (Sx(i, j, k) <= 0) return; + + Rhox(i, j, k) = (Ey(i, j, k) * ly(i, j, k) - Ey(i, j, k + 1) * ly(i, j, k + 1) + + Ez(i, j + 1, k) * lz(i, j + 1, k) - Ez(i, j, k) * lz(i, j, k)) / Sx(i, j, k); + + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + if (Sy(i, j, k) <= 0) return; + + Rhoy(i, j, k) = (Ez(i, j, k) * lz(i, j, k) - Ez(i + 1, j, k) * lz(i + 1, j, k) + + Ex(i, j, k + 1) * lx(i, j, k + 1) - Ex(i, j, k) * lx(i, j, k)) / Sy(i, j, k); + + }, + + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + if (Sz(i, j, k) <= 0) return; + + Rhoz(i, j, k) = (Ex(i, j, k) * lx(i, j, k) - Ex(i, j + 1, k) * lx(i, j + 1, k) + + Ey(i + 1, j, k) * ly(i + 1, j, k) - Ey(i, j, k) * ly(i, j, k)) / Sz(i, j, k); + + } + ); + + if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) + { + amrex::Gpu::synchronize(); + wt = amrex::second() - wt; + amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt); + } + } +#else + amrex::ignore_unused(Efield, edge_lengths, face_areas, ECTRhofield, lev); +#endif +} + +void FiniteDifferenceSolver::EvolveBCartesianECT ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl, + std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell, + std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing, + const int lev, amrex::Real const dt ) { +#ifdef AMREX_USE_EB + amrex::LayoutData<amrex::Real> *cost = WarpX::getCosts(lev); + + Venl[0]->setVal(0.); + Venl[1]->setVal(0.); + Venl[2]->setVal(0.); + + // Loop through the grids, and over the tiles within each grid +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield[0]); mfi.isValid(); ++mfi) { + + if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) { + amrex::Gpu::synchronize(); + } + Real wt = amrex::second(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Extract field data for this grid/tile + Array4<Real> const &B = Bfield[idim]->array(mfi); + Array4<Real> const &Rho = ECTRhofield[idim]->array(mfi); + Array4<Real> const &Venl_dim = Venl[idim]->array(mfi); + + amrex::Array4<int> const &flag_info_cell_dim = flag_info_cell[idim]->array(mfi); + amrex::Array4<Real> const &S = face_areas[idim]->array(mfi); + amrex::Array4<Real> const &S_mod = area_mod[idim]->array(mfi); + + auto &borrowing_dim = (*borrowing[idim])[mfi]; + auto borrowing_dim_neigh_faces = borrowing_dim.neigh_faces.data(); + auto borrowing_dim_area = borrowing_dim.area.data(); + + auto const &borrowing_inds = (*borrowing[idim])[mfi].inds.data(); + auto const &borrowing_size = (*borrowing[idim])[mfi].size.array(); + auto const &borrowing_inds_pointer = (*borrowing[idim])[mfi].inds_pointer.array(); + + // Extract tileboxes for which to loop + Box const &tb = mfi.tilebox(Bfield[idim]->ixType().toIntVect()); + + //Take care of the unstable cells + amrex::ParallelFor(tb, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + + if (S(i, j, k) <= 0) return; + + if (!(flag_info_cell_dim(i, j, k) == 0)) + return; + + Venl_dim(i, j, k) = Rho(i, j, k) * S(i, j, k); + amrex::Real rho_enl; + + // First we compute the rho of the enlarged face + for (int offset = 0; offset<borrowing_size(i, j, k); offset++) { + int ind = borrowing_inds[*borrowing_inds_pointer(i, j, k) + offset]; + auto vec = FaceInfoBox::uint8_to_inds(borrowing_dim_neigh_faces[ind]); + int ip, jp, kp; + if(idim == 0){ + ip = i; + jp = j + vec(0); + kp = k + vec(1); + }else if(idim == 1){ + ip = i + vec(0); + jp = j; + kp = k + vec(1); + }else{ + ip = i + vec(0); + jp = j + vec(1); + kp = k; + } + + Venl_dim(i, j, k) += Rho(ip, jp, kp) * borrowing_dim_area[ind]; + + } + + rho_enl = Venl_dim(i, j, k) / S_mod(i, j, k); + + for (int offset = 0; offset < borrowing_size(i, j, k); offset++) { + int ind = borrowing_inds[*borrowing_inds_pointer(i, j, k) + offset]; + auto vec = FaceInfoBox::uint8_to_inds(borrowing_dim_neigh_faces[ind]); + int ip, jp, kp; + if(idim == 0){ + ip = i; + jp = j + vec(0); + kp = k + vec(1); + }else if(idim == 1){ + ip = i + vec(0); + jp = j; + kp = k + vec(1); + }else{ + ip = i + vec(0); + jp = j + vec(1); + kp = k; + } + + Venl_dim(ip, jp, kp) += rho_enl * borrowing_dim_area[ind]; + + } + + B(i, j, k) = B(i, j, k) - dt * rho_enl; + + }); + + //Take care of the stable cells + amrex::ParallelFor(tb, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + if (S(i, j, k) <= 0) return; + + if (flag_info_cell_dim(i, j, k) == 0) { + return; + } + else if (flag_info_cell_dim(i, j, k) == 1) { + //Stable cell which hasn't been intruded + B(i, j, k) = B(i, j, k) - dt * Rho(i, j, k); + } else if (flag_info_cell_dim(i, j, k) == 2) { + //Stable cell which has been intruded + Venl_dim(i, j, k) += Rho(i, j, k) * S_mod(i, j, k); + B(i, j, k) = B(i, j, k) - dt * Venl_dim(i, j, k) / S(i, j, k); + } + + }); + + } + if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) + { + amrex::Gpu::synchronize(); + wt = amrex::second() - wt; + amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt); + } + } +#else + amrex::ignore_unused(Bfield, face_areas, area_mod, ECTRhofield, Venl, flag_info_cell, borrowing, + lev, dt); +#endif +} + #else // corresponds to ifndef WARPX_DIM_RZ template<typename T_Algo> diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp index 299f92f2f..5b182812a 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp @@ -55,7 +55,7 @@ void FiniteDifferenceSolver::EvolveBPML ( EvolveBPMLCartesian <CartesianNodalAlgorithm> (Bfield, Efield, dt, dive_cleaning); - } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) { EvolveBPMLCartesian <CartesianYeeAlgorithm> (Bfield, Efield, dt, dive_cleaning); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 241992539..ac19984cd 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -49,11 +49,21 @@ void FiniteDifferenceSolver::EvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, std::unique_ptr<amrex::MultiFab> const& Ffield, int lev, amrex::Real const dt ) { - // Select algorithm (The choice of algorithm is a runtime option, - // but we compile code for each algorithm, using templates) +#ifdef AMREX_USE_EB + if (m_fdtd_algo != MaxwellSolverAlgo::ECT) { + amrex::ignore_unused(face_areas, ECTRhofield); + } +#else + amrex::ignore_unused(face_areas, ECTRhofield); +#endif + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) #ifdef WARPX_DIM_RZ if (m_fdtd_algo == MaxwellSolverAlgo::Yee){ ignore_unused(edge_lengths); @@ -63,10 +73,14 @@ void FiniteDifferenceSolver::EvolveE ( EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); - } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) { EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); - +#ifdef AMREX_USE_EB + if (m_fdtd_algo == MaxwellSolverAlgo::ECT) { + EvolveRhoCartesianECT(Efield, edge_lengths, face_areas, ECTRhofield, lev); + } +#endif } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); @@ -157,13 +171,16 @@ void FiniteDifferenceSolver::EvolveECartesian ( // Skip field push if this cell is fully covered by embedded boundaries if (ly(i,j,k) <= 0) return; #endif + Ey(i, j, k) += c2 * dt * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k) - PhysConst::mu0 * jy(i, j, k) ); + }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ + #ifdef AMREX_USE_EB // Skip field push if this cell is fully covered by embedded boundaries if (lz(i,j,k) <= 0) return; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp index ea5ecb381..b182a836b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp @@ -61,7 +61,7 @@ void FiniteDifferenceSolver::EvolveEPML ( EvolveEPMLCartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles ); - } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) { EvolveEPMLCartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles ); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 820f2b876..2a96ae44d 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -8,6 +8,7 @@ #ifndef WARPX_FINITE_DIFFERENCE_SOLVER_H_ #define WARPX_FINITE_DIFFERENCE_SOLVER_H_ +#include "EmbeddedBoundary/WarpXFaceInfoBox.H" #include "FiniteDifferenceSolver_fwd.H" #include "BoundaryConditions/PML_fwd.H" @@ -49,12 +50,19 @@ class FiniteDifferenceSolver std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, std::unique_ptr<amrex::MultiFab> const& Gfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl, + std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell, + std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing, int lev, amrex::Real const dt ); void EvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield, std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, std::unique_ptr<amrex::MultiFab> const& Ffield, int lev, amrex::Real const dt ); @@ -200,6 +208,23 @@ class FiniteDifferenceSolver std::array<std::unique_ptr<amrex::MultiFab>,3> const& Bfield, amrex::Real const dt); + void EvolveRhoCartesianECT ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, int lev); + + void EvolveBCartesianECT ( + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& face_areas, + std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& area_mod, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& ECTRhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Venl, + std::array< std::unique_ptr<amrex::iMultiFab>, 3 >& flag_info_cell, + std::array< std::unique_ptr<amrex::LayoutData<FaceInfoBox> >, 3 >& borrowing, + int lev, amrex::Real const dt + ); + template< typename T_Algo > void ComputeDivECartesian ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp index aab4a69c5..cf6a603c7 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -65,7 +65,7 @@ FiniteDifferenceSolver::FiniteDifferenceSolver ( CartesianNodalAlgorithm::InitializeStencilCoefficients( cell_size, m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z ); - } else if (fdtd_algo == MaxwellSolverAlgo::Yee) { + } else if (fdtd_algo == MaxwellSolverAlgo::Yee || fdtd_algo == MaxwellSolverAlgo::ECT) { CartesianYeeAlgorithm::InitializeStencilCoefficients( cell_size, m_h_stencil_coefs_x, m_h_stencil_coefs_y, m_h_stencil_coefs_z ); |