From 6dfb11f5875a8d172ba11b4b4bb612f45c6d1f75 Mon Sep 17 00:00:00 2001 From: ablelly Date: Fri, 7 Jun 2019 20:42:12 +0200 Subject: Added a flag for particles in PMLs --- Source/WarpX.cpp | 1 + 1 file changed, 1 insertion(+) (limited to 'Source/WarpX.cpp') diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c2cf97f30..829f151b0 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -379,6 +379,7 @@ WarpX::ReadParameters () pp.query("do_pml", do_pml); pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); + pp.query("pml_has_particles", pml_has_particles); pp.query("dump_openpmd", dump_openpmd); pp.query("dump_plotfiles", dump_plotfiles); -- cgit v1.2.3 From 5daf8d3a3c4987a332197bcaebc931660bbb076a Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 13 Jun 2019 01:38:22 +0200 Subject: Added a flag to do j-damping in pmls --- Source/Evolve/WarpXEvolveEM.cpp | 7 ++++++- Source/WarpX.H | 1 + Source/WarpX.cpp | 1 + 3 files changed, 8 insertions(+), 1 deletion(-) (limited to 'Source/WarpX.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 30768f31f..4c0ca3c70 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -314,7 +314,12 @@ WarpX::OneStep_nosub (Real cur_time) } } } - // DAMP J!!!!!! + + if (do_pml && do_pml_j_damping){ + // damp current in pmls + amrex::Print() << "===== DAMPING IN PMLs =====" << std::endl; + } + EvolveF(0.5*dt[0], DtType::FirstHalf); FillBoundaryF(); EvolveB(0.5*dt[0]); // We now have B^{n+1/2} diff --git a/Source/WarpX.H b/Source/WarpX.H index f08955ebf..bdaa170a6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -490,6 +490,7 @@ private: int pml_ncell = 10; int pml_delta = 10; int pml_has_particles = 0; + int do_pml_j_damping = 0; amrex::Vector > pml; amrex::Real moving_window_x = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 829f151b0..5f976b8f3 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -380,6 +380,7 @@ WarpX::ReadParameters () pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); pp.query("pml_has_particles", pml_has_particles); + pp.query("do_pml_j_damping", do_pml_j_damping); pp.query("dump_openpmd", dump_openpmd); pp.query("dump_plotfiles", dump_plotfiles); -- cgit v1.2.3 From b6834f615a8085f91ce07dfa1413d0d9ea4cc24a Mon Sep 17 00:00:00 2001 From: ablelly Date: Mon, 15 Jul 2019 23:19:01 +0200 Subject: Work on the merge --- Source/BoundaryConditions/PML.H | 8 +- Source/BoundaryConditions/PML.cpp | 105 ++- Source/BoundaryConditions/PML_Cplx.cpp | 958 --------------------------- Source/BoundaryConditions/WarpXEvolvePML.cpp | 6 +- Source/Evolve/WarpXEvolveEM.cpp | 22 +- Source/WarpX.H | 1 + Source/WarpX.cpp | 5 + 7 files changed, 56 insertions(+), 1049 deletions(-) delete mode 100644 Source/BoundaryConditions/PML_Cplx.cpp (limited to 'Source/WarpX.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index fd8fe12ec..62a378986 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -71,9 +71,9 @@ struct SigmaBox using SigmaVect = std::array; - SigmaVect sigma; // sigma/epsilon - SigmaVect sigma_cum; // cumsum(sigma)/(c*epsilon) - SigmaVect sigma_star; // sigma_star/mu + SigmaVect sigma; + SigmaVect sigma_cum; + SigmaVect sigma_star; SigmaVect sigma_star_cum; SigmaVect sigma_fac; SigmaVect sigma_cum_fac; @@ -81,7 +81,7 @@ struct SigmaBox SigmaVect sigma_star_cum_fac; std::string pml_type; // if pml normal to x : if pml Hi : pml_type={1,0,0}; if pml Lo: pml_type={-1,0,0} - std::array pml_type_array; + std::array pml_type_array; }; namespace amrex { diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index b01a95f10..1f1cc80c3 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -84,9 +84,6 @@ namespace SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) { - amrex::Print()<< "===== BUILD SigmaBox ====="<Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if ( ! geom->isPeriodic(idim) ) { domain0.grow(idim, -ncell); } } - // Box domain0 = amrex::grow(geom->Domain(), -ncell); - amrex::Print() << "[" << domain0.smallEnd()[0]<<", "<< domain0.smallEnd()[1]<<", "<< domain0.smallEnd()[2]<< ", "<setVal(0.0); pml_j_fp[2]->setVal(0.0); - // if (pml_has_particles){ - // pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) - // pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) - // pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) - // pml_j_fp[0]->setVal(0.0); - // pml_j_fp[1]->setVal(0.0); - // pml_j_fp[2]->setVal(0.0); - // amrex::Print() << "PML HAS PARTICLES - fine"<< std::endl; - // - // } - if (do_dive_cleaning) { pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf)); @@ -605,16 +571,6 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_j_cp[0]->setVal(0.0); pml_j_cp[1]->setVal(0.0); pml_j_cp[2]->setVal(0.0); - // if (pml_has_particles) - // { - // pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); - // pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); - // pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); - // pml_j_cp[0]->setVal(0.0); - // pml_j_cp[1]->setVal(0.0); - // pml_j_cp[2]->setVal(0.0); - // amrex::Print() << "PML HAS PARTICLES - coarse"<< std::endl; - // } // sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); @@ -679,15 +635,6 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, BoxArray ba(bl); ba.removeOverlap(false); - // BoxList bl_2 = BoxList(ba); - // - // amrex::Print() << "Printing PML boxes AFTER cleaning" << std::endl; - // amrex::Print() << "[" << std::endl; - // for (const Box& b: bl_2) { - // amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "< 0) // Copy from pml to the ghost cells of regular data +// { +// MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); +// MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); +// if (ncp == 3) { +// MultiFab::Add(totpmlmf,pml,2,0,1,0); +// } +// +// MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); +// tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); +// +// #ifdef _OPENMP +// #pragma omp parallel +// #endif +// for (MFIter mfi(reg); mfi.isValid(); ++mfi) +// { +// const FArrayBox& src = tmpregmf[mfi]; +// FArrayBox& dst = reg[mfi]; +// const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); +// for (const Box& bx : bl) +// { +// dst.copy(src, bx, 0, bx, 0, 1); +// } +// } +// } +// +// // Copy from regular data to PML's first component +// // Zero out the second (and third) component +// MultiFab::Copy(tmpregmf,reg,0,0,1,0); +// tmpregmf.setVal(0.0, 1, ncp-1, 0); +// pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period); +// } + void PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) { @@ -924,16 +914,8 @@ PML::CopyRegInPMLs (MultiFab& pml, MultiFab& reg, const Geometry& geom) // #endif const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); - // const int ncp = pml.nComp(); const auto& period = geom.periodicity(); - // MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), 1, ngr); - // tmpregmf.setVal(0.0, 0, 1, ngr); - // MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, ngp); - // totpmlmf.setVal(0.0, 0, 1, ngp); - // realise sum of splitted fields inside pml - - pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period); } @@ -948,7 +930,6 @@ PML::CopyPMLsInReg (MultiFab& pml, MultiFab& reg, const Geometry& geom) const IntVect& ngp = pml.nGrowVect(); const auto& period = geom.periodicity(); - // reg.ParallelCopy(pml, 0, 0, 1, IntVect(0), ngr, period); reg.ParallelCopy(pml, 0, 0, 1, ngp, ngr, period); } diff --git a/Source/BoundaryConditions/PML_Cplx.cpp b/Source/BoundaryConditions/PML_Cplx.cpp deleted file mode 100644 index 0a47c5d84..000000000 --- a/Source/BoundaryConditions/PML_Cplx.cpp +++ /dev/null @@ -1,958 +0,0 @@ - -#include -#include -#include - -#include -#include - -#include - -#ifdef _OPENMP -#include -#endif - -using namespace amrex; - -namespace -{ - static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star, - const Box& overlap, const Box& grid, Real fac) - { - int glo = grid.smallEnd(idim); - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; - for (int i = olo; i <= ohi+1; ++i) - { - Real offset = static_cast(glo-i); - sigma[i-slo] = fac*(offset*offset); - } - for (int i = olo; i <= ohi; ++i) - { - Real offset = static_cast(glo-i) - 0.5; - sigma_star[i-sslo] = fac*(offset*offset); - } - } - - static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star, - const Box& overlap, const Box& grid, Real fac) - { - int ghi = grid.bigEnd(idim); - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; - for (int i = olo; i <= ohi+1; ++i) - { - Real offset = static_cast(i-ghi-1); - sigma[i-slo] = fac*(offset*offset); - } - for (int i = olo; i <= ohi; ++i) - { - Real offset = static_cast(i-ghi) - 0.5; - sigma_star[i-sslo] = fac*(offset*offset); - } - } - - static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap) - { - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; - std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0); - std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0); - } -} - -SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) -{ - BL_ASSERT(box.cellCentered()); - - const IntVect& sz = box.size(); - const int* lo = box.loVect(); - const int* hi = box.hiVect(); - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { - sigma [idim].resize(sz[idim]+1); - sigma_star [idim].resize(sz[idim] ); - sigma_fac [idim].resize(sz[idim]+1); - sigma_star_fac[idim].resize(sz[idim] ); - - sigma [idim].m_lo = lo[idim]; - sigma [idim].m_hi = hi[idim]+1; - sigma_star [idim].m_lo = lo[idim]; - sigma_star [idim].m_hi = hi[idim]; - sigma_fac [idim].m_lo = lo[idim]; - sigma_fac [idim].m_hi = hi[idim]+1; - sigma_star_fac[idim].m_lo = lo[idim]; - sigma_star_fac[idim].m_hi = hi[idim]; - } - - Array fac; - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - fac[idim] = 4.0*PhysConst::c/(dx[idim]*static_cast(delta*delta)); - } - - const std::vector >& isects = grids.intersections(box, false, ncell); - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { - int jdim = (idim+1) % AMREX_SPACEDIM; -#if (AMREX_SPACEDIM == 3) - int kdim = (idim+2) % AMREX_SPACEDIM; -#endif - - Vector direct_faces, side_faces, direct_side_edges, side_side_edges, corners; - for (const auto& kv : isects) - { - const Box& grid_box = grids[kv.first]; - - if (amrex::grow(grid_box, idim, ncell).intersects(box)) - { - direct_faces.push_back(kv.first); - } - else if (amrex::grow(grid_box, jdim, ncell).intersects(box)) - { - side_faces.push_back(kv.first); - } -#if (AMREX_SPACEDIM == 3) - else if (amrex::grow(grid_box, kdim, ncell).intersects(box)) - { - side_faces.push_back(kv.first); - } - else if (amrex::grow(amrex::grow(grid_box,idim,ncell), - jdim,ncell).intersects(box)) - { - direct_side_edges.push_back(kv.first); - } - else if (amrex::grow(amrex::grow(grid_box,idim,ncell), - kdim,ncell).intersects(box)) - { - direct_side_edges.push_back(kv.first); - } - else if (amrex::grow(amrex::grow(grid_box,jdim,ncell), - kdim,ncell).intersects(box)) - { - side_side_edges.push_back(kv.first); - } -#endif - else - { - corners.push_back(kv.first); - } - } - - for (auto gid : corners) - { - const Box& grid_box = grids[gid]; - - Box lobox = amrex::adjCellLo(grid_box, idim, ncell); - lobox.grow(jdim,ncell); -#if (AMREX_SPACEDIM == 3) - lobox.grow(kdim,ncell); -#endif - Box looverlap = lobox & box; - if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); - } - - Box hibox = amrex::adjCellHi(grid_box, idim, ncell); - hibox.grow(jdim,ncell); -#if (AMREX_SPACEDIM == 3) - hibox.grow(kdim,ncell); -#endif - Box hioverlap = hibox & box; - if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); - } - - if (!looverlap.ok() && !hioverlap.ok()) { - amrex::Abort("SigmaBox::SigmaBox(): corners, how did this happen?\n"); - } - } - -#if (AMREX_SPACEDIM == 3) - for (auto gid : side_side_edges) - { - const Box& grid_box = grids[gid]; - const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; - if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); - } else { - amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n"); - } - } - - for (auto gid : direct_side_edges) - { - const Box& grid_box = grids[gid]; - - Box lobox = amrex::adjCellLo(grid_box, idim, ncell); - Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box; - if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); - } - - Box hibox = amrex::adjCellHi(grid_box, idim, ncell); - Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box; - if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); - } - - if (!looverlap.ok() && !hioverlap.ok()) { - amrex::Abort("SigmaBox::SigmaBox(): direct_side_edges, how did this happen?\n"); - } - } -#endif - - for (auto gid : side_faces) - { - const Box& grid_box = grids[gid]; -#if (AMREX_SPACEDIM == 2) - const Box& overlap = amrex::grow(grid_box,jdim,ncell) & box; -#else - const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; -#endif - if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); - } else { - amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); - } - } - - for (auto gid : direct_faces) - { - const Box& grid_box = grids[gid]; - - const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell); - Box looverlap = lobox & box; - if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); - } - - const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); - Box hioverlap = hibox & box; - if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); - } - - if (!looverlap.ok() && !hioverlap.ok()) { - amrex::Abort("SigmaBox::SigmaBox(): direct faces, how did this happen?\n"); - } - } - - if (direct_faces.size() > 1) { - amrex::Abort("SigmaBox::SigmaBox(): direct_faces.size() > 1, Box gaps not wide enough?\n"); - } - } -} - -void -SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) -{ - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { - for (int i = 0, N = sigma_star[idim].size(); i < N; ++i) - { - if (sigma_star[idim][i] == 0.0) - { - sigma_star_fac[idim][i] = 1.0; - } - else - { - sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt); - } - } - } -} - -void -SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) -{ - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { - for (int i = 0, N = sigma[idim].size(); i < N; ++i) - { - if (sigma[idim][i] == 0.0) - { - sigma_fac[idim][i] = 1.0; - } - else - { - sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt); - } - } - } -} - -MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm, - const BoxArray& grid_ba, const Real* dx, int ncell, int delta) - : FabArray(ba,dm,1,0,MFInfo(), - FabFactory(grid_ba,dx,ncell,delta)) -{} - -void -MultiSigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) -{ - if (dt == dt_B) return; - - dt_B = dt; - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(*this); mfi.isValid(); ++mfi) - { - (*this)[mfi].ComputePMLFactorsB(dx, dt); - } -} - -void -MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) -{ - if (dt == dt_E) return; - - dt_E = dt; - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(*this); mfi.isValid(); ++mfi) - { - (*this)[mfi].ComputePMLFactorsE(dx, dt); - } -} - -PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, - const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window) - : m_geom(geom), - m_cgeom(cgeom) -{ - const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); - if (ba.size() == 0) { - m_ok = false; - return; - } else { - m_ok = true; - } - - DistributionMapping dm{ba}; - - int nge = 2; - int ngb = 2; - int ngf = (do_moving_window) ? 2 : 0; - if (WarpX::maxwell_fdtd_solver_id == 1) ngf = std::max( ngf, 1 ); - - pml_E_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::Ex_nodal_flag), dm, 3, nge)); - pml_E_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::Ey_nodal_flag), dm, 3, nge)); - pml_E_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Ez_nodal_flag), dm, 3, nge)); - pml_B_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::Bx_nodal_flag), dm, 2, ngb)); - pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb)); - pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb)); - - pml_E_fp[0]->setVal(0.0); - pml_E_fp[1]->setVal(0.0); - pml_E_fp[2]->setVal(0.0); - pml_B_fp[0]->setVal(0.0); - pml_B_fp[1]->setVal(0.0); - pml_B_fp[2]->setVal(0.0); - - if (do_dive_cleaning) - { - pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf)); - pml_F_fp->setVal(0.0); - } - - sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); - - if (cgeom) - { - - nge = 1; - ngb = 1; - - BoxArray grid_cba = grid_ba; - grid_cba.coarsen(ref_ratio); - const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell); - - DistributionMapping cdm{cba}; - - pml_E_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::Ex_nodal_flag), cdm, 3, nge)); - pml_E_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::Ey_nodal_flag), cdm, 3, nge)); - pml_E_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::Ez_nodal_flag), cdm, 3, nge)); - pml_B_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::Bx_nodal_flag), cdm, 2, ngb)); - pml_B_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::By_nodal_flag), cdm, 2, ngb)); - pml_B_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::Bz_nodal_flag), cdm, 2, ngb)); - - pml_E_cp[0]->setVal(0.0); - pml_E_cp[1]->setVal(0.0); - pml_E_cp[2]->setVal(0.0); - pml_B_cp[0]->setVal(0.0); - pml_B_cp[1]->setVal(0.0); - pml_B_cp[2]->setVal(0.0); - - if (do_dive_cleaning) - { - pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf)); - pml_F_cp->setVal(0.0); - } - - sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); - } - -} - -// BoxArray -// PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) -// { -// Box domain = geom.Domain(); -// for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { -// if ( ! geom.isPeriodic(idim) ) { -// domain.grow(idim, ncell); -// } -// } -// -// BoxList bl; -// for (int i = 0, N = grid_ba.size(); i < N; ++i) -// { -// const Box& grid_bx = grid_ba[i]; -// const IntVect& grid_bx_sz = grid_bx.size(); -// AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, -// "Consider using larger amr.blocking_factor"); -// -// Box bx = grid_bx; -// bx.grow(ncell); -// bx &= domain; -// -// Vector bndryboxes; -// #if (AMREX_SPACEDIM == 3) -// int kbegin = -1, kend = 1; -// #else -// int kbegin = 0, kend = 0; -// #endif -// for (int kk = kbegin; kk <= kend; ++kk) { -// for (int jj = -1; jj <= 1; ++jj) { -// for (int ii = -1; ii <= 1; ++ii) { -// if (ii != 0 || jj != 0 || kk != 0) { -// Box b = grid_bx; -// b.shift(grid_bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); -// b &= bx; -// if (b.ok()) { -// bndryboxes.push_back(b); -// } -// } -// } -// } -// } -// -// const BoxList& noncovered = grid_ba.complementIn(bx); -// for (const Box& b : noncovered) { -// for (const auto& bb : bndryboxes) { -// Box ib = b & bb; -// if (ib.ok()) { -// bl.push_back(ib); -// } -// } -// } -// } -// -// BoxArray ba(bl); -// ba.removeOverlap(false); -// -// return ba; -// } - -BoxArray -PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) -{ - Box domain = geom.Domain(); - - IntVect limInfDomain = domain.smallEnd(); - IntVect limSupDomain = domain.bigEnd(); - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if ( ! geom.isPeriodic(idim) ) { - // Create the simulation area (inside the whole domain) (only if non-periodical boundary conditions) - domain.grow(idim, -ncell); // I don't know if that works - // domain.growHi(idim, -ncell); - // domain.growLo(idim, -ncell); - } - } - - BoxList bl; - for (int i = 0, N = grid_ba.size(); i < N; ++i) - { - const Box& grid_bx = grid_ba[i]; - // const IntVect& grid_bx_sz = grid_bx.size(); - // AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, - // "Consider using larger amr.blocking_factor"); - - Box bx = grid_bx; - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - bx.grow(idim, -ncell); - } - bx &= domain; // is this step necessary? We're always inside the domain by doing so, even if there are some periodical boundary conditions... - const IntVect& bx_sz = bx.size(); - - Vector bndryboxes; -#if (AMREX_SPACEDIM == 3) - int kbegin = -1, kend = 1; -#else - int kbegin = 0, kend = 0; -#endif - for (int kk = kbegin; kk <= kend; ++kk) { - for (int jj = -1; jj <= 1; ++jj) { - for (int ii = -1; ii <= 1; ++ii) { - if (ii != 0 || jj != 0 || kk != 0) { - // if ((ii != 0 && jj == 0 && kk == 0)||(ii == 0 && jj != 0 && kk == 0)||(ii == 0 && jj == 0 && kk != 0)) { - if(ii!=0 && jj!=0){ - int xlimg, xlimd, ylimg, ylimd; - if (ii == -1){ - xlimg = grid_bx.smallEnd()[0]; - xlimd = limInfDomain[0]; - } - else if (ii = 1) { - xlimg = grid_bx.bigEnd()[0]; - xlimd = limSupDomain[0]; - } - - if (jj == -1){ - ylimg = grid_bx.smallEnd()[1]; - ylimd = limInfDomain[1]; - } - else if (jj = 1){ - ylimg = grid_bx.bigEnd()[1]; - ylimd = limSupDomain[1]; - } - - if (xlimd==xlimg && ylimd==ylimg){ - amrex::Print() << "INDICES : i = " << ii << " j = " << jj << std::endl; - Box b = grid_bx; //grid_bx; - b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); - b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); - b &= grid_bx;// grid_bx; - amrex::Print() << "BOITE = [" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<ComputePMLFactorsB(m_geom->CellSize(), dt); - sigba_fp->ComputePMLFactorsE(m_geom->CellSize(), dt); - } - if (sigba_cp) { - sigba_cp->ComputePMLFactorsB(m_cgeom->CellSize(), dt); - sigba_cp->ComputePMLFactorsE(m_cgeom->CellSize(), dt); - } -} - -std::array -PML::GetE_fp () -{ - return {pml_E_fp[0].get(), pml_E_fp[1].get(), pml_E_fp[2].get()}; -} - -std::array -PML::GetB_fp () -{ - return {pml_B_fp[0].get(), pml_B_fp[1].get(), pml_B_fp[2].get()}; -} - -std::array -PML::GetE_cp () -{ - return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()}; -} - -std::array -PML::GetB_cp () -{ - return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()}; -} - -MultiFab* -PML::GetF_fp () -{ - return pml_F_fp.get(); -} - -MultiFab* -PML::GetF_cp () -{ - return pml_F_cp.get(); -} - -void -PML::ExchangeB (const std::array& B_fp, - const std::array& B_cp) -{ - ExchangeB(PatchType::fine, B_fp); - ExchangeB(PatchType::coarse, B_cp); -} - -void -PML::ExchangeB (PatchType patch_type, - const std::array& Bp) -{ - if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) - { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom); - } - else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) - { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); - } -} - -void -PML::ExchangeE (const std::array& E_fp, - const std::array& E_cp) -{ - ExchangeE(PatchType::fine, E_fp); - ExchangeE(PatchType::coarse, E_cp); -} - -void -PML::ExchangeE (PatchType patch_type, - const std::array& Ep) -{ - if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) - { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom); - } - else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) - { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); - } -} - -void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) -{ - ExchangeF(PatchType::fine, F_fp); - ExchangeF(PatchType::coarse, F_cp); -} - -void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp) -{ - if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom); - } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom); - } -} - -void -PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) -{ - const IntVect& ngr = reg.nGrowVect(); - const IntVect& ngp = pml.nGrowVect(); - const int ncp = pml.nComp(); - const auto& period = geom.periodicity(); - - MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr); - - if (ngp.max() > 0) // Copy from pml to the ghost cells of regular data - { - MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); - MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); - if (ncp == 3) { - MultiFab::Add(totpmlmf,pml,2,0,1,0); - } - - MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); - tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(reg); mfi.isValid(); ++mfi) - { - const FArrayBox& src = tmpregmf[mfi]; - FArrayBox& dst = reg[mfi]; - const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); - for (const Box& bx : bl) - { - dst.copy(src, bx, 0, bx, 0, 1); - } - } - } - - // Copy from regular data to PML's first component - // Zero out the second (and third) component - MultiFab::Copy(tmpregmf,reg,0,0,1,0); - tmpregmf.setVal(0.0, 1, ncp-1, 0); - pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period); -} - -void -PML::FillBoundary () -{ - FillBoundaryE(); - FillBoundaryB(); - FillBoundaryF(); -} - -void -PML::FillBoundaryE () -{ - FillBoundaryE(PatchType::fine); - FillBoundaryE(PatchType::coarse); -} - -void -PML::FillBoundaryE (PatchType patch_type) -{ - if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0) - { - const auto& period = m_geom->periodicity(); - Vector mf{pml_E_fp[0].get(),pml_E_fp[1].get(),pml_E_fp[2].get()}; - amrex::FillBoundary(mf, period); - } - else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) - { - const auto& period = m_cgeom->periodicity(); - Vector mf{pml_E_cp[0].get(),pml_E_cp[1].get(),pml_E_cp[2].get()}; - amrex::FillBoundary(mf, period); - } -} - -void -PML::FillBoundaryB () -{ - FillBoundaryB(PatchType::fine); - FillBoundaryB(PatchType::coarse); -} - -void -PML::FillBoundaryB (PatchType patch_type) -{ - if (patch_type == PatchType::fine && pml_B_fp[0]) - { - const auto& period = m_geom->periodicity(); - Vector mf{pml_B_fp[0].get(),pml_B_fp[1].get(),pml_B_fp[2].get()}; - amrex::FillBoundary(mf, period); - } - else if (patch_type == PatchType::coarse && pml_B_cp[0]) - { - const auto& period = m_cgeom->periodicity(); - Vector mf{pml_B_cp[0].get(),pml_B_cp[1].get(),pml_B_cp[2].get()}; - amrex::FillBoundary(mf, period); - } -} - -void -PML::FillBoundaryF () -{ - FillBoundaryF(PatchType::fine); - FillBoundaryF(PatchType::coarse); -} - -void -PML::FillBoundaryF (PatchType patch_type) -{ - if (patch_type == PatchType::fine && pml_F_fp && pml_F_fp->nGrowVect().max() > 0) - { - const auto& period = m_geom->periodicity(); - pml_F_fp->FillBoundary(period); - } - else if (patch_type == PatchType::coarse && pml_F_cp && pml_F_cp->nGrowVect().max() > 0) - { - const auto& period = m_cgeom->periodicity(); - pml_F_cp->FillBoundary(period); - } -} - -void -PML::CheckPoint (const std::string& dir) const -{ - if (pml_E_fp[0]) - { - VisMF::Write(*pml_E_fp[0], dir+"_Ex_fp"); - VisMF::Write(*pml_E_fp[1], dir+"_Ey_fp"); - VisMF::Write(*pml_E_fp[2], dir+"_Ez_fp"); - VisMF::Write(*pml_B_fp[0], dir+"_Bx_fp"); - VisMF::Write(*pml_B_fp[1], dir+"_By_fp"); - VisMF::Write(*pml_B_fp[2], dir+"_Bz_fp"); - } - - if (pml_E_cp[0]) - { - VisMF::Write(*pml_E_cp[0], dir+"_Ex_cp"); - VisMF::Write(*pml_E_cp[1], dir+"_Ey_cp"); - VisMF::Write(*pml_E_cp[2], dir+"_Ez_cp"); - VisMF::Write(*pml_B_cp[0], dir+"_Bx_cp"); - VisMF::Write(*pml_B_cp[1], dir+"_By_cp"); - VisMF::Write(*pml_B_cp[2], dir+"_Bz_cp"); - } -} - -void -PML::Restart (const std::string& dir) -{ - if (pml_E_fp[0]) - { - VisMF::Read(*pml_E_fp[0], dir+"_Ex_fp"); - VisMF::Read(*pml_E_fp[1], dir+"_Ey_fp"); - VisMF::Read(*pml_E_fp[2], dir+"_Ez_fp"); - VisMF::Read(*pml_B_fp[0], dir+"_Bx_fp"); - VisMF::Read(*pml_B_fp[1], dir+"_By_fp"); - VisMF::Read(*pml_B_fp[2], dir+"_Bz_fp"); - } - - if (pml_E_cp[0]) - { - VisMF::Read(*pml_E_cp[0], dir+"_Ex_cp"); - VisMF::Read(*pml_E_cp[1], dir+"_Ey_cp"); - VisMF::Read(*pml_E_cp[2], dir+"_Ez_cp"); - VisMF::Read(*pml_B_cp[0], dir+"_Bx_cp"); - VisMF::Read(*pml_B_cp[1], dir+"_By_cp"); - VisMF::Read(*pml_B_cp[2], dir+"_Bz_cp"); - } -} diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index 08173e424..dd67bfee5 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -111,21 +111,19 @@ WarpX::DampJPML (int lev, PatchType patch_type) #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - // amrex::Print()<<"===== DAMP PML ====="<ok()){ @@ -321,14 +317,7 @@ WarpX::OneStep_nosub (Real cur_time) } if (do_pml && do_pml_j_damping){ - // amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<ok()){ - // amrex::Print()<< "[AV COPY] max_Jx = "<< current_fp[lev][0].get()->min(0) << std::endl; - // amrex::Print()<< "[AV COPY] max_Jx_pml = "<< pml[lev]->Getj_fp()[0]->min(0) << std::endl; - // amrex::Print()<< "===== Copy J from PML to Reg ====="<< std::endl; pml[lev]->CopyJinReg({ current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get() }, { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }); - // amrex::Print()<< "[AP COPY] max_Jx = "<< current_fp[lev][0].get()->min(0) << std::endl; - } } } EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { - // amrex::Print()<< "WarpXEvolveEM.cpp : Damp "< > pml; amrex::Real moving_window_x = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 77fd83f7a..08a8b7b88 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -382,6 +382,11 @@ WarpX::ReadParameters () pp.query("pml_delta", pml_delta); pp.query("pml_has_particles", pml_has_particles); pp.query("do_pml_j_damping", do_pml_j_damping); + pp.query("do_pml_in_domain", do_pml_in_domain); + + if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ + amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)"); + } pp.query("dump_openpmd", dump_openpmd); pp.query("dump_plotfiles", dump_plotfiles); -- cgit v1.2.3 From 065b2d9eaaaf1dc6edfc62cfbc01291178426e9f Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 25 Jul 2019 01:45:29 +0200 Subject: [WIP] Adding flags do_pml_Lo, do_pml_hi --- Source/BoundaryConditions/PML.H | 14 +++---- Source/BoundaryConditions/PML.cpp | 59 ++++++++++++++++-------------- Source/BoundaryConditions/PML_routines.F90 | 2 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- Source/Parallelization/WarpXComm.cpp | 18 ++++----- Source/WarpX.H | 4 +- Source/WarpX.cpp | 15 ++++++++ 7 files changed, 67 insertions(+), 47 deletions(-) (limited to 'Source/WarpX.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index c64eb8193..3f64db558 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -149,25 +149,25 @@ public: { return *sigba_cp; } void ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain); + const std::array& B_cp, int do_pml_in_domain, int ncell); void ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain); + const std::array& E_cp, int do_pml_in_domain, int ncell); void CopyJinPMLs (const std::array& j_fp, const std::array& j_cp); void CopyJinReg (const std::array& j_fp, const std::array& j_cp); void ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain); + const std::array& Bp, int do_pml_in_domain, int ncell); void ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain); + const std::array& Ep, int do_pml_in_domain, int ncell); void CopyJinPMLs (PatchType patch_type, const std::array& jp); void CopyJinReg (PatchType patch_type, const std::array& jp); - void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain); - void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain); + void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain, int ncell); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain, int ncell); void FillBoundary (); void FillBoundaryE (); @@ -205,7 +205,7 @@ private: static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); - static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain); + static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain, int ncell); static void CopyRegInPMLs (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); static void CopyPMLsInReg (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 51107b498..0e6c24131 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -708,53 +708,53 @@ PML::GetF_cp () void PML::ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain) + const std::array& B_cp, int do_pml_in_domain, int ncell) { - ExchangeB(PatchType::fine, B_fp, do_pml_in_domain); - ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain); + ExchangeB(PatchType::fine, B_fp, do_pml_in_domain, ncell); + ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain, ncell); } void PML::ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain) + const std::array& Bp, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain, ncell); } } void PML::ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain) + const std::array& E_cp, int do_pml_in_domain, int ncell) { - ExchangeE(PatchType::fine, E_fp, do_pml_in_domain); - ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain); + ExchangeE(PatchType::fine, E_fp, do_pml_in_domain, ncell); + ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain, ncell); } void PML::ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain) + const std::array& Ep, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain, ncell); } } @@ -811,29 +811,32 @@ PML::CopyJinReg (const std::array& j_fp, } void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) +PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain, int ncell) { - ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); - ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); + ExchangeF(PatchType::fine, F_fp, do_pml_in_domain, ncell); + ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain, ncell); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) +PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); + Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain); + Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain, ncell); } } void -PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, int do_pml_in_domain) +PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, int do_pml_in_domain, int ncell) { if (do_pml_in_domain){ const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); + // amrex::Print()<< "##### PML EXCHANGE #####"<ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain); + if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain, pml_ncell); const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 7c00a5297..17d87c73f 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,7 +18,7 @@ WarpX::ExchangeWithPmlB (int lev) Bfield_fp[lev][2].get() }, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, do_pml_in_domain); + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); } } @@ -31,7 +31,7 @@ WarpX::ExchangeWithPmlE (int lev) Efield_fp[lev][2].get() }, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }, do_pml_in_domain); + Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); } } @@ -40,7 +40,7 @@ WarpX::ExchangeWithPmlF (int lev) { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), - F_cp[lev].get(), do_pml_in_domain); + F_cp[lev].get(), do_pml_in_domain, pml_ncell); } } @@ -252,7 +252,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }, do_pml_in_domain); + Efield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryE(patch_type); } @@ -267,7 +267,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }, do_pml_in_domain); + Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryE(patch_type); } @@ -294,7 +294,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }, do_pml_in_domain); + Bfield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryB(patch_type); } const auto& period = Geom(lev).periodicity(); @@ -308,7 +308,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, do_pml_in_domain); + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryB(patch_type); } const auto& cperiod = Geom(lev-1).periodicity(); @@ -331,7 +331,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain); + pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryF(patch_type); } @@ -342,7 +342,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain); + pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryF(patch_type); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 54cad9ae0..20db5c75d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -489,10 +489,12 @@ private: // PML int do_pml = 1; int pml_ncell = 10; - int pml_delta = 10; + int pml_delta = 5; int pml_has_particles = 0; int do_pml_j_damping = 0; int do_pml_in_domain = 0; + amrex::IntVect do_pml_Lo; + amrex::IntVect do_pml_Hi; amrex::Vector > pml; amrex::Real moving_window_x = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 08a8b7b88..46b1926a0 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -384,6 +384,21 @@ WarpX::ReadParameters () pp.query("do_pml_j_damping", do_pml_j_damping); pp.query("do_pml_in_domain", do_pml_in_domain); + Vector parse_do_pml_Lo(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Lo", parse_do_pml_Lo); + do_pml_Lo[0] = parse_do_pml_Lo[0]; + do_pml_Lo[1] = parse_do_pml_Lo[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Lo[2] = parse_do_pml_Lo[2]; +#endif + Vector parse_do_pml_Hi(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Hi", parse_do_pml_Hi); + do_pml_Hi[0] = parse_do_pml_Hi[0]; + do_pml_Hi[1] = parse_do_pml_Hi[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Hi[2] = parse_do_pml_Hi[2]; +#endif + if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)"); } -- cgit v1.2.3 From 7375c64a9a8c4ff2782392f04a423e4519f0095a Mon Sep 17 00:00:00 2001 From: ablelly Date: Mon, 5 Aug 2019 17:11:49 +0200 Subject: Added load balancing at first iteration only. --- Source/Evolve/WarpXEvolveEM.cpp | 2 +- Source/WarpX.H | 1 + Source/WarpX.cpp | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) (limited to 'Source/WarpX.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index cd5388ec9..ce27c67ff 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -54,7 +54,7 @@ WarpX::EvolveEM (int numsteps) amrex::Abort("LoadBalance for PSATD: TODO"); #endif - if (step > 0 && (step+1) % load_balance_int == 0) + if (step > 0 && step < load_balance_max_step && (step+1) % load_balance_int == 0) { LoadBalance(); // Reset the costs to 0 diff --git a/Source/WarpX.H b/Source/WarpX.H index 24fd61914..2512285ef 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -511,6 +511,7 @@ private: amrex::Real const_dt = 0.5e-11; int load_balance_int = -1; + int load_balance_max_step = std::numeric_limits::max(); amrex::Vector > costs; int load_balance_with_sfc = 0; amrex::Real load_balance_knapsack_factor = 1.24; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 46b1926a0..ea2af4330 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -473,6 +473,7 @@ WarpX::ReadParameters () } pp.query("load_balance_int", load_balance_int); + pp.query("load_balance_max_step", load_balance_max_step); pp.query("load_balance_with_sfc", load_balance_with_sfc); pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); -- cgit v1.2.3 From 7a5adaa22b1e1830b2046aba6ba38de90c96138c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 21 Aug 2019 13:13:23 -0700 Subject: Code cleanup + add test --- Examples/Tests/particles_in_PML/analysis_2d.py | 26 ++++++ Examples/Tests/particles_in_PML/inputs2d | 56 ++++++++++++ Source/BoundaryConditions/PML.cpp | 115 ++----------------------- Source/Evolve/WarpXEvolveEM.cpp | 18 +--- Source/FieldSolver/WarpXPushFieldsEM.cpp | 1 - Source/Parallelization/WarpXComm.cpp | 12 +-- Source/Particles/PhysicalParticleContainer.cpp | 32 +++---- Source/WarpX.H | 1 - Source/WarpX.cpp | 1 - 9 files changed, 113 insertions(+), 149 deletions(-) create mode 100755 Examples/Tests/particles_in_PML/analysis_2d.py create mode 100644 Examples/Tests/particles_in_PML/inputs2d (limited to 'Source/WarpX.cpp') diff --git a/Examples/Tests/particles_in_PML/analysis_2d.py b/Examples/Tests/particles_in_PML/analysis_2d.py new file mode 100755 index 000000000..5b0279a9b --- /dev/null +++ b/Examples/Tests/particles_in_PML/analysis_2d.py @@ -0,0 +1,26 @@ +""" +This script tests the absorption of particles in the PML. + +The input file inputs_2d is used: it features a positive and a +negative particle, going in opposite direction and eventually +leaving the box. This script tests that the field in the box +is close to 0 once the particles have left. With regular +PML, this test fails, since the particles leave a spurious +charge, with associated fields, behind them. +""" +import sys +import yt +import numpy as np +yt.funcs.mylog.setLevel(0) + +# Open plotfile specified in command line +filename = sys.argv[1] +ds = yt.load( filename ) + +# Check that the field is low enough +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray() +Ey_array = ad0['Ey'].to_ndarray() +Ez_array = ad0['Ez'].to_ndarray() +max_Efield = max(Ex_array.max(), Ey_array.max(), Ez_array.max()) +assert max_Efield < 0.0003 \ No newline at end of file diff --git a/Examples/Tests/particles_in_PML/inputs2d b/Examples/Tests/particles_in_PML/inputs2d new file mode 100644 index 000000000..733032c53 --- /dev/null +++ b/Examples/Tests/particles_in_PML/inputs2d @@ -0,0 +1,56 @@ +max_step = 200 +amr.plot_int = 30 +amr.n_cell = 224 224 + +amr.blocking_factor = 8 +amr.max_grid_size = 1024 +amr.max_level = 0 + +warpx.dump_plotfiles = 1 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 + +# PML +warpx.do_pml = 1 +warpx.pml_ncell = 12 +warpx.pml_delta = 6 +warpx.pml_has_particles = 0 +warpx.do_pml_in_domain = 0 +warpx.do_pml_j_damping = 0 + + +# Algorithms +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = vectorized +algo.particle_pusher = vay +algo.maxwell_fdtd_solver = ckc +warpx.cfl = 1.0 +warpx.use_filter = 1 + +# Particle species +particles.nspecies = 2 +particles.species_names = electron proton + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "singleparticle" +electron.single_particle_pos = 0. 0. 0. +electron.single_particle_vel = 2. 0. 0. +electron.single_particle_weight = 1. + +proton.charge = q_e +proton.mass = m_p # Very heavy ; should not move +proton.injection_style = "singleparticle" +proton.single_particle_pos = 0. 0. 0. +proton.single_particle_vel = -2. 0. 0. +proton.single_particle_weight = 1. + +# Particle shape factor in each direction +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 0b73859b4..15f18201b 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -24,22 +24,19 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real x = 10.0; - Real theta = 10.0; - Real coeff_damp = 1.; //std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast(glo-i); sigma[i-slo] = fac*(offset*offset); - sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_cum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast(glo-i) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); - sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_star_cum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } } @@ -52,20 +49,17 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real x = 10.0; - Real theta = 10.0; - Real coeff_damp = 1.;//std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast(i-ghi-1); sigma[i-slo] = fac*(offset*offset); - sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/PhysConst::c; } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast(i-ghi) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); - sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/PhysConst::c; } } @@ -119,8 +113,6 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n sigma_star_cum_fac[idim].m_hi = hi[idim]; } - amrex::Print()< 0){pml_type[2]=idim+'0';} //std::to_string(idim) + if (idim > 0){pml_type[2]=idim+'0';} } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox & box; if (hioverlap.ok()) { - amrex::Print()<<"["< 1, Box gaps not wide enough?\n"); } } - amrex::Print()< 0) // Copy from pml to the ghost cells of regular data { - MultiFab::Copy(tmpregmf, reg, 0, 0, 1, IntVect(0)); //ngr); + MultiFab::Copy(tmpregmf, reg, 0, 0, 1, IntVect(0)); MultiFab::Copy(totpmlmf, pml, 0, 0, ncp, ngp); tmpregmf.setVal(0.0, 1, ncp-1, 0); - // ce qui foncitonne a peu pres pour le moment tmpregmf.ParallelCopy(totpmlmf,0, 0, ncp, IntVect(0), IntVect(0), period); totpmlmf.ParallelCopy(tmpregmf,0, 0, ncp, IntVect(0), ngp, period); - MultiFab::Copy(pml, totpmlmf, 0, 0, ncp, ngp); //test -// #ifdef _OPENMP -// #pragma omp parallel -// #endif - // // amrex::Print()<<"####### EXCHANGE INFORMATIONS GHOST CELLS PML #######"<>> EXCHANGE BOXES <<<<<" <>> EXCHANGE BOXES <<<<<" < 0 && step < load_balance_max_step && (step+1) % load_balance_int == 0) + if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); // Reset the costs to 0 @@ -329,22 +329,6 @@ WarpX::OneStep_nosub (Real cur_time) EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); EvolveF(0.5*dt[0], DtType::SecondHalf); - - - // if (do_pml && pml_has_particles){ - // for (int lev = 0; lev <= finest_level; ++lev) - // { - // if (pml[lev]->ok()){ - // pml[lev]->CopyJinReg({ current_fp[lev][0].get(), - // current_fp[lev][1].get(), - // current_fp[lev][2].get() }, - // { current_cp[lev][0].get(), - // current_cp[lev][1].get(), - // current_cp[lev][2].get() }); - // } - // } - // } - EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { DampPML(); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 18c98b8f8..011d9cfcc 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -241,7 +241,6 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - // const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 06cbf157a..c539cd89a 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -256,7 +256,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_fp[lev][0].get(), + { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), Efield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -273,7 +273,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_cp[lev][0].get(), + { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -302,7 +302,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_fp[lev][0].get(), + { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -318,9 +318,9 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_cp[lev][0].get(), - Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, + { Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell, do_pml_Lo, do_pml_Hi); pml[lev]->FillBoundaryB(patch_type); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 20819a78c..d10390204 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -33,7 +33,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -42,7 +42,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -59,9 +59,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -130,7 +130,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -203,7 +203,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); - + particle_tile.push_back_real(particle_comps["uxold"], u[0]); particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); @@ -931,11 +931,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PICSAR::FieldGather", blp_pxr_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -949,7 +949,7 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1227,13 +1227,13 @@ PhysicalParticleContainer::Evolve (int lev, eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); ceyfab = &filtered_Ey; - + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); cbxfab = &filtered_Bx; - + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); @@ -1453,7 +1453,7 @@ PhysicalParticleContainer::SplitParticles(int lev) } // Add local arrays psplit_x etc. to the temporary // particle container pctmp_split. Split particles - // are tagged with p.id()=NoSplitParticleID so that + // are tagged with p.id()=NoSplitParticleID so that // they are not re-split when entering a higher level // AddNParticles calls Redistribute, so that particles // in pctmp_split are in the proper grids and tiles @@ -1550,7 +1550,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1694,7 +1694,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1776,7 +1776,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); diff --git a/Source/WarpX.H b/Source/WarpX.H index b4dfe9ae4..3e8d45f49 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -524,7 +524,6 @@ private: amrex::Real const_dt = 0.5e-11; int load_balance_int = -1; - int load_balance_max_step = std::numeric_limits::max(); amrex::Vector > costs; int load_balance_with_sfc = 0; amrex::Real load_balance_knapsack_factor = 1.24; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6abf27216..305fedbb4 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -498,7 +498,6 @@ WarpX::ReadParameters () } pp.query("load_balance_int", load_balance_int); - pp.query("load_balance_max_step", load_balance_max_step); pp.query("load_balance_with_sfc", load_balance_with_sfc); pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); -- cgit v1.2.3