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..58dad1bb7fa --- /dev/null +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -0,0 +1,86 @@ +#! /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+0.5)*dy + lo[1] + z = k*dz + lo[2] + + 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)) + + 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) + +rel_tol_err = 1e-1 + +# Compute relative l^2 error on By +By_sim = data['By'].to_ndarray() +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() +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] + +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..4817075405f --- /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 = diag1 +diag1.intervals = 1000 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez Bx By Bz 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..69fffed2c6a --- /dev/null +++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube.json @@ -0,0 +1,11 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.0066283741197868335, + "Bz": 0.006628374119786833, + "Ex": 5102618.47115243, + "Ey": 0.0, + "Ez": 0.0 + } +} + diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index ad964f7f028..cfc937220e1 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2002,3 +2002,19 @@ compileTest = 0 doVis = 0 compareParticles = 0 analysisRoutine = Examples/Modules/resampling/analysis_leveling_thinning.py + +[embedded_boundary_cube] +buildDir = . +inputFile = Examples/Modules/embedded_boundary_cube/inputs_3d +runtime_params = +dim = 3 +addToCompileString = USE_EB=TRUE +restartTest = 0 +useMPI = 1 +numprocs = 4 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0s +analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py \ No newline at end of file diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp index ec17d9efccc..8214eec9b16 100644 --- a/Source/EmbeddedBoundary/WarpXInitEB.cpp +++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp @@ -15,8 +15,162 @@ 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()); + +#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 + 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(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_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + 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 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) = 0.; + }); + } 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 +} + +/** + * \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 + 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); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + 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 +} + +/** + * \brief Scale the edges lengths by the mesh width to obtain the real lengths. + */ +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(); + + 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_lengths_dim = m_edge_lengths[maxLevel()][idim]->array(mfi); + 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]; + }); + } + } +#endif +} + +/** + * \brief Scale the edges areas by the mesh width to obtain the real areas. + */ +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(); + + for (amrex::MFIter mfi(flags); mfi.isValid(); ++mfi) { + amrex::Box const &box = mfi.validbox(); + 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_areas_dim = m_face_areas[maxLevel()][idim]->array(mfi); + 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 } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 9b79c0a8409..c445e48c52f 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -44,6 +44,7 @@ WarpX::Evolve (int numsteps) bool early_params_checked = false; // check typos in inputs after step 1 Real walltime, walltime_start = amrex::second(); + 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..c24b869e047 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 + 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 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,28 @@ void FiniteDifferenceSolver::EvolveBCartesian ( amrex::ParallelFor(tbx, tby, tbz, [=] 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) - 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 + // 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) - 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 + // 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) - 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..3cb433e9776 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 + 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 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,10 @@ void FiniteDifferenceSolver::EvolveECartesian ( amrex::ParallelFor(tex, tey, tez, [=] 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 * ( - 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 +130,10 @@ 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 * ( - 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 +141,10 @@ 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 * ( - 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..a991c539dc7 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], m_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], 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], 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], 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/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 87c70fba01b..0489dfc8b3d 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 b63bc6e5650..bc035f38358 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -908,6 +908,11 @@ private: amrex::Vector, 3 > > Bfield_fp; amrex::Vector, 3 > > Efield_avg_fp; amrex::Vector, 3 > > Bfield_avg_fp; + + //EB grid info + amrex::Vector, 3 > > m_edge_lengths; + amrex::Vector, 3 > > m_face_areas; + // store fine patch amrex::Vector, 3 > > current_store; @@ -1069,6 +1074,11 @@ private: void InitEB (); + void ComputeEdgeLengths (); + void ComputeFaceAreas (); + void ScaleEdges (); + void ScaleAreas (); + private: // void EvolvePSATD (int numsteps); void PushPSATD (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index de5a8acc4d0..460fe9c6191 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -179,7 +179,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. @@ -230,6 +229,9 @@ WarpX::WarpX () Efield_avg_fp.resize(nlevs_max); Bfield_avg_fp.resize(nlevs_max); + m_edge_lengths.resize(nlevs_max); + m_face_areas.resize(nlevs_max); + current_store.resize(nlevs_max); if (do_current_centering) @@ -1319,6 +1321,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()) + { + 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 + 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)