From 9e90cf995a1b68d77f94cf00983da6d8b1fe6710 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel <47607756+lgiacome@users.noreply.github.com> Date: Mon, 30 Aug 2021 21:38:04 +0200 Subject: ECT conformal solver (#1923) * adding base implementation of the conformal solver * adding some preprocessor directives * qualifying the isnan's correctly * getting rid of some unused fields * computing S_stab on the fly * using empty in length check * Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp * replaced a looponcpu with a parallelfor * trying to not pass lending/borrowing info by reference in evolveB * passing data using dataPtr * converting inds into a device vector * simplifying some += * extracting the inds BaseFabs as DataPtr * Revert "extracting the inds BaseFabs as DataPtr" This reverts commit bc709d2fcaa347c119514de79a3f57169a05af65. * implementing new way ogf handling cell extensions (gpu compatible) * fixing some white spaces * removed a typo * made a function static * tidying up * tidying up * getting rid of the rho vector * moving the rho update to evolveE * refectoring the cell extension * small chenge in input parameters * updating WarpX.H * bug fix * fixing another bug * rotating the cube * updated gitignore * using the conformal soler in the cube test * fixing the signature of a function * trying to fix the setVal problem * Revert "trying to fix the setVal problem" This reverts commit c7d0e5e3709b730ff570183b2a6df5f587ca4640. * trying to fix the setVal problem * cleaning some comments * removing an abort in device code * Including geometric information in the external field initializer * cleaned up the test * moving the geometric data initialization before the external fields initialization * changing way of saving info about intruded cells * fixed some white spaces * Using S_mod instead of S_red and S_enl * substituting a std::vector with amrex::Array1D * bug fix in the uint8_t -> coords maps * Condensed ect cell info into one single MultiFab * fixing some includes * fixing some more includes * fixed a typo * making some functions available in gpu code * using tilebox instead of validbox for the geometry initialization * communicating the geometric info in the guard cells * fixing the guard cells initialization for ect * fixing an unused parameter * fixing an unused parameter * ignored some unused variables when EB is not supported * cleaning up * cleaning up * ignored some unused variables when EB is not supported * evolving rho just for ect * handling some more unused variables * removing some white spaces * adding the rotated cube test * removed some white spaces * change cells into faces * small fix in unused parameters * fixed a comment * for safety for now just initialize the geometric data when lev==maxLevel * adding some preprocessor directives to isolate EB code * updating Makefile stuff * Improving the rotated cube test * bug fix in mesh refinement * fixing boundary conditions in rotated cube test * deactivating MPI in rotated cube test * activating mpi in test * improving outputs of cells extension * improved the documentation * Update Source/EmbeddedBoundary/WarpXInitEB.cpp Co-authored-by: Weiqun Zhang * Update Source/EmbeddedBoundary/WarpXInitEB.cpp Co-authored-by: Weiqun Zhang * Removed some obsolete isnan's * undone change to gitignore * adding some missing AMREX_GPU_DEVICEs * Changing some amrex::Real(0) into 0. Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Making some variables const Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Adding amrex:: to some ignore_unused Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Enhancing readability of some parallelfor's Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Removed some unecessary includes * Fixing some tags Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Fixing the function CountExtFaces * added a comment for Rhofield * Fixed a typo in CountExtCells * trying to remove accesses to private members in ComputeEdgeLengths * making some functions public * undoing some useless changes * undoing some useless changes * adding some AMREX_GPU_DEVICEs * adding some unused variables * adding some AMREX_GPU_DEVICEs to fix warnings * fixing relative error * Making several variables const Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Fixed a comment Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * changing some double to amrex::Real * removing commented lines * changing some double to amrex::Real * removing commented lines * inizialing nelems_x,y,z to avoid a warning * Fixing number of cells * Adding missing analysis routine in rotated cube test Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> * Making some variables constexpr Co-authored-by: Weiqun Zhang * Improving some reduce operations Co-authored-by: Weiqun Zhang * reading time directly from the output * fixed a few data types * fixing another int * Replacing or->||, and->&&, not->! * Adding license info * Adding a white space Co-authored-by: Axel Huebl * removed some unused imports * Moving the doxygen documentation * Including AMReX_LayoutData.H * removing some useless parameters * adding ect solver to the doc * Removing some useless reductions * Small change for consistency * Handling the resizing of arrays * Handling correctly the resing of arrays * Fixing some whitespaces * Fixing an indentation * Improving a comment * Removed a useless initialization * Renamed Rhofield to ECTRhofield * Renamed Rhofield to ECTRhofield * Added some if conditions to isolate some EB related code * Improved some comments Co-authored-by: Remi Lehe * Made some functions not member of WarpX anymore * Isolated some EB-only code * Isolated some EB-only code * Using the square brackets to access vectors * Ignoring some unused variables * Bug fix in initialization * Removed a redundant if * Using enumeration for connectivity and improved names * Added a comment * added a few comments * Removed some useless template parameters * Printing some info only in verbose mode * Revert "Removed some useless template parameters" This reverts commit 2c527980e6872c0212767c27f00e2b53ecbcfd0a. * Fixed a bug with the neighbours enumeration * Initializing geom data alsoo in the ghost cells Co-authored-by: Weiqun Zhang Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Co-authored-by: Axel Huebl Co-authored-by: Remi Lehe --- .../FieldSolver/FiniteDifferenceSolver/EvolveB.cpp | 233 +++++++++++++++++++++ 1 file changed, 233 insertions(+) (limited to 'Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp') 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 #include #include +#include #include #include @@ -49,8 +50,17 @@ void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr, 3 > const& Efield, std::unique_ptr const& Gfield, std::array< std::unique_ptr, 3 > const& face_areas, + std::array< std::unique_ptr, 3 > const& area_mod, + std::array< std::unique_ptr, 3 >& ECTRhofield, + std::array< std::unique_ptr, 3 >& Venl, + std::array< std::unique_ptr, 3 >& flag_info_cell, + std::array< std::unique_ptr >, 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 ( 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, 3 > const& Efield, + std::array< std::unique_ptr, 3 > const& edge_lengths, + std::array< std::unique_ptr, 3 > const& face_areas, + std::array< std::unique_ptr, 3 >& ECTRhofield, const int lev ) { +#ifdef AMREX_USE_EB + amrex::LayoutData* 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 const &Ex = Efield[0]->array(mfi); + Array4 const &Ey = Efield[1]->array(mfi); + Array4 const &Ez = Efield[2]->array(mfi); + Array4 const &Rhox = ECTRhofield[0]->array(mfi); + Array4 const &Rhoy = ECTRhofield[1]->array(mfi); + Array4 const &Rhoz = ECTRhofield[2]->array(mfi); + amrex::Array4 const &lx = edge_lengths[0]->array(mfi); + amrex::Array4 const &ly = edge_lengths[1]->array(mfi); + amrex::Array4 const &lz = edge_lengths[2]->array(mfi); + amrex::Array4 const &Sx = face_areas[0]->array(mfi); + amrex::Array4 const &Sy = face_areas[1]->array(mfi); + amrex::Array4 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, 3 >& Bfield, + std::array< std::unique_ptr, 3 > const& face_areas, + std::array< std::unique_ptr, 3 > const& area_mod, + std::array< std::unique_ptr, 3 >& ECTRhofield, + std::array< std::unique_ptr, 3 >& Venl, + std::array< std::unique_ptr, 3 >& flag_info_cell, + std::array< std::unique_ptr >, 3 >& borrowing, + const int lev, amrex::Real const dt ) { +#ifdef AMREX_USE_EB + amrex::LayoutData *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 const &B = Bfield[idim]->array(mfi); + Array4 const &Rho = ECTRhofield[idim]->array(mfi); + Array4 const &Venl_dim = Venl[idim]->array(mfi); + + amrex::Array4 const &flag_info_cell_dim = flag_info_cell[idim]->array(mfi); + amrex::Array4 const &S = face_areas[idim]->array(mfi); + amrex::Array4 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 -- cgit v1.2.3