From 2c30b530e439787c733dd922cf95daf384dd94de Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 12:41:59 +0200 Subject: [PATCH 01/35] Added staircased embedded boundaris to the YEE solver --- Source/EmbeddedBoundary/CMakeLists.txt | 1 + Source/EmbeddedBoundary/WarpXInitEB.cpp | 166 ++++++++++++++++++ Source/Evolve/WarpXEvolve.cpp | 13 ++ .../FiniteDifferenceSolver/EvolveB.cpp | 23 ++- .../FiniteDifferenceSolver/EvolveE.cpp | 23 ++- .../FiniteDifferenceSolver.H | 4 + Source/FieldSolver/WarpXPushFieldsEM.cpp | 8 +- Source/Initialization/WarpXInitData.cpp | 7 + Source/WarpX.H | 17 ++ Source/WarpX.cpp | 19 +- 10 files changed, 270 insertions(+), 11 deletions(-) diff --git a/Source/EmbeddedBoundary/CMakeLists.txt b/Source/EmbeddedBoundary/CMakeLists.txt index 65a11610d18..7e11e2930f9 100644 --- a/Source/EmbeddedBoundary/CMakeLists.txt +++ b/Source/EmbeddedBoundary/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(WarpX PRIVATE WarpXInitEB.cpp + WarpXInitEMField.cpp ) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index ec17d9efccc..440f80bc794 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -18,5 +18,171 @@ WarpX::InitEB () pp_eb2.add("geom_type", "all_regular"); // use all_regular by default } amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel()); + #endif } + +void +WarpX::ComputeEdgeLengths() { + BL_PROFILE("ComputeEdgeLengths"); + + auto const eb_fact = fieldEBFactory(maxLevel()); + + auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &edge_centroid = eb_fact.getEdgeCent(); +// for ( amrex::MFIter mfi(*Bfield_fp[maxLevel()][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + amrex::Box const &box = mfi.validbox(); + amrex::FabType fab_type = flags[mfi].getType(box); + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + } else { + // mix of regular, cut and covered + // Note that edge_centroid[idim] is a MultiCutFab that does not have + // data if the fab type regular or covered. + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const &edge = edge_centroid[idim]->const_array(mfi); + auto const &edge_lengths_dim = edge_lengths[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), + [=](int i, int j, int k) { + if (edge(i, j, k) == amrex::Real(-1.0)) { + // This edge is all covered + edge_lengths_dim(i, j, k) = 0.; + } else if (edge(i, j, k) == amrex::Real(1.0)) { + // This edge is all open + edge_lengths_dim(i, j, k) = 1.; + } else { + // This edge is cut. + amrex::Real edge_cent = edge(i, j, k); // edge centroid: (-0.5,0.5) + if (edge_cent < amrex::Real(0.0)) { + // The right side is covered. + edge_lengths_dim(i, j, k) = + amrex::Real(2.0) * edge_cent + amrex::Real(0.5) + 0.5; // (0, 1) + } else { + // The left side is covered + // TODO: rewrite well + edge_lengths_dim(i, j, k) = 1 - + (amrex::Real(2.0) * edge_cent - amrex::Real(0.5) + 0.5); // (0, 1) + } + } + }); + } + } + + } + +} + +void +WarpX::ComputeFaceAreas() { + BL_PROFILE("ComputeFaceAreas"); + + auto const eb_fact = fieldEBFactory(maxLevel()); + + auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &area_frac = eb_fact.getAreaFrac(); + + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + amrex::Box const &box = mfi.validbox(); + amrex::FabType fab_type = flags[mfi].getType(box); + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + } else { + // mix of regular, cut and covered + // Note that edge_centroid[idim] is a MultiCutFab that does not have + // data if the fab type regular or covered. + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = face_areas[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = face(i, j, k); + }); + } + } + + } + +} + +void +WarpX::ScaleEdges() { + auto const& cell_size = CellSize(maxLevel()); + + auto const eb_fact = fieldEBFactory(maxLevel()); + + auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &edge_centroid = eb_fact.getEdgeCent(); + + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + amrex::Box const &box = mfi.validbox(); + const auto lo = amrex::lbound(box); + const auto hi = amrex::ubound(box); + amrex::FabType fab_type = flags[mfi].getType(box); + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + } else { + // mix of regular, cut and covered + // Note that edge_centroid[idim] is a MultiCutFab that does not have + // data if the fab type regular or covered. + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const &edge = edge_centroid[1]->const_array(mfi); + auto const &edge_lengths_dim = edge_lengths[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; + } + } + } + } + } +} + +void +WarpX::ScaleAreas() { + auto const& cell_size = CellSize(maxLevel()); + amrex::Real full_area; + + auto const eb_fact = fieldEBFactory(maxLevel()); + + auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &area_frac = eb_fact.getAreaFrac(); + + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + amrex::Box const &box = mfi.validbox(); + amrex::FabType fab_type = flags[mfi].getType(box); + const auto lo = amrex::lbound(box); + const auto hi = amrex::ubound(box); + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + } else { + // mix of regular, cut and covered + // Note that edge_centroid[idim] is a MultiCutFab that does not have + // data if the fab type regular or covered. + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (idim == 0) { + full_area = cell_size[1]*cell_size[2]; + } else if (idim == 1) { + full_area = cell_size[0]*cell_size[2]; + } else { + full_area = cell_size[0]*cell_size[1]; + } + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = face_areas[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + face_areas_dim(i,j,hi.z/2) *= full_area; + } + } + } + } + } +} diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 39b51b26140..a92bf569289 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -45,6 +45,19 @@ WarpX::Evolve (int numsteps) bool early_params_checked = false; // check typos in inputs after step 1 Real walltime, walltime_start = amrex::second(); + +#ifdef AMREX_USE_EB + //If necessary initialize the em fields + std::string geom_type; + amrex::ParmParse pp_eb2("eb2"); + bool init_em_field; + pp_eb2.get("geom_type", geom_type); + pp_eb2.get("init_em_field", init_em_field); + if(init_em_field and geom_type=="sphere") { + InitEMFieldSphere(); + } +#endif + for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) { Real walltime_beg_step = amrex::second(); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 2f66796f3e1..b68450a37ef 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -25,6 +25,7 @@ using namespace amrex; void FiniteDifferenceSolver::EvolveB ( std::array< std::unique_ptr, 3 >& Bfield, std::array< std::unique_ptr, 3 > const& Efield, + std::array< std::unique_ptr, 3 > const& face_areas, int lev, amrex::Real const dt ) { // Select algorithm (The choice of algorithm is a runtime option, @@ -37,15 +38,15 @@ void FiniteDifferenceSolver::EvolveB ( #else if (m_do_nodal) { - EvolveBCartesian ( Bfield, Efield, lev, dt ); + EvolveBCartesian ( Bfield, Efield, face_areas, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { - EvolveBCartesian ( Bfield, Efield, lev, dt ); + EvolveBCartesian ( Bfield, Efield, face_areas, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - EvolveBCartesian ( Bfield, Efield, lev, dt ); + EvolveBCartesian ( Bfield, Efield, face_areas, lev, dt ); #endif } else { @@ -61,6 +62,7 @@ template void FiniteDifferenceSolver::EvolveBCartesian ( std::array< std::unique_ptr, 3 >& Bfield, std::array< std::unique_ptr, 3 > const& Efield, + std::array< std::unique_ptr, 3 > const& face_areas, int lev, amrex::Real const dt ) { amrex::LayoutData* cost = WarpX::getCosts(lev); @@ -84,6 +86,12 @@ void FiniteDifferenceSolver::EvolveBCartesian ( Array4 const& Ey = Efield[1]->array(mfi); Array4 const& Ez = Efield[2]->array(mfi); +#ifdef AMREX_USE_EB + Array4 const& Sx = face_areas[0]->array(mfi); + Array4 const& Sy = face_areas[1]->array(mfi); + Array4 const& Sz = face_areas[2]->array(mfi); +#endif + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); int const n_coefs_x = m_stencil_coefs_x.size(); @@ -101,16 +109,25 @@ void FiniteDifferenceSolver::EvolveBCartesian ( amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + if (Sx(i, j, k) <= 0) return; +#endif Bx(i, j, k) += dt * T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k) - dt * T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + if (Sy(i, j, k) <= 0) return; +#endif By(i, j, k) += dt * T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k) - dt * T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + if (Sz(i, j, k) <= 0) return; +#endif Bz(i, j, k) += dt * T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k) - dt * T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index cd1e3bdf9d3..52ad487caf0 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -28,6 +28,7 @@ void FiniteDifferenceSolver::EvolveE ( std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ) { @@ -41,15 +42,15 @@ void FiniteDifferenceSolver::EvolveE ( #else if (m_do_nodal) { - EvolveECartesian ( Efield, Bfield, Jfield, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { - EvolveECartesian ( Efield, Bfield, Jfield, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) { - EvolveECartesian ( Efield, Bfield, Jfield, Ffield, lev, dt ); + EvolveECartesian ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt ); #endif } else { @@ -66,6 +67,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ) { @@ -94,6 +96,12 @@ void FiniteDifferenceSolver::EvolveECartesian ( Array4 const& jy = Jfield[1]->array(mfi); Array4 const& jz = Jfield[2]->array(mfi); +#ifdef AMREX_USE_EB + Array4 const& lx = edge_lengths[0]->array(mfi); + Array4 const& ly = edge_lengths[1]->array(mfi); + Array4 const& lz = edge_lengths[2]->array(mfi); +#endif + // Extract stencil coefficients Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr(); int const n_coefs_x = m_stencil_coefs_x.size(); @@ -111,6 +119,9 @@ void FiniteDifferenceSolver::EvolveECartesian ( amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + if (lx(i, j, k) <= 0) return; +#endif Ex(i, j, k) += c2 * dt * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k) + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k) @@ -118,6 +129,9 @@ void FiniteDifferenceSolver::EvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + 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) @@ -125,6 +139,9 @@ void FiniteDifferenceSolver::EvolveECartesian ( }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#ifdef AMREX_USE_EB + if (lz(i,j,k) <= 0) return; +#endif Ez(i, j, k) += c2 * dt * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k) + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index c5041a53f1d..f434ad3d04d 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -38,11 +38,13 @@ class FiniteDifferenceSolver void EvolveB ( std::array< std::unique_ptr, 3 >& Bfield, std::array< std::unique_ptr, 3 > const& Efield, + std::array< std::unique_ptr, 3 > const& face_areas, int lev, amrex::Real const dt ); void EvolveE ( std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ); @@ -150,6 +152,7 @@ class FiniteDifferenceSolver void EvolveBCartesian ( std::array< std::unique_ptr, 3 >& Bfield, std::array< std::unique_ptr, 3 > const& Efield, + std::array< std::unique_ptr, 3 > const& face_areas, int lev, amrex::Real const dt ); template< typename T_Algo > @@ -157,6 +160,7 @@ class FiniteDifferenceSolver std::array< std::unique_ptr, 3 >& Efield, std::array< std::unique_ptr, 3 > const& Bfield, std::array< std::unique_ptr, 3 > const& Jfield, + std::array< std::unique_ptr, 3 > const& edge_lengths, std::unique_ptr const& Ffield, int lev, amrex::Real const dt ); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 120dc8b5c29..1960a5471b2 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -192,9 +192,9 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) // Evolve B field in regular cells if (patch_type == PatchType::fine) { - m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], lev, a_dt ); + m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], face_areas[lev], lev, a_dt ); } else { - m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], lev, a_dt ); + m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], face_areas[lev], lev, a_dt ); } // Evolve B field in PML cells @@ -243,10 +243,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) // Evolve E field in regular cells if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveE( Efield_fp[lev], Bfield_fp[lev], - current_fp[lev], F_fp[lev], lev, a_dt ); + current_fp[lev], edge_lengths[lev], F_fp[lev], lev, a_dt ); } else { m_fdtd_solver_cp[lev]->EvolveE( Efield_cp[lev], Bfield_cp[lev], - current_cp[lev], F_cp[lev], lev, a_dt ); + current_cp[lev], edge_lengths[lev], F_cp[lev], lev, a_dt ); } // Evolve E field in PML cells diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 1852d871f9d..68e137480f5 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -86,6 +86,13 @@ WarpX::InitData () PostRestart(); } +#ifdef AMREX_USE_EB + ComputeEdgeLengths(); + ComputeFaceAreas(); + ScaleEdges(); + ScaleAreas(); +#endif + ComputePMLFactors(); if (WarpX::use_fdtd_nci_corr) { diff --git a/Source/WarpX.H b/Source/WarpX.H index 2c4b9c20e6f..9a0c28c6ac5 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -802,6 +802,13 @@ private: amrex::Vector, 3 > > Bfield_fp; amrex::Vector, 3 > > Efield_avg_fp; amrex::Vector, 3 > > Bfield_avg_fp; + +#ifdef AMREX_USE_EB + //EB grid info + amrex::Vector, 3 > > edge_lengths; + amrex::Vector, 3 > > face_areas; +#endif + // store fine patch amrex::Vector, 3 > > current_store; @@ -958,6 +965,16 @@ private: void InitEB (); +#ifdef AMREX_USE_EB + void ComputeEdgeLengths(); + void ComputeFaceAreas(); + void ScaleEdges(); + void ScaleAreas(); + + static amrex::RealArray AnalyticSolSphere(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere); + void InitEMFieldSphere(); +#endif + private: // void EvolvePSATD (int numsteps); void PushPSATD (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 92a2efb6113..022a4520667 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -173,7 +173,6 @@ WarpX::WarpX () InitEB(); // Geometry on all levels has been defined already. - // No valid BoxArray and DistributionMapping have been defined. // But the arrays for them have been resized. @@ -227,6 +226,11 @@ WarpX::WarpX () Efield_avg_fp.resize(nlevs_max); Bfield_avg_fp.resize(nlevs_max); +#ifdef AMREX_USE_EB + edge_lengths.resize(nlevs_max); + face_areas.resize(nlevs_max); +#endif + current_store.resize(nlevs_max); F_cp.resize(nlevs_max); @@ -1217,6 +1221,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm Efield_avg_fp[lev][1] = std::make_unique(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE,tag("Efield_avg_fp[y]")); Efield_avg_fp[lev][2] = std::make_unique(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE,tag("Efield_avg_fp[z]")); +#ifdef AMREX_USE_EB + // EB info are needed only at the coarsest level + if (lev == maxLevel()) + { + edge_lengths[lev][0] = std::make_unique(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[x]")); + edge_lengths[lev][1] = std::make_unique(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[y]")); + edge_lengths[lev][2] = std::make_unique(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[z]")); + face_areas[lev][0] = std::make_unique(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE,tag("face_areas[x]")); + face_areas[lev][1] = std::make_unique(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE,tag("face_areas[y]")); + face_areas[lev][2] = std::make_unique(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE,tag("face_areas[z]")); + } +#endif + bool deposit_charge = do_dive_cleaning || (plot_rho && do_back_transformed_diagnostics); if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { deposit_charge = do_dive_cleaning || (plot_rho && do_back_transformed_diagnostics) From 9d165116f4fdc81bcc066cbcaaeda4d4e7e79e01 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 12:47:38 +0200 Subject: [PATCH 02/35] adding spherical resonating cavity test --- .../embedded_boundary_sphere/inputs_3d | 23 ++++++++ .../plots/plot_fields.py | 58 +++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 Examples/Modules/embedded_boundary_sphere/inputs_3d create mode 100644 Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py diff --git a/Examples/Modules/embedded_boundary_sphere/inputs_3d b/Examples/Modules/embedded_boundary_sphere/inputs_3d new file mode 100644 index 00000000000..a9513b62ce6 --- /dev/null +++ b/Examples/Modules/embedded_boundary_sphere/inputs_3d @@ -0,0 +1,23 @@ +max_step = 1000 +amr.n_cell = 48 48 48 +amr.max_grid_size = 128 +amr.max_level = 0 + +geometry.coord_sys = 0 +geometry.is_periodic = 0 0 0 +geometry.prob_lo = -0.5 -0.5 -0.5 +geometry.prob_hi = 0.5 0.5 0.5 +warpx.do_pml = 0 +warpx.const_dt = 1e-6 +warpx.cfl = 1 + +eb2.geom_type = sphere +eb2.sphere_center = 0. 0. 0. +eb2.sphere_radius = 0.3 +eb2.sphere_has_fluid_inside = true +eb2.init_em_field = true + +diagnostics.diags_names = diag +diag.intervals = 1 +diag.diag_type = Full +diag.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho diff --git a/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py b/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py new file mode 100644 index 00000000000..7228456a055 --- /dev/null +++ b/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py @@ -0,0 +1,58 @@ +import yt +import os +from scipy.constants import mu_0 +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import AxesGrid + +fig = plt.figure(1, (1,15)) + +grid = AxesGrid(fig, (0.075,0.075,0.85,0.85), + nrows_ncols = (1, 3), + axes_pad = 1.0, + label_mode = "1", + share_all = True, + cbar_location="right", + cbar_mode="each", + cbar_size="3%", + cbar_pad="0%") + + +path, dirs, files = next(os.walk("../diags/")) +N_dirs = len(dirs) + +for i in range(0, N_dirs): + filename = '../diags/diag'+str(int(i)).zfill(5)+'/' + print(filename) + ds = yt.load(filename) + comp = 'Ez' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -5000, zmax = 5000) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[0].axes + plot.cax = grid.cbar_axes[0] + p._setup_plots() + + comp = 'Bx' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[1].axes + plot.cax = grid.cbar_axes[1] + p._setup_plots() + + comp = 'By' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[2].axes + plot.cax = grid.cbar_axes[2] + p._setup_plots() + + fig.set_size_inches(13,5) + p.save(str(i)) From 6917ee0e5330a45fd543405570020fcbfce16deb Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 13:26:41 +0200 Subject: [PATCH 03/35] adding functions for fields initialization --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 76 ++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 Source/EmbeddedBoundary/WarpXInitEMField.cpp diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp new file mode 100644 index 00000000000..febf41ac173 --- /dev/null +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -0,0 +1,76 @@ + +#include "WarpX.H" +#include + +amrex::RealArray +WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { + amrex::Real k=2.7437/r_sphere; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = PhysConst::mu0; + + amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); + + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + + return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; +} + +void +WarpX::InitEMFieldSphere (){ + amrex::Real r_sphere; + amrex::ParmParse pp_eb2("eb2"); + amrex::ParmParse pp_geometry("geometry"); + pp_eb2.get("sphere_radius", r_sphere); + amrex::RealArray prob_lo; + pp_geometry.get("prob_lo", prob_lo); + auto &Bfield = Bfield_fp[maxLevel()]; + auto &Efield = Efield_fp[maxLevel()]; + + auto &face_areas_0 = m_face_areas[maxLevel()]; + for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + amrex::Array4 const &Bx = Bfield[0]->array(mfi); + amrex::Array4 const &By = Bfield[1]->array(mfi); + amrex::Array4 const &Bz = Bfield[2]->array(mfi); + amrex::Array4 const &Sx = face_areas_0[0]->array(mfi); + amrex::Array4 const &Sy = face_areas_0[1]->array(mfi); + amrex::Array4 const &Sz = face_areas_0[2]->array(mfi); + amrex::Box const &box = mfi.validbox(); + const auto hi = amrex::ubound(box); + amrex::Real x, y, z; + for (int i = 0; i < hi.x; i++) { + for (int j = 0; j < hi.y; j++) { + for (int k = 0; k < hi.z; k++) { + if (Sx(i, j, k) > 0) { + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = k * CellSize(maxLevel())[2] + prob_lo[2]; + + Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; + } + + if (Sy(i, j, k) > 0) { + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = j * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + + By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; + } + + if (Sz(i, j, k) > 0){ + x = i * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + + Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; + } + } + } + } + } +} \ No newline at end of file From 418dee3d86bf7719541595ef75b5bdfa73861342 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 13:30:54 +0200 Subject: [PATCH 04/35] style adjustments --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 167 +++++++----------- .../FiniteDifferenceSolver/EvolveB.cpp | 6 +- .../FiniteDifferenceSolver/EvolveE.cpp | 6 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 14 +- Source/WarpX.H | 16 +- Source/WarpX.cpp | 16 +- 6 files changed, 90 insertions(+), 135 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 440f80bc794..abe53be3907 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -23,60 +23,45 @@ WarpX::InitEB () } void -WarpX::ComputeEdgeLengths() { +WarpX::ComputeEdgeLengths () { BL_PROFILE("ComputeEdgeLengths"); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); auto const &edge_centroid = eb_fact.getEdgeCent(); -// for ( amrex::MFIter mfi(*Bfield_fp[maxLevel()][0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){ amrex::Box const &box = mfi.validbox(); - amrex::FabType fab_type = flags[mfi].getType(box); - if (fab_type == amrex::FabType::regular) { - // every cell in box is all regular - } else if (fab_type == amrex::FabType::covered) { - // every cell in box is all covered - } else { - // mix of regular, cut and covered - // Note that edge_centroid[idim] is a MultiCutFab that does not have - // data if the fab type regular or covered. - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge = edge_centroid[idim]->const_array(mfi); - auto const &edge_lengths_dim = edge_lengths[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), - [=](int i, int j, int k) { - if (edge(i, j, k) == amrex::Real(-1.0)) { - // This edge is all covered - edge_lengths_dim(i, j, k) = 0.; - } else if (edge(i, j, k) == amrex::Real(1.0)) { - // This edge is all open - edge_lengths_dim(i, j, k) = 1.; - } else { - // This edge is cut. - amrex::Real edge_cent = edge(i, j, k); // edge centroid: (-0.5,0.5) - if (edge_cent < amrex::Real(0.0)) { - // The right side is covered. - edge_lengths_dim(i, j, k) = - amrex::Real(2.0) * edge_cent + amrex::Real(0.5) + 0.5; // (0, 1) - } else { - // The left side is covered - // TODO: rewrite well - edge_lengths_dim(i, j, k) = 1 - - (amrex::Real(2.0) * edge_cent - amrex::Real(0.5) + 0.5); // (0, 1) - } - } - }); - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){ + auto const &edge = edge_centroid[idim]->const_array(mfi); + auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), + [=](int i, int j, int k) { + if (edge(i, j, k) == amrex::Real(-1.0)) { + // This edge is all covered + edge_lengths_dim(i, j, k) = 0.; + } else if (edge(i, j, k) == amrex::Real(1.0)) { + // This edge is all open + edge_lengths_dim(i, j, k) = 1.; + } else { + // This edge is cut. + amrex::Real edge_cent = edge(i, j, k); // edge centroid: (-0.5,0.5) + if (edge_cent < amrex::Real(0.0)) { + // The right side is covered. + edge_lengths_dim(i, j, k) = amrex::Real(2.0) * edge_cent + amrex::Real(0.5) + 0.5; // (0, 1) + } else { + // The left side is covered + edge_lengths_dim(i, j, k) = 0.5 - amrex::Real(2.0) * edge_cent - amrex::Real(0.5); // (0, 1) + } + } + }); } - } - } + void -WarpX::ComputeFaceAreas() { +WarpX::ComputeFaceAreas () { BL_PROFILE("ComputeFaceAreas"); auto const eb_fact = fieldEBFactory(maxLevel()); @@ -86,35 +71,21 @@ WarpX::ComputeFaceAreas() { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); - amrex::FabType fab_type = flags[mfi].getType(box); - if (fab_type == amrex::FabType::regular) { - // every cell in box is all regular - } else if (fab_type == amrex::FabType::covered) { - // every cell in box is all covered - } else { - // mix of regular, cut and covered - // Note that edge_centroid[idim] is a MultiCutFab that does not have - // data if the fab type regular or covered. - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &face = area_frac[idim]->const_array(mfi); - auto const &face_areas_dim = face_areas[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), - [=](int i, int j, int k) { - face_areas_dim(i, j, k) = face(i, j, k); - }); - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = face(i, j, k); + }); } - } - } void -WarpX::ScaleEdges() { +WarpX::ScaleEdges () { auto const& cell_size = CellSize(maxLevel()); - auto const eb_fact = fieldEBFactory(maxLevel()); - auto const &flags = eb_fact.getMultiEBCellFlagFab(); auto const &edge_centroid = eb_fact.getEdgeCent(); @@ -122,28 +93,19 @@ WarpX::ScaleEdges() { amrex::Box const &box = mfi.validbox(); const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); - amrex::FabType fab_type = flags[mfi].getType(box); - if (fab_type == amrex::FabType::regular) { - // every cell in box is all regular - } else if (fab_type == amrex::FabType::covered) { - // every cell in box is all covered - } else { - // mix of regular, cut and covered - // Note that edge_centroid[idim] is a MultiCutFab that does not have - // data if the fab type regular or covered. - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge = edge_centroid[1]->const_array(mfi); - auto const &edge_lengths_dim = edge_lengths[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; - } - } - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const &edge = edge_centroid[1]->const_array(mfi); + auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; + } + } } } } + void WarpX::ScaleAreas() { auto const& cell_size = CellSize(maxLevel()); @@ -156,33 +118,24 @@ WarpX::ScaleAreas() { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); - amrex::FabType fab_type = flags[mfi].getType(box); const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); - if (fab_type == amrex::FabType::regular) { - // every cell in box is all regular - } else if (fab_type == amrex::FabType::covered) { - // every cell in box is all covered - } else { - // mix of regular, cut and covered - // Note that edge_centroid[idim] is a MultiCutFab that does not have - // data if the fab type regular or covered. - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (idim == 0) { - full_area = cell_size[1]*cell_size[2]; - } else if (idim == 1) { - full_area = cell_size[0]*cell_size[2]; - } else { - full_area = cell_size[0]*cell_size[1]; - } - auto const &face = area_frac[idim]->const_array(mfi); - auto const &face_areas_dim = face_areas[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - face_areas_dim(i,j,hi.z/2) *= full_area; - } - } - } + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (idim == 0) { + full_area = cell_size[1]*cell_size[2]; + } else if (idim == 1) { + full_area = cell_size[0]*cell_size[2]; + } else { + full_area = cell_size[0]*cell_size[1]; + } + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + face_areas_dim(i,j,hi.z/2) *= full_area; + } + } } } } + diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index b68450a37ef..670effb2fc4 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -87,9 +87,9 @@ void FiniteDifferenceSolver::EvolveBCartesian ( Array4 const& Ez = Efield[2]->array(mfi); #ifdef AMREX_USE_EB - Array4 const& Sx = face_areas[0]->array(mfi); - Array4 const& Sy = face_areas[1]->array(mfi); - Array4 const& Sz = face_areas[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); #endif // Extract stencil coefficients diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 52ad487caf0..77d431c9f81 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -97,9 +97,9 @@ void FiniteDifferenceSolver::EvolveECartesian ( Array4 const& jz = Jfield[2]->array(mfi); #ifdef AMREX_USE_EB - Array4 const& lx = edge_lengths[0]->array(mfi); - Array4 const& ly = edge_lengths[1]->array(mfi); - Array4 const& lz = edge_lengths[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); #endif // Extract stencil coefficients diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1960a5471b2..1a898c1e673 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -192,9 +192,9 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) // Evolve B field in regular cells if (patch_type == PatchType::fine) { - m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], face_areas[lev], lev, a_dt ); + m_fdtd_solver_fp[lev]->EvolveB(Bfield_fp[lev], Efield_fp[lev], m_face_areas[lev], lev, a_dt ); } else { - m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], face_areas[lev], lev, a_dt ); + m_fdtd_solver_cp[lev]->EvolveB(Bfield_cp[lev], Efield_cp[lev], m_face_areas[lev], lev, a_dt ); } // Evolve B field in PML cells @@ -242,11 +242,13 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { // Evolve E field in regular cells if (patch_type == PatchType::fine) { - m_fdtd_solver_fp[lev]->EvolveE( Efield_fp[lev], Bfield_fp[lev], - current_fp[lev], edge_lengths[lev], F_fp[lev], lev, a_dt ); + m_fdtd_solver_fp[lev]->EvolveE(Efield_fp[lev], Bfield_fp[lev], + current_fp[lev], m_edge_lengths[lev], + F_fp[lev], lev, a_dt ); } else { - m_fdtd_solver_cp[lev]->EvolveE( Efield_cp[lev], Bfield_cp[lev], - current_cp[lev], edge_lengths[lev], F_cp[lev], lev, a_dt ); + m_fdtd_solver_cp[lev]->EvolveE(Efield_cp[lev], Bfield_cp[lev], + current_cp[lev], m_edge_lengths[lev], + F_cp[lev], lev, a_dt ); } // Evolve E field in PML cells diff --git a/Source/WarpX.H b/Source/WarpX.H index 3ad6ca0051d..4f39e26b251 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -869,8 +869,8 @@ private: #ifdef AMREX_USE_EB //EB grid info - amrex::Vector, 3 > > edge_lengths; - amrex::Vector, 3 > > face_areas; + amrex::Vector, 3 > > m_edge_lengths; + amrex::Vector, 3 > > m_face_areas; #endif // store fine patch @@ -1033,13 +1033,13 @@ private: void InitEB (); #ifdef AMREX_USE_EB - void ComputeEdgeLengths(); - void ComputeFaceAreas(); - void ScaleEdges(); - void ScaleAreas(); + void ComputeEdgeLengths (); + void ComputeFaceAreas (); + void ScaleEdges (); + void ScaleAreas (); - static amrex::RealArray AnalyticSolSphere(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere); - void InitEMFieldSphere(); + static amrex::RealArray AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere); + void InitEMFieldSphere (); #endif private: diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c10634381ff..83917a24bc2 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -230,8 +230,8 @@ WarpX::WarpX () Bfield_avg_fp.resize(nlevs_max); #ifdef AMREX_USE_EB - edge_lengths.resize(nlevs_max); - face_areas.resize(nlevs_max); + m_edge_lengths.resize(nlevs_max); + m_face_areas.resize(nlevs_max); #endif current_store.resize(nlevs_max); @@ -1254,12 +1254,12 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // EB info are needed only at the coarsest level if (lev == maxLevel()) { - edge_lengths[lev][0] = std::make_unique(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[x]")); - edge_lengths[lev][1] = std::make_unique(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[y]")); - edge_lengths[lev][2] = std::make_unique(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE,tag("edge_lengths[z]")); - face_areas[lev][0] = std::make_unique(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE,tag("face_areas[x]")); - face_areas[lev][1] = std::make_unique(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE,tag("face_areas[y]")); - face_areas[lev][2] = std::make_unique(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE,tag("face_areas[z]")); + m_edge_lengths[lev][0] = std::make_unique(amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[x]")); + m_edge_lengths[lev][1] = std::make_unique(amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[y]")); + m_edge_lengths[lev][2] = std::make_unique(amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[z]")); + m_face_areas[lev][0] = std::make_unique(amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[x]")); + m_face_areas[lev][1] = std::make_unique(amrex::convert(ba, By_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[y]")); + m_face_areas[lev][2] = std::make_unique(amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[z]")); } #endif From 01742abc32ee973241d083e35736c1bac4aa6cde Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 13:47:14 +0200 Subject: [PATCH 05/35] fixing tabs --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 52 +++++------ Source/EmbeddedBoundary/WarpXInitEMField.cpp | 90 ++++++++++---------- Source/Evolve/WarpXEvolve.cpp | 18 ++-- Source/FieldSolver/WarpXPushFieldsEM.cpp | 8 +- Source/WarpX.cpp | 12 +-- 5 files changed, 90 insertions(+), 90 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index abe53be3907..2b57b7cdd89 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -72,12 +72,12 @@ WarpX::ComputeFaceAreas () { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &face = area_frac[idim]->const_array(mfi); - auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), - [=](int i, int j, int k) { + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + [=](int i, int j, int k) { face_areas_dim(i, j, k) = face(i, j, k); - }); + }); } } } @@ -94,13 +94,13 @@ WarpX::ScaleEdges () { const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge = edge_centroid[1]->const_array(mfi); - auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; - } - } + auto const &edge = edge_centroid[1]->const_array(mfi); + auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; + } + } } } } @@ -121,20 +121,20 @@ WarpX::ScaleAreas() { const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (idim == 0) { - full_area = cell_size[1]*cell_size[2]; - } else if (idim == 1) { - full_area = cell_size[0]*cell_size[2]; - } else { - full_area = cell_size[0]*cell_size[1]; - } - auto const &face = area_frac[idim]->const_array(mfi); - auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - face_areas_dim(i,j,hi.z/2) *= full_area; - } - } + if (idim == 0) { + full_area = cell_size[1]*cell_size[2]; + } else if (idim == 1) { + full_area = cell_size[0]*cell_size[2]; + } else { + full_area = cell_size[0]*cell_size[1]; + } + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + for(int i = lo.x; i <= hi.x; ++i) { + for(int j = lo.y; j <= hi.y; ++j) { + face_areas_dim(i,j,hi.z/2) *= full_area; + } + } } } } diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index febf41ac173..408e511220b 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -23,54 +23,54 @@ WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Re void WarpX::InitEMFieldSphere (){ - amrex::Real r_sphere; - amrex::ParmParse pp_eb2("eb2"); - amrex::ParmParse pp_geometry("geometry"); - pp_eb2.get("sphere_radius", r_sphere); - amrex::RealArray prob_lo; - pp_geometry.get("prob_lo", prob_lo); - auto &Bfield = Bfield_fp[maxLevel()]; - auto &Efield = Efield_fp[maxLevel()]; + amrex::Real r_sphere; + amrex::ParmParse pp_eb2("eb2"); + amrex::ParmParse pp_geometry("geometry"); + pp_eb2.get("sphere_radius", r_sphere); + amrex::RealArray prob_lo; + pp_geometry.get("prob_lo", prob_lo); + auto &Bfield = Bfield_fp[maxLevel()]; + auto &Efield = Efield_fp[maxLevel()]; - auto &face_areas_0 = m_face_areas[maxLevel()]; - for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - amrex::Array4 const &Bx = Bfield[0]->array(mfi); - amrex::Array4 const &By = Bfield[1]->array(mfi); - amrex::Array4 const &Bz = Bfield[2]->array(mfi); - amrex::Array4 const &Sx = face_areas_0[0]->array(mfi); - amrex::Array4 const &Sy = face_areas_0[1]->array(mfi); - amrex::Array4 const &Sz = face_areas_0[2]->array(mfi); - amrex::Box const &box = mfi.validbox(); - const auto hi = amrex::ubound(box); - amrex::Real x, y, z; - for (int i = 0; i < hi.x; i++) { - for (int j = 0; j < hi.y; j++) { - for (int k = 0; k < hi.z; k++) { - if (Sx(i, j, k) > 0) { - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = k * CellSize(maxLevel())[2] + prob_lo[2]; + auto &face_areas_0 = m_face_areas[maxLevel()]; + for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + amrex::Array4 const &Bx = Bfield[0]->array(mfi); + amrex::Array4 const &By = Bfield[1]->array(mfi); + amrex::Array4 const &Bz = Bfield[2]->array(mfi); + amrex::Array4 const &Sx = face_areas_0[0]->array(mfi); + amrex::Array4 const &Sy = face_areas_0[1]->array(mfi); + amrex::Array4 const &Sz = face_areas_0[2]->array(mfi); + amrex::Box const &box = mfi.validbox(); + const auto hi = amrex::ubound(box); + amrex::Real x, y, z; + for (int i = 0; i < hi.x; i++) { + for (int j = 0; j < hi.y; j++) { + for (int k = 0; k < hi.z; k++) { + if (Sx(i, j, k) > 0) { + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = k * CellSize(maxLevel())[2] + prob_lo[2]; - Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; - } + Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; + } - if (Sy(i, j, k) > 0) { - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = j * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + if (Sy(i, j, k) > 0) { + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = j * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; - } + By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; + } - if (Sz(i, j, k) > 0){ - x = i * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + if (Sz(i, j, k) > 0){ + x = i * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; - } - } - } - } - } -} \ No newline at end of file + Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; + } + } + } + } + } +} diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 3f09de5fad8..c680876bc03 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -46,15 +46,15 @@ WarpX::Evolve (int numsteps) Real walltime, walltime_start = amrex::second(); #ifdef AMREX_USE_EB - //If necessary initialize the em fields - std::string geom_type; - amrex::ParmParse pp_eb2("eb2"); - bool init_em_field; - pp_eb2.get("geom_type", geom_type); - pp_eb2.get("init_em_field", init_em_field); - if(init_em_field and geom_type=="sphere") { - InitEMFieldSphere(); - } + //If necessary initialize the em fields + std::string geom_type; + amrex::ParmParse pp_eb2("eb2"); + bool init_em_field; + pp_eb2.get("geom_type", geom_type); + pp_eb2.get("init_em_field", init_em_field); + if(init_em_field and geom_type=="sphere") { + InitEMFieldSphere(); + } #endif for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1a898c1e673..a991c539dc7 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -243,12 +243,12 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) // Evolve E field in regular cells if (patch_type == PatchType::fine) { m_fdtd_solver_fp[lev]->EvolveE(Efield_fp[lev], Bfield_fp[lev], - current_fp[lev], m_edge_lengths[lev], - F_fp[lev], lev, a_dt ); + current_fp[lev], m_edge_lengths[lev], + F_fp[lev], lev, a_dt ); } else { m_fdtd_solver_cp[lev]->EvolveE(Efield_cp[lev], Bfield_cp[lev], - current_cp[lev], m_edge_lengths[lev], - F_cp[lev], lev, a_dt ); + current_cp[lev], m_edge_lengths[lev], + F_cp[lev], lev, a_dt ); } // Evolve E field in PML cells diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 83917a24bc2..bae2d32d345 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1254,12 +1254,12 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // EB info are needed only at the coarsest level if (lev == maxLevel()) { - m_edge_lengths[lev][0] = std::make_unique(amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[x]")); - m_edge_lengths[lev][1] = std::make_unique(amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[y]")); - m_edge_lengths[lev][2] = std::make_unique(amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[z]")); - m_face_areas[lev][0] = std::make_unique(amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[x]")); - m_face_areas[lev][1] = std::make_unique(amrex::convert(ba, By_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[y]")); - m_face_areas[lev][2] = std::make_unique(amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[z]")); + m_edge_lengths[lev][0] = std::make_unique(amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[x]")); + m_edge_lengths[lev][1] = std::make_unique(amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[y]")); + m_edge_lengths[lev][2] = std::make_unique(amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngE, tag("m_edge_lengths[z]")); + m_face_areas[lev][0] = std::make_unique(amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[x]")); + m_face_areas[lev][1] = std::make_unique(amrex::convert(ba, By_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[y]")); + m_face_areas[lev][2] = std::make_unique(amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngE, tag("m_face_areas[z]")); } #endif From 13072662a282d3ea447206aa1cdce470ceeb6b06 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 13:50:13 +0200 Subject: [PATCH 06/35] fixed name of analysis script --- .../plots/analysis_fields.py | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py diff --git a/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py b/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py new file mode 100644 index 00000000000..7228456a055 --- /dev/null +++ b/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py @@ -0,0 +1,58 @@ +import yt +import os +from scipy.constants import mu_0 +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import AxesGrid + +fig = plt.figure(1, (1,15)) + +grid = AxesGrid(fig, (0.075,0.075,0.85,0.85), + nrows_ncols = (1, 3), + axes_pad = 1.0, + label_mode = "1", + share_all = True, + cbar_location="right", + cbar_mode="each", + cbar_size="3%", + cbar_pad="0%") + + +path, dirs, files = next(os.walk("../diags/")) +N_dirs = len(dirs) + +for i in range(0, N_dirs): + filename = '../diags/diag'+str(int(i)).zfill(5)+'/' + print(filename) + ds = yt.load(filename) + comp = 'Ez' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -5000, zmax = 5000) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[0].axes + plot.cax = grid.cbar_axes[0] + p._setup_plots() + + comp = 'Bx' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[1].axes + plot.cax = grid.cbar_axes[1] + p._setup_plots() + + comp = 'By' + p = yt.SlicePlot( ds, 2, comp, origin='native' ) + p.set_log(comp, False) + p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) + plot = p.plots[comp] + plot.figure = fig + plot.axes = grid[2].axes + plot.cax = grid.cbar_axes[2] + p._setup_plots() + + fig.set_size_inches(13,5) + p.save(str(i)) From a1fbfaa30e3e37a2f6dc8355125ed88523813669 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 13:54:20 +0200 Subject: [PATCH 07/35] fixed name of analysis script --- .../plots/plot_fields.py | 58 ------------------- 1 file changed, 58 deletions(-) delete mode 100644 Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py diff --git a/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py b/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py deleted file mode 100644 index 7228456a055..00000000000 --- a/Examples/Modules/embedded_boundary_sphere/plots/plot_fields.py +++ /dev/null @@ -1,58 +0,0 @@ -import yt -import os -from scipy.constants import mu_0 -import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import AxesGrid - -fig = plt.figure(1, (1,15)) - -grid = AxesGrid(fig, (0.075,0.075,0.85,0.85), - nrows_ncols = (1, 3), - axes_pad = 1.0, - label_mode = "1", - share_all = True, - cbar_location="right", - cbar_mode="each", - cbar_size="3%", - cbar_pad="0%") - - -path, dirs, files = next(os.walk("../diags/")) -N_dirs = len(dirs) - -for i in range(0, N_dirs): - filename = '../diags/diag'+str(int(i)).zfill(5)+'/' - print(filename) - ds = yt.load(filename) - comp = 'Ez' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -5000, zmax = 5000) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[0].axes - plot.cax = grid.cbar_axes[0] - p._setup_plots() - - comp = 'Bx' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[1].axes - plot.cax = grid.cbar_axes[1] - p._setup_plots() - - comp = 'By' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[2].axes - plot.cax = grid.cbar_axes[2] - p._setup_plots() - - fig.set_size_inches(13,5) - p.save(str(i)) From 417e137ef7516404dc34ac7704a1371c93816b0e Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 15:13:36 +0200 Subject: [PATCH 08/35] fixed a few wrong preprocessor directives --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 8 ++++++++ Source/WarpX.H | 4 ---- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 2b57b7cdd89..cc66fdf9bec 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -24,6 +24,7 @@ WarpX::InitEB () void WarpX::ComputeEdgeLengths () { +#ifdef AMREX_USE_EB BL_PROFILE("ComputeEdgeLengths"); auto const eb_fact = fieldEBFactory(maxLevel()); @@ -57,11 +58,13 @@ WarpX::ComputeEdgeLengths () { }); } } +#endif } void WarpX::ComputeFaceAreas () { +#ifdef AMREX_USE_EB BL_PROFILE("ComputeFaceAreas"); auto const eb_fact = fieldEBFactory(maxLevel()); @@ -80,10 +83,12 @@ WarpX::ComputeFaceAreas () { }); } } +#endif } void WarpX::ScaleEdges () { +#ifdef AMREX_USE_EB auto const& cell_size = CellSize(maxLevel()); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); @@ -103,11 +108,13 @@ WarpX::ScaleEdges () { } } } +#endif } void WarpX::ScaleAreas() { +#ifdef AMREX_USE_EB auto const& cell_size = CellSize(maxLevel()); amrex::Real full_area; @@ -137,5 +144,6 @@ WarpX::ScaleAreas() { } } } +#endif } diff --git a/Source/WarpX.H b/Source/WarpX.H index 4f39e26b251..252d64b3f2a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -867,11 +867,9 @@ private: amrex::Vector, 3 > > Efield_avg_fp; amrex::Vector, 3 > > Bfield_avg_fp; -#ifdef AMREX_USE_EB //EB grid info amrex::Vector, 3 > > m_edge_lengths; amrex::Vector, 3 > > m_face_areas; -#endif // store fine patch amrex::Vector, 3 > > current_store; @@ -1032,7 +1030,6 @@ private: void InitEB (); -#ifdef AMREX_USE_EB void ComputeEdgeLengths (); void ComputeFaceAreas (); void ScaleEdges (); @@ -1040,7 +1037,6 @@ private: static amrex::RealArray AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere); void InitEMFieldSphere (); -#endif private: // void EvolvePSATD (int numsteps); From 601f9eb2ec6f8c2100304379b2bea1c6cf9d1851 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 15:30:43 +0200 Subject: [PATCH 09/35] workaround for missing boost --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 31 ++++++++++++-------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index 408e511220b..adb424ad49c 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -1,29 +1,35 @@ #include "WarpX.H" + +#ifdef AMREX_USE_EB #include +#endif amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { - amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = PhysConst::mu0; +#ifdef AMREX_USE_EB + amrex::Real k=2.7437/r_sphere; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = PhysConst::mu0; - amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); + amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; +#endif } void WarpX::InitEMFieldSphere (){ - amrex::Real r_sphere; +#ifdef AMREX_USE_EB + amrex::Real r_sphere; amrex::ParmParse pp_eb2("eb2"); amrex::ParmParse pp_geometry("geometry"); pp_eb2.get("sphere_radius", r_sphere); @@ -73,4 +79,5 @@ WarpX::InitEMFieldSphere (){ } } } +#endif } From bcab86f576dc8bc7c8bfcc8f5447d69fe80d1af6 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 15:43:17 +0200 Subject: [PATCH 10/35] Revert "workaround for missing boost" This reverts commit 601f9eb2ec6f8c2100304379b2bea1c6cf9d1851. --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 31 ++++++++------------ 1 file changed, 12 insertions(+), 19 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index adb424ad49c..408e511220b 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -1,35 +1,29 @@ #include "WarpX.H" - -#ifdef AMREX_USE_EB #include -#endif amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { -#ifdef AMREX_USE_EB - amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = PhysConst::mu0; + amrex::Real k=2.7437/r_sphere; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = PhysConst::mu0; - amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); + amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; -#endif } void WarpX::InitEMFieldSphere (){ -#ifdef AMREX_USE_EB - amrex::Real r_sphere; + amrex::Real r_sphere; amrex::ParmParse pp_eb2("eb2"); amrex::ParmParse pp_geometry("geometry"); pp_eb2.get("sphere_radius", r_sphere); @@ -79,5 +73,4 @@ WarpX::InitEMFieldSphere (){ } } } -#endif } From 9dd2bd22e6e48a643f65e18d7c14680e01b39cb7 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 7 Apr 2021 15:45:54 +0200 Subject: [PATCH 11/35] another workaround for missing boost --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 37 ++++++++++++-------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index 408e511220b..5911be3208b 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -1,29 +1,35 @@ #include "WarpX.H" + +#ifdef AMREX_USE_EB #include +#endif amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { - amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = PhysConst::mu0; - - amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + amrex::Real k=2.7437/r_sphere; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = PhysConst::mu0; +#ifdef AMREX_USE_EB + amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); +#else + amrex::Real H_phi = 0; +#endif + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); - return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; + return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; } void WarpX::InitEMFieldSphere (){ - amrex::Real r_sphere; +#ifdef AMREX_USE_EB + amrex::Real r_sphere; amrex::ParmParse pp_eb2("eb2"); amrex::ParmParse pp_geometry("geometry"); pp_eb2.get("sphere_radius", r_sphere); @@ -73,4 +79,5 @@ WarpX::InitEMFieldSphere (){ } } } +#endif } From a53fdba0206046aaa19bc67b0ba6ed4834170d52 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 10:50:43 +0200 Subject: [PATCH 12/35] getting rid of boost by depending on c++17 --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index 5911be3208b..4e065beff83 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -1,9 +1,6 @@ #include "WarpX.H" - -#ifdef AMREX_USE_EB -#include -#endif +#include amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { @@ -14,11 +11,8 @@ WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Re amrex::Real H_r = 0; amrex::Real H_theta = 0; amrex::Real mu_r = PhysConst::mu0; -#ifdef AMREX_USE_EB - amrex::Real H_phi = k / (r_sphere * mu_r) * boost::math::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); -#else - amrex::Real H_phi = 0; -#endif + amrex::Real H_phi = k / (r_sphere * mu_r) * std::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); From 00f442b41de1e4a15df108228673cc7347303091 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 10:55:48 +0200 Subject: [PATCH 13/35] Removed a few unused variables --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 4 ---- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 1 - 2 files changed, 5 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index cc66fdf9bec..e2acdfe8786 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -92,14 +92,12 @@ WarpX::ScaleEdges () { auto const& cell_size = CellSize(maxLevel()); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); - auto const &edge_centroid = eb_fact.getEdgeCent(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge = edge_centroid[1]->const_array(mfi); auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); for(int i = lo.x; i <= hi.x; ++i) { for(int j = lo.y; j <= hi.y; ++j) { @@ -121,7 +119,6 @@ WarpX::ScaleAreas() { auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); - auto const &area_frac = eb_fact.getAreaFrac(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); @@ -135,7 +132,6 @@ WarpX::ScaleAreas() { } else { full_area = cell_size[0]*cell_size[1]; } - auto const &face = area_frac[idim]->const_array(mfi); auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); for(int i = lo.x; i <= hi.x; ++i) { for(int j = lo.y; j <= hi.y; ++j) { diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index 4e065beff83..ebd1f4a1caa 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -30,7 +30,6 @@ WarpX::InitEMFieldSphere (){ amrex::RealArray prob_lo; pp_geometry.get("prob_lo", prob_lo); auto &Bfield = Bfield_fp[maxLevel()]; - auto &Efield = Efield_fp[maxLevel()]; auto &face_areas_0 = m_face_areas[maxLevel()]; for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { From b7971ee5a77a0723bb6a3524ff78683e22e7a7c9 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 10:57:13 +0200 Subject: [PATCH 14/35] adding USE_EB to addToCompileString for EB testing --- Regression/WarpX-tests.ini | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 6bd1ce58069..c25796eeedb 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -1915,3 +1915,19 @@ compileTest = 0 doVis = 0 compareParticles = 0 analysisRoutine = Examples/Modules/resampling/analysis_leveling_thinning.py + +[embedded_boundary_sphere] +buildDir = . +inputFile = Examples/Modules/embedded_boundary_sphere +runtime_params = +dim = 3 +addToCompileString = USE_EB=TRUE +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py \ No newline at end of file From a60de3cd78bb37a087068570288310f7914beca5 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 11:01:02 +0200 Subject: [PATCH 15/35] removed tabs --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 22 ++++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index ebd1f4a1caa..bcc46b26277 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -4,18 +4,18 @@ amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { - amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = PhysConst::mu0; + amrex::Real k=2.7437/r_sphere; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = PhysConst::mu0; amrex::Real H_phi = k / (r_sphere * mu_r) * std::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; } @@ -23,7 +23,7 @@ WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Re void WarpX::InitEMFieldSphere (){ #ifdef AMREX_USE_EB - amrex::Real r_sphere; + amrex::Real r_sphere; amrex::ParmParse pp_eb2("eb2"); amrex::ParmParse pp_geometry("geometry"); pp_eb2.get("sphere_radius", r_sphere); From b28b862729117dd3b52391522665f1a191e38235 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 11:14:17 +0200 Subject: [PATCH 16/35] fixing the inputs name for EB sphere test --- Regression/WarpX-tests.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c25796eeedb..e605e1a87f8 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -1918,7 +1918,7 @@ analysisRoutine = Examples/Modules/resampling/analysis_leveling_thinning.py [embedded_boundary_sphere] buildDir = . -inputFile = Examples/Modules/embedded_boundary_sphere +inputFile = Examples/Modules/embedded_boundary_sphere/inputs_3d runtime_params = dim = 3 addToCompileString = USE_EB=TRUE @@ -1929,5 +1929,5 @@ useOMP = 1 numthreads = 1 compileTest = 0 doVis = 0 -compareParticles = 0 +compareParticles = 0s analysisRoutine = Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py \ No newline at end of file From 7e129a62bfac5224b657eef38b5f326b1408a140 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 11:40:36 +0200 Subject: [PATCH 17/35] shortened the test --- Examples/Modules/embedded_boundary_sphere/inputs_3d | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Modules/embedded_boundary_sphere/inputs_3d b/Examples/Modules/embedded_boundary_sphere/inputs_3d index a9513b62ce6..975b6468bb9 100644 --- a/Examples/Modules/embedded_boundary_sphere/inputs_3d +++ b/Examples/Modules/embedded_boundary_sphere/inputs_3d @@ -1,4 +1,4 @@ -max_step = 1000 +max_step = 100 amr.n_cell = 48 48 48 amr.max_grid_size = 128 amr.max_level = 0 From 66fac2168968fd4435132b9865674aa94dcb5390 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 8 Apr 2021 11:41:28 +0200 Subject: [PATCH 18/35] zero padding the names of the images --- .../Modules/embedded_boundary_sphere/plots/analysis_fields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py b/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py index 7228456a055..284d95affea 100644 --- a/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py @@ -55,4 +55,4 @@ p._setup_plots() fig.set_size_inches(13,5) - p.save(str(i)) + p.save(str(i).zfill(4)) From e8ea5e1da14f0fe164201df969b2e734ee409fde Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 9 Apr 2021 10:45:21 +0200 Subject: [PATCH 19/35] adjusted two for loops --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 50 ++++++++++++------------- 1 file changed, 24 insertions(+), 26 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index e2acdfe8786..03f6c12fcff 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -37,7 +37,7 @@ WarpX::ComputeEdgeLengths () { auto const &edge = edge_centroid[idim]->const_array(mfi); auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), - [=](int i, int j, int k) { + [=](int i, int j, int k) { if (edge(i, j, k) == amrex::Real(-1.0)) { // This edge is all covered edge_lengths_dim(i, j, k) = 0.; @@ -68,7 +68,6 @@ WarpX::ComputeFaceAreas () { BL_PROFILE("ComputeFaceAreas"); auto const eb_fact = fieldEBFactory(maxLevel()); - auto const &flags = eb_fact.getMultiEBCellFlagFab(); auto const &area_frac = eb_fact.getAreaFrac(); @@ -78,8 +77,8 @@ WarpX::ComputeFaceAreas () { auto const &face = area_frac[idim]->const_array(mfi); auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), - [=](int i, int j, int k) { - face_areas_dim(i, j, k) = face(i, j, k); + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = face(i, j, k); }); } } @@ -92,18 +91,17 @@ WarpX::ScaleEdges () { auto const& cell_size = CellSize(maxLevel()); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &edge_centroid = eb_fact.getEdgeCent(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); - const auto lo = amrex::lbound(box); - const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - edge_lengths_dim(i,j,hi.z/2) *= cell_size[idim]; - } - } + auto const &edge = edge_centroid[idim]->const_array(mfi); + auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), + [=](int i, int j, int k) { + edge_lengths_dim(i, j, k) *= cell_size[idim]; + }); } } #endif @@ -117,29 +115,29 @@ WarpX::ScaleAreas() { amrex::Real full_area; auto const eb_fact = fieldEBFactory(maxLevel()); - auto const &flags = eb_fact.getMultiEBCellFlagFab(); + auto const &area_frac = eb_fact.getAreaFrac(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); const auto lo = amrex::lbound(box); const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (idim == 0) { + if (idim == 0) { full_area = cell_size[1]*cell_size[2]; - } else if (idim == 1) { - full_area = cell_size[0]*cell_size[2]; - } else { - full_area = cell_size[0]*cell_size[1]; - } - auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); - for(int i = lo.x; i <= hi.x; ++i) { - for(int j = lo.y; j <= hi.y; ++j) { - face_areas_dim(i,j,hi.z/2) *= full_area; - } - } + } else if (idim == 1) { + full_area = cell_size[0]*cell_size[2]; + } else { + full_area = cell_size[0]*cell_size[1]; + } + auto const &face = area_frac[idim]->const_array(mfi); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) *= full_area; + }); + } } #endif } - From 983e6fb79e6778cf7e18afd8024038703cbaa7f9 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 9 Apr 2021 11:16:28 +0200 Subject: [PATCH 20/35] removed some unused variables --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 03f6c12fcff..d43bac4df17 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -120,8 +120,6 @@ WarpX::ScaleAreas() { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); - const auto lo = amrex::lbound(box); - const auto hi = amrex::ubound(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if (idim == 0) { full_area = cell_size[1]*cell_size[2]; From 49cfb370f06ae5d6bce12c39d2330a2dde65c0c0 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 21 Apr 2021 15:25:07 +0200 Subject: [PATCH 21/35] improving the fields initialization --- Source/EmbeddedBoundary/WarpXInitEMField.cpp | 88 ++++++++++---------- 1 file changed, 42 insertions(+), 46 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp index bcc46b26277..1fb53fce1e2 100644 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEMField.cpp @@ -5,34 +5,34 @@ amrex::RealArray WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = PhysConst::mu0; + amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); + amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); + amrex::Real phi = atan2(y, x); + amrex::Real H_r = 0; + amrex::Real H_theta = 0; + amrex::Real mu_r = 1; //PhysConst::mu0; amrex::Real H_phi = k / (r_sphere * mu_r) * std::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); + amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); + amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); + amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); - return {PhysConst::mu0*H_x/1e6, PhysConst::mu0*H_y/1e6, PhysConst::mu0*H_z/1e6}; + return {PhysConst::mu0*H_x, PhysConst::mu0*H_y, PhysConst::mu0*H_z}; } void WarpX::InitEMFieldSphere (){ #ifdef AMREX_USE_EB - amrex::Real r_sphere; - amrex::ParmParse pp_eb2("eb2"); - amrex::ParmParse pp_geometry("geometry"); - pp_eb2.get("sphere_radius", r_sphere); - amrex::RealArray prob_lo; - pp_geometry.get("prob_lo", prob_lo); - auto &Bfield = Bfield_fp[maxLevel()]; + amrex::Real r_sphere; + amrex::ParmParse pp_eb2("eb2"); + amrex::ParmParse pp_geometry("geometry"); + pp_eb2.get("sphere_radius", r_sphere); + amrex::RealArray prob_lo; + pp_geometry.get("prob_lo", prob_lo); + auto &Bfield = Bfield_fp[maxLevel()]; - auto &face_areas_0 = m_face_areas[maxLevel()]; - for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto &face_areas_0 = m_face_areas[maxLevel()]; + for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { amrex::Array4 const &Bx = Bfield[0]->array(mfi); amrex::Array4 const &By = Bfield[1]->array(mfi); amrex::Array4 const &Bz = Bfield[2]->array(mfi); @@ -40,37 +40,33 @@ WarpX::InitEMFieldSphere (){ amrex::Array4 const &Sy = face_areas_0[1]->array(mfi); amrex::Array4 const &Sz = face_areas_0[2]->array(mfi); amrex::Box const &box = mfi.validbox(); - const auto hi = amrex::ubound(box); - amrex::Real x, y, z; - for (int i = 0; i < hi.x; i++) { - for (int j = 0; j < hi.y; j++) { - for (int k = 0; k < hi.z; k++) { - if (Sx(i, j, k) > 0) { - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = k * CellSize(maxLevel())[2] + prob_lo[2]; - Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; - } + amrex::LoopOnCpu(box, [=] (int i, int j, int k) { + amrex::Real x, y, z; + if (Sx(i, j, k) > 0) { + x = i * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - if (Sy(i, j, k) > 0) { - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = j * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; + } - By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; - } + if (Sy(i, j, k) > 0) { + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = j * CellSize(maxLevel())[1] + prob_lo[1]; + z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - if (Sz(i, j, k) > 0){ - x = i * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; + By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; + } - Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; - } - } - } - } - } + if (Sz(i, j, k) > 0){ + x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; + y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; + z = k * CellSize(maxLevel())[2] + prob_lo[2]; + + Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; + } + }); + } #endif } From b478097a667f9800c5f2164284967fb029412a77 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 22 Apr 2021 19:01:35 +0200 Subject: [PATCH 22/35] removed the sphere test and implemented the cube test --- .../embedded_boundary_sphere/inputs_3d | 23 ------ .../plots/analysis_fields.py | 58 --------------- .../embedded_boundary_cube.json | 10 +++ Regression/WarpX-tests.ini | 8 +-- Source/EmbeddedBoundary/CMakeLists.txt | 1 - Source/EmbeddedBoundary/WarpXInitEMField.cpp | 72 ------------------- Source/Evolve/WarpXEvolve.cpp | 12 ---- Source/WarpX.H | 3 - 8 files changed, 14 insertions(+), 173 deletions(-) delete mode 100644 Examples/Modules/embedded_boundary_sphere/inputs_3d delete mode 100644 Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py create mode 100644 Regression/Checksum/benchmarks_json/embedded_boundary_cube.json delete mode 100644 Source/EmbeddedBoundary/WarpXInitEMField.cpp diff --git a/Examples/Modules/embedded_boundary_sphere/inputs_3d b/Examples/Modules/embedded_boundary_sphere/inputs_3d deleted file mode 100644 index 975b6468bb9..00000000000 --- a/Examples/Modules/embedded_boundary_sphere/inputs_3d +++ /dev/null @@ -1,23 +0,0 @@ -max_step = 100 -amr.n_cell = 48 48 48 -amr.max_grid_size = 128 -amr.max_level = 0 - -geometry.coord_sys = 0 -geometry.is_periodic = 0 0 0 -geometry.prob_lo = -0.5 -0.5 -0.5 -geometry.prob_hi = 0.5 0.5 0.5 -warpx.do_pml = 0 -warpx.const_dt = 1e-6 -warpx.cfl = 1 - -eb2.geom_type = sphere -eb2.sphere_center = 0. 0. 0. -eb2.sphere_radius = 0.3 -eb2.sphere_has_fluid_inside = true -eb2.init_em_field = true - -diagnostics.diags_names = diag -diag.intervals = 1 -diag.diag_type = Full -diag.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho diff --git a/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py b/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py deleted file mode 100644 index 284d95affea..00000000000 --- a/Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py +++ /dev/null @@ -1,58 +0,0 @@ -import yt -import os -from scipy.constants import mu_0 -import matplotlib.pyplot as plt -from mpl_toolkits.axes_grid1 import AxesGrid - -fig = plt.figure(1, (1,15)) - -grid = AxesGrid(fig, (0.075,0.075,0.85,0.85), - nrows_ncols = (1, 3), - axes_pad = 1.0, - label_mode = "1", - share_all = True, - cbar_location="right", - cbar_mode="each", - cbar_size="3%", - cbar_pad="0%") - - -path, dirs, files = next(os.walk("../diags/")) -N_dirs = len(dirs) - -for i in range(0, N_dirs): - filename = '../diags/diag'+str(int(i)).zfill(5)+'/' - print(filename) - ds = yt.load(filename) - comp = 'Ez' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -5000, zmax = 5000) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[0].axes - plot.cax = grid.cbar_axes[0] - p._setup_plots() - - comp = 'Bx' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[1].axes - plot.cax = grid.cbar_axes[1] - p._setup_plots() - - comp = 'By' - p = yt.SlicePlot( ds, 2, comp, origin='native' ) - p.set_log(comp, False) - p.set_zlim(field = comp, zmin = -8*mu_0, zmax = 8*mu_0) - plot = p.plots[comp] - plot.figure = fig - plot.axes = grid[2].axes - plot.cax = grid.cbar_axes[2] - p._setup_plots() - - fig.set_size_inches(13,5) - p.save(str(i).zfill(4)) diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json new file mode 100644 index 00000000000..bf0f1e525fb --- /dev/null +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json @@ -0,0 +1,10 @@ +{ + "lev=0": { + "Bx": 8.949499297060736e-05, + "By": 0.0012339029977727672, + "Bz": 0.0012339029977727674, + "Ex": 5538900.045484711, + "Ey": 216461.93007927085, + "Ez": 216461.9300792709 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index e605e1a87f8..91070190772 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -1916,18 +1916,18 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Modules/resampling/analysis_leveling_thinning.py -[embedded_boundary_sphere] +[embedded_boundary_cube] buildDir = . -inputFile = Examples/Modules/embedded_boundary_sphere/inputs_3d +inputFile = Examples/Modules/embedded_boundary_cube/inputs_3d runtime_params = dim = 3 addToCompileString = USE_EB=TRUE restartTest = 0 useMPI = 1 -numprocs = 1 +numprocs = 4 useOMP = 1 numthreads = 1 compileTest = 0 doVis = 0 compareParticles = 0s -analysisRoutine = Examples/Modules/embedded_boundary_sphere/plots/analysis_fields.py \ No newline at end of file +analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py \ No newline at end of file diff --git a/Source/EmbeddedBoundary/CMakeLists.txt b/Source/EmbeddedBoundary/CMakeLists.txt index 7e11e2930f9..65a11610d18 100644 --- a/Source/EmbeddedBoundary/CMakeLists.txt +++ b/Source/EmbeddedBoundary/CMakeLists.txt @@ -1,5 +1,4 @@ target_sources(WarpX PRIVATE WarpXInitEB.cpp - WarpXInitEMField.cpp ) diff --git a/Source/EmbeddedBoundary/WarpXInitEMField.cpp b/Source/EmbeddedBoundary/WarpXInitEMField.cpp deleted file mode 100644 index 1fb53fce1e2..00000000000 --- a/Source/EmbeddedBoundary/WarpXInitEMField.cpp +++ /dev/null @@ -1,72 +0,0 @@ - -#include "WarpX.H" -#include - -amrex::RealArray -WarpX::AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere) { - amrex::Real k=2.7437/r_sphere; - amrex::Real r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); - amrex::Real theta = atan2(sqrt(pow(x, 2) + pow(y, 2)), z); - amrex::Real phi = atan2(y, x); - amrex::Real H_r = 0; - amrex::Real H_theta = 0; - amrex::Real mu_r = 1; //PhysConst::mu0; - amrex::Real H_phi = k / (r_sphere * mu_r) * std::sph_bessel(1, k*r) * sin(theta) * cos(k * PhysConst::c * t); - - amrex::Real H_x = H_r * sin(theta) * cos(phi) + H_theta * cos(theta) * cos(phi) - H_phi * sin(phi); - amrex::Real H_y = H_r * cos(theta) * cos(phi) + H_theta * cos(theta) * sin(phi) + H_phi * cos(phi); - amrex::Real H_z = H_r * cos(theta) - H_theta * sin(theta); - - return {PhysConst::mu0*H_x, PhysConst::mu0*H_y, PhysConst::mu0*H_z}; -} - -void -WarpX::InitEMFieldSphere (){ -#ifdef AMREX_USE_EB - amrex::Real r_sphere; - amrex::ParmParse pp_eb2("eb2"); - amrex::ParmParse pp_geometry("geometry"); - pp_eb2.get("sphere_radius", r_sphere); - amrex::RealArray prob_lo; - pp_geometry.get("prob_lo", prob_lo); - auto &Bfield = Bfield_fp[maxLevel()]; - - auto &face_areas_0 = m_face_areas[maxLevel()]; - for (amrex::MFIter mfi(*Bfield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - amrex::Array4 const &Bx = Bfield[0]->array(mfi); - amrex::Array4 const &By = Bfield[1]->array(mfi); - amrex::Array4 const &Bz = Bfield[2]->array(mfi); - amrex::Array4 const &Sx = face_areas_0[0]->array(mfi); - amrex::Array4 const &Sy = face_areas_0[1]->array(mfi); - amrex::Array4 const &Sz = face_areas_0[2]->array(mfi); - amrex::Box const &box = mfi.validbox(); - - amrex::LoopOnCpu(box, [=] (int i, int j, int k) { - amrex::Real x, y, z; - if (Sx(i, j, k) > 0) { - x = i * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - - Bx(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[0]; - } - - if (Sy(i, j, k) > 0) { - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = j * CellSize(maxLevel())[1] + prob_lo[1]; - z = (k + 0.5) * CellSize(maxLevel())[2] + prob_lo[2]; - - By(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[1]; - } - - if (Sz(i, j, k) > 0){ - x = (i + 0.5) * CellSize(maxLevel())[0] + prob_lo[0]; - y = (j + 0.5) * CellSize(maxLevel())[1] + prob_lo[1]; - z = k * CellSize(maxLevel())[2] + prob_lo[2]; - - Bz(i, j, k) = AnalyticSolSphere(x, y, z, 0., r_sphere)[2]; - } - }); - } -#endif -} diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index c680876bc03..d37411f289f 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -45,18 +45,6 @@ WarpX::Evolve (int numsteps) Real walltime, walltime_start = amrex::second(); -#ifdef AMREX_USE_EB - //If necessary initialize the em fields - std::string geom_type; - amrex::ParmParse pp_eb2("eb2"); - bool init_em_field; - pp_eb2.get("geom_type", geom_type); - pp_eb2.get("init_em_field", init_em_field); - if(init_em_field and geom_type=="sphere") { - InitEMFieldSphere(); - } -#endif - for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) { Real walltime_beg_step = amrex::second(); diff --git a/Source/WarpX.H b/Source/WarpX.H index 252d64b3f2a..0754876e5af 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1035,9 +1035,6 @@ private: void ScaleEdges (); void ScaleAreas (); - static amrex::RealArray AnalyticSolSphere (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t, amrex::Real r_sphere); - void InitEMFieldSphere (); - private: // void EvolvePSATD (int numsteps); void PushPSATD (amrex::Real dt); From 230cc571396c6dcaa9dcd9f484a86c5ce676202e Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 22 Apr 2021 19:02:02 +0200 Subject: [PATCH 23/35] fixed edges lengths computation and added comments --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 35 ++++++++++++++----------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index d43bac4df17..500786c5247 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -22,6 +22,10 @@ WarpX::InitEB () #endif } +/** + * \brief Compute the length of the mesh edges. Here the length is a value in [0, 1]. + * An edge of length 0 is fully covered. + */ void WarpX::ComputeEdgeLengths () { #ifdef AMREX_USE_EB @@ -34,26 +38,19 @@ WarpX::ComputeEdgeLengths () { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){ amrex::Box const &box = mfi.validbox(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){ - auto const &edge = edge_centroid[idim]->const_array(mfi); + auto const &edge_cent = edge_centroid[idim]->const_array(mfi); auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_cent).ixType()), [=](int i, int j, int k) { - if (edge(i, j, k) == amrex::Real(-1.0)) { + if (edge_cent(i, j, k) == amrex::Real(-1.0)) { // This edge is all covered edge_lengths_dim(i, j, k) = 0.; - } else if (edge(i, j, k) == amrex::Real(1.0)) { + } else if (edge_cent(i, j, k) == amrex::Real(1.0)) { // This edge is all open edge_lengths_dim(i, j, k) = 1.; } else { // This edge is cut. - amrex::Real edge_cent = edge(i, j, k); // edge centroid: (-0.5,0.5) - if (edge_cent < amrex::Real(0.0)) { - // The right side is covered. - edge_lengths_dim(i, j, k) = amrex::Real(2.0) * edge_cent + amrex::Real(0.5) + 0.5; // (0, 1) - } else { - // The left side is covered - edge_lengths_dim(i, j, k) = 0.5 - amrex::Real(2.0) * edge_cent - amrex::Real(0.5); // (0, 1) - } + edge_lengths_dim(i, j, k) = 1 - abs(amrex::Real(2.0)*edge_cent(i, j, k)); } }); } @@ -61,7 +58,10 @@ WarpX::ComputeEdgeLengths () { #endif } - +/** + * \brief Compute the ara of the mesh faces. Here the area is a value in [0, 1]. + * An edge of area 0 is fully covered. + */ void WarpX::ComputeFaceAreas () { #ifdef AMREX_USE_EB @@ -85,10 +85,13 @@ WarpX::ComputeFaceAreas () { #endif } +/** + * \brief Scale the edges lengths by the mesh width to obtain the real lengths. + */ void WarpX::ScaleEdges () { #ifdef AMREX_USE_EB - auto const& cell_size = CellSize(maxLevel()); + auto const &cell_size = CellSize(maxLevel()); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); auto const &edge_centroid = eb_fact.getEdgeCent(); @@ -107,7 +110,9 @@ WarpX::ScaleEdges () { #endif } - +/** + * \brief Scale the edges areas by the mesh width to obtain the real areas. + */ void WarpX::ScaleAreas() { #ifdef AMREX_USE_EB From 002a1f4afb510480c593a88ef6caea3f4164e115 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Thu, 22 Apr 2021 23:46:22 +0200 Subject: [PATCH 24/35] Fixed the case of all_regular geometries --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 88 +++++++++++++++++-------- 1 file changed, 59 insertions(+), 29 deletions(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 500786c5247..6e0df5279d4 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -15,7 +15,8 @@ WarpX::InitEB () amrex::ParmParse pp_eb2("eb2"); if (!pp_eb2.contains("geom_type")) { - pp_eb2.add("geom_type", "all_regular"); // use all_regular by default + std::string geom_type = "all_regular"; + pp_eb2.add("geom_type", geom_type); // use all_regular by default } amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel()); @@ -37,22 +38,37 @@ WarpX::ComputeEdgeLengths () { auto const &edge_centroid = eb_fact.getEdgeCent(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi){ amrex::Box const &box = mfi.validbox(); + amrex::FabType fab_type = flags[mfi].getType(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim){ - auto const &edge_cent = edge_centroid[idim]->const_array(mfi); auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_cent).ixType()), - [=](int i, int j, int k) { - if (edge_cent(i, j, k) == amrex::Real(-1.0)) { - // This edge is all covered - edge_lengths_dim(i, j, k) = 0.; - } else if (edge_cent(i, j, k) == amrex::Real(1.0)) { - // This edge is all open + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()), + [=](int i, int j, int k) { edge_lengths_dim(i, j, k) = 1.; - } else { - // This edge is cut. - edge_lengths_dim(i, j, k) = 1 - abs(amrex::Real(2.0)*edge_cent(i, j, k)); - } - }); + }); + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()), + [=](int i, int j, int k) { + edge_lengths_dim(i, j, k) = 1.; + }); + } else { + auto const &edge_cent = edge_centroid[idim]->const_array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_cent).ixType()), + [=](int i, int j, int k) { + if (edge_cent(i, j, k) == amrex::Real(-1.0)) { + // This edge is all covered + edge_lengths_dim(i, j, k) = 0.; + } else if (edge_cent(i, j, k) == amrex::Real(1.0)) { + // This edge is all open + edge_lengths_dim(i, j, k) = 1.; + } else { + // This edge is cut. + edge_lengths_dim(i, j, k) = 1 - abs(amrex::Real(2.0)*edge_cent(i, j, k)); + } + }); + } } } #endif @@ -73,13 +89,28 @@ WarpX::ComputeFaceAreas () { for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); + amrex::FabType fab_type = flags[mfi].getType(box); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &face = area_frac[idim]->const_array(mfi); - auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), - [=](int i, int j, int k) { - face_areas_dim(i, j, k) = face(i, j, k); - }); + auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + if (fab_type == amrex::FabType::regular) { + // every cell in box is all regular + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = amrex::Real(1.); + }); + } else if (fab_type == amrex::FabType::covered) { + // every cell in box is all covered + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = amrex::Real(0.); + }); + } else { + auto const &face = area_frac[idim]->const_array(mfi); + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + [=](int i, int j, int k) { + face_areas_dim(i, j, k) = face(i, j, k); + }); + } } } #endif @@ -91,20 +122,20 @@ WarpX::ComputeFaceAreas () { void WarpX::ScaleEdges () { #ifdef AMREX_USE_EB + BL_PROFILE("ScaleEdges"); + auto const &cell_size = CellSize(maxLevel()); auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); - auto const &edge_centroid = eb_fact.getEdgeCent(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - auto const &edge = edge_centroid[idim]->const_array(mfi); auto const &edge_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge).ixType()), + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()), [=](int i, int j, int k) { - edge_lengths_dim(i, j, k) *= cell_size[idim]; - }); + edge_lengths_dim(i, j, k) *= cell_size[idim]; + }); } } #endif @@ -116,12 +147,13 @@ WarpX::ScaleEdges () { void WarpX::ScaleAreas() { #ifdef AMREX_USE_EB + BL_PROFILE("ScaleAreas"); + auto const& cell_size = CellSize(maxLevel()); amrex::Real full_area; auto const eb_fact = fieldEBFactory(maxLevel()); auto const &flags = eb_fact.getMultiEBCellFlagFab(); - auto const &area_frac = eb_fact.getAreaFrac(); for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { amrex::Box const &box = mfi.validbox(); @@ -133,13 +165,11 @@ WarpX::ScaleAreas() { } else { full_area = cell_size[0]*cell_size[1]; } - auto const &face = area_frac[idim]->const_array(mfi); auto const &face_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); - amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face).ixType()), + amrex::LoopOnCpu(amrex::convert(box, amrex::Box(face_areas_dim).ixType()), [=](int i, int j, int k) { face_areas_dim(i, j, k) *= full_area; }); - } } #endif From c46c6bb470bc7c715eebe667fcf14b7fde90cd3e Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 08:28:11 +0200 Subject: [PATCH 25/35] fixing a bug that was breaking some tests --- Source/WarpX.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e060d42071c..460fe9c6191 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -229,10 +229,8 @@ WarpX::WarpX () Efield_avg_fp.resize(nlevs_max); Bfield_avg_fp.resize(nlevs_max); -#ifdef AMREX_USE_EB m_edge_lengths.resize(nlevs_max); m_face_areas.resize(nlevs_max); -#endif current_store.resize(nlevs_max); From 0cc72dd9fafaa5f366a2e3ce1e55cb656344b39f Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 08:29:45 +0200 Subject: [PATCH 26/35] adding test folder --- .../embedded_boundary_cube/analysis_fields.py | 110 ++++++++++++++++++ .../Modules/embedded_boundary_cube/inputs_3d | 37 ++++++ 2 files changed, 147 insertions(+) create mode 100644 Examples/Modules/embedded_boundary_cube/analysis_fields.py create mode 100644 Examples/Modules/embedded_boundary_cube/inputs_3d diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py new file mode 100644 index 00000000000..37a0aa5ae6f --- /dev/null +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -0,0 +1,110 @@ +#! /usr/bin/env python + +import yt +import os, sys +from scipy.constants import mu_0, pi, c +import numpy as np +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# This is a script that analyses the simulation results from +# the script `inputs_3d`. This simulates a TMmnp mode in a PEC cubic resonator. +# The magnetic field in the simulation is given (in theory) by: +# $$ B_x = \frac{-2\mu}{h^2}\, k_x k_z \sin(k_x x)\cos(k_y y)\cos(k_z z)\cos( \omega_p t)$$ +# $$ B_y = \frac{-2\mu}{h^2}\, k_y k_z \cos(k_x x)\sin(k_y y)\cos(k_z z)\cos( \omega_p t)$$ +# $$ B_z = \cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$ +# with +# $$ h^2 = k_x^2 + k_y^2 + k_z^2$$ +# $$ k_x = \frac{m\pi}{L}$$ +# $$ k_y = \frac{n\pi}{L}$$ +# $$ k_z = \frac{p\pi}{L}$$ + +hi = [0.8, 0.8, 0.8] +lo = [-0.8, -0.8, -0.8] +ncells = [48, 48, 48] +dx = (hi[0] - lo[0])/ncells[0] +dy = (hi[1] - lo[1])/ncells[1] +dz = (hi[2] - lo[2])/ncells[2] +m = 0 +n = 1 +p = 1 +Lx = 1 +Ly = 1 +Lz = 1 +h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2 +t = 1.3342563807926085e-08 + +# Compute the analytic solution +Bx_th = np.zeros(ncells) +By_th = np.zeros(ncells) +Bz_th = np.zeros(ncells) +for i in range(ncells[0]): + for j in range(ncells[1]): + for k in range(ncells[2]): + x = i*dx + lo[0] + y = j*dy + lo[1] + z = k*dz + lo[2] + + Bx_th[i, j, k] = -2/h_2*mu_0*(m * pi/Lx)*(p * pi/Lz) * (np.sin(m * pi/Lx * (x - Lx/2)) * + np.sin(n * pi/Ly * (y - Ly/2)) * + np.cos(p * pi/Lz * (z - Lz/2)) * + (-Lx/2 < x < Lx/2) * + (-Ly/2 < y < Ly/2) * + (-Lz/2 < z < Lz/2) * + np.cos(np.sqrt(2) * + np.pi / Lx * c * t)) + + By_th[i, j, k] = -2/h_2*mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * + np.sin(n * pi/Ly * (y - Ly/2)) * + np.cos(p * pi/Lz * (z - Lz/2)) * + (-Lx/2 < x < Lx/2) * + (-Ly/2 < y < Ly/2) * + (-Lz/2 < z < Lz/2) * + np.cos(np.sqrt(2) * + np.pi / Lx * c * t)) + + Bz_th[i, j, k] = mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * + np.cos(n * pi/Ly * (y - Ly/2)) * + np.sin(p * pi/Lz * (z - Lz/2)) * + (-Lx/2 < x < Lx/2) * + (-Ly/2 < y < Ly/2) * + (-Lz/2 < z < Lz/2) * + np.cos(np.sqrt(2) * + np.pi / Lx * c * t)) + +# Compute the analytic solution +path, dirs, files = next(os.walk("diags/")) +tsteps = np.zeros_like(dirs, dtype=int) +for i, d in enumerate(dirs): + # Strip off diags + d = d[4:] + # Strip off leading zeros + while d[0] == 0: + d = d[1:] + tsteps[i] = int(d) + +# Open the right plot file +filename = 'diags/diag00208/' + +ds = yt.load(filename) +data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) + +tol_err = 1e-5 + +# Compute l^2 error on Bx +Bx_sim = data['Bx'].to_ndarray() +err_x = np.sqrt(np.sum(np.square(Bx_sim - Bx_th))*dx*dy*dz) +assert(err_x < tol_err) +# Compute l^2 error on By +By_sim = data['By'].to_ndarray() +err_y = np.sqrt(np.sum(np.square(By_sim - By_th))*dx*dy*dz) +assert(err_y < tol_err) + +# Compute l^2 error on Bz +Bz_sim = data['Bz'].to_ndarray() +err_z = np.sqrt(np.sum(np.square(Bz_sim - Bz_th))*dx*dy*dz) +assert(err_z < tol_err) + +test_name = os.path.split(os.getcwd())[1] + +checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Examples/Modules/embedded_boundary_cube/inputs_3d b/Examples/Modules/embedded_boundary_cube/inputs_3d new file mode 100644 index 00000000000..dacf2efa200 --- /dev/null +++ b/Examples/Modules/embedded_boundary_cube/inputs_3d @@ -0,0 +1,37 @@ +stop_time = 1.3342563807926085e-08 +amr.n_cell = 48 48 48 +amr.max_grid_size = 128 +amr.max_level = 0 + +geometry.coord_sys = 0 +geometry.is_periodic = 0 0 0 +geometry.prob_lo = -0.8 -0.8 -0.8 +geometry.prob_hi = 0.8 0.8 0.8 +warpx.do_pml = 0 +warpx.const_dt = 1e-6 +warpx.cfl = 1 + +eb2.geom_type = box +eb2.box_lo = -0.5 -0.5 -0.5 +eb2.box_hi = 0.5 0.5 0.5 +eb2.box_has_fluid_inside = true + +warpx.B_ext_grid_init_style = parse_B_ext_grid_function +my_constants.m = 0 +my_constants.n = 1 +my_constants.p = 1 +my_constants.pi = 3.141592653589793 +my_constants.Lx = 1 +my_constants.Ly = 1 +my_constants.Lz = 1 +my_constants.h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2 +my_constants.mu_0 = 1.25663706212e-06 + +warpx.By_external_grid_function(x,y,z) = -2/h_2 * (n * pi / Ly) * (p * pi / Lz) * cos(m * pi / Lx * (x - Lx / 2)) * sin(n * pi / Ly * (y - Ly / 2)) * cos(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-Lx/2)*(x-Ly/2)*(y-Lz/2)*(z-Lx/2)*(x-Ly/2)*(y-Lz/2)*(z-0.5)*(x<0.5)*(y>-0.5)*(y<0.5)*(z>-0.5)*(z<0.5) + +diagnostics.diags_names = diag +diag.intervals = 1 +diag.diag_type = Full +diag.fields_to_plot = Ex Ey Ez Bx By Bz From 2f595d8c123fdf07b669b236359606a1bfd5aa1a Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 10:24:18 +0200 Subject: [PATCH 27/35] fixed the default values for the EB cube test --- .../benchmarks_json/embedded_boundary_cube.json | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json index bf0f1e525fb..adbd1fbf782 100644 --- a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json @@ -1,10 +1,10 @@ { "lev=0": { - "Bx": 8.949499297060736e-05, - "By": 0.0012339029977727672, - "Bz": 0.0012339029977727674, - "Ex": 5538900.045484711, - "Ey": 216461.93007927085, - "Ez": 216461.9300792709 + "Bx": 0.0, + "By": 0.0066283741197868335, + "Bz": 0.006628374119786833, + "Ex": 5102618.47115243, + "Ey": 0.0, + "Ez": 0.0 } -} \ No newline at end of file +} From c0396be3d6217334d34c4d525bb162d1524bc755 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 11:44:55 +0200 Subject: [PATCH 28/35] simplified the analysis script --- .../Modules/embedded_boundary_cube/analysis_fields.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py index 37a0aa5ae6f..ed39d4a51ec 100644 --- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -72,16 +72,6 @@ np.cos(np.sqrt(2) * np.pi / Lx * c * t)) -# Compute the analytic solution -path, dirs, files = next(os.walk("diags/")) -tsteps = np.zeros_like(dirs, dtype=int) -for i, d in enumerate(dirs): - # Strip off diags - d = d[4:] - # Strip off leading zeros - while d[0] == 0: - d = d[1:] - tsteps[i] = int(d) # Open the right plot file filename = 'diags/diag00208/' From 529d0705b5709a30616f771a34670678938baf77 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 11:49:11 +0200 Subject: [PATCH 29/35] fixed cubic resonator default results --- Regression/Checksum/benchmarks_json/embedded_boundary_cube.json | 1 + 1 file changed, 1 insertion(+) diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json index adbd1fbf782..69fffed2c6a 100644 --- a/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json @@ -8,3 +8,4 @@ "Ez": 0.0 } } + From b8587dcc16452d2e67324b3326ccaad3d4b5fe7f Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 15:11:10 +0200 Subject: [PATCH 30/35] inputting the plot file name from command line --- Examples/Modules/embedded_boundary_cube/analysis_fields.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py index ed39d4a51ec..9b8ef5bbb50 100644 --- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -74,8 +74,7 @@ # Open the right plot file -filename = 'diags/diag00208/' - +filename = sys.argv[1] ds = yt.load(filename) data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) From 39e76afe9dbb8e082f90406d30fce1ec133bafc1 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Fri, 23 Apr 2021 16:22:11 +0200 Subject: [PATCH 31/35] fixing the diag name --- Examples/Modules/embedded_boundary_cube/inputs_3d | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Examples/Modules/embedded_boundary_cube/inputs_3d b/Examples/Modules/embedded_boundary_cube/inputs_3d index dacf2efa200..f61f51fdfa4 100644 --- a/Examples/Modules/embedded_boundary_cube/inputs_3d +++ b/Examples/Modules/embedded_boundary_cube/inputs_3d @@ -31,7 +31,7 @@ warpx.By_external_grid_function(x,y,z) = -2/h_2 * (n * pi / Ly) * (p * pi / Lz) warpx.Bx_external_grid_function(x,y,z) = -2/h_2 * (m * pi / Lx) * (p * pi / Lz) * sin(m * pi / Lx * (x - Lx / 2)) * cos(n * pi / Ly * (y - Ly / 2)) * cos(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-Lx/2)*(x-Ly/2)*(y-Lz/2)*(z-0.5)*(x<0.5)*(y>-0.5)*(y<0.5)*(z>-0.5)*(z<0.5) -diagnostics.diags_names = diag -diag.intervals = 1 -diag.diag_type = Full -diag.fields_to_plot = Ex Ey Ez Bx By Bz +diagnostics.diags_names = diag1 +diag1.intervals = 1 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez Bx By Bz From c6f86383d7ff1b5a3dc2d50bf4ecdeb120292ff0 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel <47607756+lgiacome@users.noreply.github.com> Date: Tue, 27 Apr 2021 09:25:31 +0200 Subject: [PATCH 32/35] Fixed a bug in edges initialization Co-authored-by: Remi Lehe --- Source/EmbeddedBoundary/WarpXInitEB.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index 6e0df5279d4..8214eec9b16 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -51,7 +51,7 @@ WarpX::ComputeEdgeLengths () { // every cell in box is all covered amrex::LoopOnCpu(amrex::convert(box, amrex::Box(edge_lengths_dim).ixType()), [=](int i, int j, int k) { - edge_lengths_dim(i, j, k) = 1.; + edge_lengths_dim(i, j, k) = 0.; }); } else { auto const &edge_cent = edge_centroid[idim]->const_array(mfi); From e56107e80338655c520013de8ca3001a8320ff0d Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel <47607756+lgiacome@users.noreply.github.com> Date: Tue, 27 Apr 2021 09:27:46 +0200 Subject: [PATCH 33/35] Adding comments to the staircased yee solver (thanks Remi) Co-authored-by: Remi Lehe --- Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp | 3 +++ Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp | 3 +++ 2 files changed, 6 insertions(+) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 670effb2fc4..c24b869e047 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -110,6 +110,7 @@ void FiniteDifferenceSolver::EvolveBCartesian ( [=] 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 (Sx(i, j, k) <= 0) return; #endif Bx(i, j, k) += dt * T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k) @@ -118,6 +119,7 @@ void FiniteDifferenceSolver::EvolveBCartesian ( [=] 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 (Sy(i, j, k) <= 0) return; #endif By(i, j, k) += dt * T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k) @@ -126,6 +128,7 @@ void FiniteDifferenceSolver::EvolveBCartesian ( [=] 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 (Sz(i, j, k) <= 0) return; #endif Bz(i, j, k) += dt * T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k) diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index 77d431c9f81..3cb433e9776 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -120,6 +120,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] 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 (lx(i, j, k) <= 0) return; #endif Ex(i, j, k) += c2 * dt * ( @@ -130,6 +131,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] 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 (ly(i,j,k) <= 0) return; #endif Ey(i, j, k) += c2 * dt * ( @@ -140,6 +142,7 @@ void FiniteDifferenceSolver::EvolveECartesian ( [=] 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; #endif Ez(i, j, k) += c2 * dt * ( From 2f85151a908b460ba09798b22b517b983146e9e9 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Tue, 27 Apr 2021 10:39:03 +0200 Subject: [PATCH 34/35] fixed the cube resonator test --- .../embedded_boundary_cube/analysis_fields.py | 56 ++++++++----------- .../Modules/embedded_boundary_cube/inputs_3d | 2 +- 2 files changed, 23 insertions(+), 35 deletions(-) diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py index 9b8ef5bbb50..83d1685d286 100644 --- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -6,6 +6,7 @@ import numpy as np sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI +import matplotlib.pyplot as plt # This is a script that analyses the simulation results from # the script `inputs_3d`. This simulates a TMmnp mode in a PEC cubic resonator. @@ -42,57 +43,44 @@ for j in range(ncells[1]): for k in range(ncells[2]): x = i*dx + lo[0] - y = j*dy + lo[1] + y = (j+0.5)*dy + lo[1] z = k*dz + lo[2] - Bx_th[i, j, k] = -2/h_2*mu_0*(m * pi/Lx)*(p * pi/Lz) * (np.sin(m * pi/Lx * (x - Lx/2)) * - np.sin(n * pi/Ly * (y - Ly/2)) * - np.cos(p * pi/Lz * (z - Lz/2)) * - (-Lx/2 < x < Lx/2) * - (-Ly/2 < y < Ly/2) * - (-Lz/2 < z < Lz/2) * - np.cos(np.sqrt(2) * - np.pi / Lx * c * t)) - By_th[i, j, k] = -2/h_2*mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * np.sin(n * pi/Ly * (y - Ly/2)) * np.cos(p * pi/Lz * (z - Lz/2)) * - (-Lx/2 < x < Lx/2) * - (-Ly/2 < y < Ly/2) * - (-Lz/2 < z < Lz/2) * + (-Lx/2 <= x < Lx/2) * + (-Ly/2 <= y < Ly/2) * + (-Lz/2 <= z < Lz/2) * np.cos(np.sqrt(2) * np.pi / Lx * c * t)) - Bz_th[i, j, k] = mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * - np.cos(n * pi/Ly * (y - Ly/2)) * - np.sin(p * pi/Lz * (z - Lz/2)) * - (-Lx/2 < x < Lx/2) * - (-Ly/2 < y < Ly/2) * - (-Lz/2 < z < Lz/2) * - np.cos(np.sqrt(2) * - np.pi / Lx * c * t)) - + x = i*dx + lo[0] + y = j*dy + lo[1] + z = (k+0.5)*dz + lo[2] + Bz_th[i, j, k] = mu_0*(np.cos(m * pi/Lx * (x - Lx/2)) * + np.cos(n * pi/Ly * (y - Ly/2)) * + np.sin(p * pi/Lz * (z - Lz/2)) * + (-Lx/2 <= x < Lx/2) * + (-Ly/2 <= y < Ly/2) * + (-Lz/2 <= z < Lz/2) * + np.cos(np.sqrt(2) * np.pi / Lx * c * t)) # Open the right plot file filename = sys.argv[1] ds = yt.load(filename) data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) -tol_err = 1e-5 +rel_tol_err = 1e-1 -# Compute l^2 error on Bx -Bx_sim = data['Bx'].to_ndarray() -err_x = np.sqrt(np.sum(np.square(Bx_sim - Bx_th))*dx*dy*dz) -assert(err_x < tol_err) -# Compute l^2 error on By +# Compute relative l^2 error on By By_sim = data['By'].to_ndarray() -err_y = np.sqrt(np.sum(np.square(By_sim - By_th))*dx*dy*dz) -assert(err_y < tol_err) - -# Compute l^2 error on Bz +rel_err_y = np.sqrt( np.sum(np.square(By_sim - By_th)) / np.sum(np.square(By_th))) +assert(rel_err_y < rel_tol_err) +# Compute relative l^2 error on Bz Bz_sim = data['Bz'].to_ndarray() -err_z = np.sqrt(np.sum(np.square(Bz_sim - Bz_th))*dx*dy*dz) -assert(err_z < tol_err) +rel_err_z = np.sqrt( np.sum(np.square(Bz_sim - Bz_th)) / np.sum(np.square(Bz_th))) +assert(rel_err_z < rel_tol_err) test_name = os.path.split(os.getcwd())[1] diff --git a/Examples/Modules/embedded_boundary_cube/inputs_3d b/Examples/Modules/embedded_boundary_cube/inputs_3d index f61f51fdfa4..4817075405f 100644 --- a/Examples/Modules/embedded_boundary_cube/inputs_3d +++ b/Examples/Modules/embedded_boundary_cube/inputs_3d @@ -32,6 +32,6 @@ warpx.Bx_external_grid_function(x,y,z) = -2/h_2 * (m * pi / Lx) * (p * pi / Lz) warpx.Bz_external_grid_function(x,y,z) = cos(m * pi / Lx * (x - Lx / 2)) * cos(n * pi / Ly * (y - Ly / 2)) * sin(p * pi / Lz * (z - Lz / 2))*mu_0*(x>-0.5)*(x<0.5)*(y>-0.5)*(y<0.5)*(z>-0.5)*(z<0.5) diagnostics.diags_names = diag1 -diag1.intervals = 1 +diag1.intervals = 1000 diag1.diag_type = Full diag1.fields_to_plot = Ex Ey Ez Bx By Bz From 211bd5199081b9be3877e9cfd1e440fb4bc724df Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Tue, 27 Apr 2021 11:30:45 +0200 Subject: [PATCH 35/35] removed an unused import --- Examples/Modules/embedded_boundary_cube/analysis_fields.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py index 83d1685d286..58dad1bb7fa 100644 --- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -6,7 +6,6 @@ import numpy as np sys.path.insert(1, '../../../../warpx/Regression/Checksum/') import checksumAPI -import matplotlib.pyplot as plt # This is a script that analyses the simulation results from # the script `inputs_3d`. This simulates a TMmnp mode in a PEC cubic resonator.