diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/PhysicalParticleContainer.cpp | 8 | ||||
| -rw-r--r-- | Source/RigidInjectedParticleContainer.cpp | 2 | ||||
| -rw-r--r-- | Source/WarpX.H | 8 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 56 | ||||
| -rw-r--r-- | Source/WarpXIO.cpp | 173 | ||||
| -rw-r--r-- | Source/WarpXParticleContainer.cpp | 35 | ||||
| -rw-r--r-- | Source/WarpX_f.H | 2 | ||||
| -rw-r--r-- | Source/WarpX_picsar.F90 | 11 |
8 files changed, 189 insertions, 106 deletions
diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 564724f49..85644760c 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -983,7 +983,7 @@ PhysicalParticleContainer::FieldGather (int lev, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); if (cost) { @@ -1288,7 +1288,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*bxfab), BL_TO_FORTRAN_ANYD(*byfab), BL_TO_FORTRAN_ANYD(*bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); if (np_gather < np) @@ -1385,7 +1385,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*cbxfab), BL_TO_FORTRAN_ANYD(*cbyfab), BL_TO_FORTRAN_ANYD(*cbzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); } @@ -1560,7 +1560,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); warpx_particle_pusher_momenta(&np, diff --git a/Source/RigidInjectedParticleContainer.cpp b/Source/RigidInjectedParticleContainer.cpp index eb2c1c4a8..db3623705 100644 --- a/Source/RigidInjectedParticleContainer.cpp +++ b/Source/RigidInjectedParticleContainer.cpp @@ -427,7 +427,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &l_lower_order_in_v, + &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); // Save the position and momenta, making copies diff --git a/Source/WarpX.H b/Source/WarpX.H index 08fe657b4..bf82bceca 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -114,6 +114,9 @@ public: static int n_field_gather_buffer; static int n_current_deposition_buffer; + // do nodal + static int do_nodal; + const amrex::MultiFab& getcurrent (int lev, int direction) {return *current_fp[lev][direction];} const amrex::MultiFab& getEfield (int lev, int direction) {return *Efield_aux[lev][direction];} const amrex::MultiFab& getBfield (int lev, int direction) {return *Bfield_aux[lev][direction];} @@ -389,6 +392,10 @@ private: return gather_buffer_masks[lev].get(); } + void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, + const amrex::IntVect& ngE, const amrex::IntVect& ngJ, + const amrex::IntVect& ngRho, int ngF); + amrex::Vector<int> istep; // which step? amrex::Vector<int> nsubsteps; // how many substeps on each level? @@ -443,7 +450,6 @@ private: amrex::Vector<std::unique_ptr<amrex::iMultiFab> > rho_fp_owner_masks; amrex::Vector<std::unique_ptr<amrex::iMultiFab> > rho_cp_owner_masks; - // div E cleaning int do_dive_cleaning = 0; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 1828ebcec..ba54acee0 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -94,6 +94,8 @@ IntVect WarpX::jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX int WarpX::n_field_gather_buffer = 0; int WarpX::n_current_deposition_buffer = -1; +int WarpX::do_nodal = false; + WarpX* WarpX::m_instance = nullptr; @@ -410,6 +412,21 @@ WarpX::ReadParameters () pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); pp.query("do_dynamic_scheduling", do_dynamic_scheduling); + + pp.query("do_nodal", do_nodal); + if (do_nodal) { + Bx_nodal_flag = IntVect::TheNodeVector(); + By_nodal_flag = IntVect::TheNodeVector(); + Bz_nodal_flag = IntVect::TheNodeVector(); + Ex_nodal_flag = IntVect::TheNodeVector(); + Ey_nodal_flag = IntVect::TheNodeVector(); + Ez_nodal_flag = IntVect::TheNodeVector(); + jx_nodal_flag = IntVect::TheNodeVector(); + jy_nodal_flag = IntVect::TheNodeVector(); + jz_nodal_flag = IntVect::TheNodeVector(); + // Use same shape factors in all directions, for gathering + l_lower_order_in_v = false; + } } { @@ -430,22 +447,22 @@ WarpX::ReadParameters () pp.query("particle_pusher", particle_pusher_algo); std::string s_solver = ""; pp.query("maxwell_fdtd_solver", s_solver); - std::transform(s_solver.begin(), - s_solver.end(), - s_solver.begin(), - ::tolower); + std::transform(s_solver.begin(), + s_solver.end(), + s_solver.begin(), + ::tolower); // if maxwell_fdtd_solver is specified, set the value // of maxwell_fdtd_solver_id accordingly. - // Otherwise keep the default value maxwell_fdtd_solver_id=0 - if (s_solver != "") { - if (s_solver == "yee") { - maxwell_fdtd_solver_id = 0; - } else if (s_solver == "ckc") { - maxwell_fdtd_solver_id = 1; - } else { - amrex::Abort("Unknown FDTD Solver type " + s_solver); + // Otherwise keep the default value maxwell_fdtd_solver_id=0 + if (s_solver != "") { + if (s_solver == "yee") { + maxwell_fdtd_solver_id = 0; + } else if (s_solver == "ckc") { + maxwell_fdtd_solver_id = 1; + } else { + amrex::Abort("Unknown FDTD Solver type " + s_solver); + } } - } } #ifdef WARPX_USE_PSATD @@ -619,6 +636,13 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d // CKC solver requires one additional guard cell if (maxwell_fdtd_solver_id == 1) ngF = std::max( ngF, 1 ); + AllocLevelMFs(lev, ba, dm, ngE, ngJ, ngRho, ngF); +} + +void +WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, + const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF) +{ // // The fine patch // @@ -638,7 +662,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_fp_owner_masks[lev][0] = std::move(current_fp[lev][0]->OwnerMask(period)); current_fp_owner_masks[lev][1] = std::move(current_fp[lev][1]->OwnerMask(period)); current_fp_owner_masks[lev][2] = std::move(current_fp[lev][2]->OwnerMask(period)); - + if (do_dive_cleaning || plot_rho) { rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); @@ -708,7 +732,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); - const auto& cperiod = Geom(lev).periodicity(); + const auto& cperiod = Geom(lev-1).periodicity(); current_cp_owner_masks[lev][0] = std::move(current_cp[lev][0]->OwnerMask(cperiod)); current_cp_owner_masks[lev][1] = std::move(current_cp[lev][1]->OwnerMask(cperiod)); current_cp_owner_masks[lev][2] = std::move(current_cp[lev][2]->OwnerMask(cperiod)); @@ -770,6 +794,8 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d if (load_balance_int > 0) { costs[lev].reset(new MultiFab(ba, dm, 1, 0)); } + + } std::array<Real,3> diff --git a/Source/WarpXIO.cpp b/Source/WarpXIO.cpp index 93fb90b8f..e68b4a4e8 100644 --- a/Source/WarpXIO.cpp +++ b/Source/WarpXIO.cpp @@ -782,13 +782,24 @@ WarpX::WritePlotFile () const mf[lev].reset(new MultiFab(grids[lev], dmap[lev], ncomp, ngrow)); Vector<const MultiFab*> srcmf(AMREX_SPACEDIM); - PackPlotDataPtrs(srcmf, current_fp[lev]); int dcomp = 0; - amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); + + // j + if (do_nodal) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *current_fp[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *current_fp[lev][1], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *current_fp[lev][2], 0, 1); + } + else + { + PackPlotDataPtrs(srcmf, current_fp[lev]); + amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *current_fp[lev][1], 0, 1); + MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *current_fp[lev][1], 0, 1); #endif + } if (lev == 0) { varnames.push_back("jx"); @@ -797,12 +808,22 @@ WarpX::WritePlotFile () const } dcomp += 3; - PackPlotDataPtrs(srcmf, Efield_aux[lev]); - amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); + // E + if (do_nodal) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *Efield_aux[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_aux[lev][1], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *Efield_aux[lev][2], 0, 1); + } + else + { + PackPlotDataPtrs(srcmf, Efield_aux[lev]); + amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_aux[lev][1], 0, 1); + MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_aux[lev][1], 0, 1); #endif + } if (lev == 0) { varnames.push_back("Ex"); @@ -811,12 +832,22 @@ WarpX::WritePlotFile () const } dcomp += 3; - PackPlotDataPtrs(srcmf, Bfield_aux[lev]); - amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); + // B + if (do_nodal) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *Bfield_aux[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Bfield_aux[lev][1], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *Bfield_aux[lev][2], 0, 1); + } + else + { + PackPlotDataPtrs(srcmf, Bfield_aux[lev]); + amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - MultiFab::Copy(*mf[lev], *Bfield_aux[lev][1], 0, dcomp+1, 1, ngrow); + MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); + MultiFab::Copy(*mf[lev], *Bfield_aux[lev][1], 0, dcomp+1, 1, ngrow); #endif + } if (lev == 0) { varnames.push_back("Bx"); @@ -897,6 +928,7 @@ WarpX::WritePlotFile () const if (plot_divb) { + if (do_nodal) amrex::Abort("TODO: do_nodal && plot_divb"); ComputeDivB(*mf[lev], dcomp, {Bfield_aux[lev][0].get(),Bfield_aux[lev][1].get(),Bfield_aux[lev][2].get()}, WarpX::CellSize(lev)); @@ -909,6 +941,7 @@ WarpX::WritePlotFile () const if (plot_dive) { + if (do_nodal) amrex::Abort("TODO: do_nodal && plot_dive"); const BoxArray& ba = amrex::convert(boxArray(lev),IntVect::TheUnitVector()); MultiFab dive(ba,DistributionMap(lev),1,0); ComputeDivE(dive, 0, @@ -944,12 +977,21 @@ WarpX::WritePlotFile () const if (plot_finepatch) { - PackPlotDataPtrs(srcmf, Efield_fp[lev]); - amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); + if (do_nodal) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *Efield_fp[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_fp[lev][1], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *Efield_fp[lev][2], 0, 1); + } + else + { + PackPlotDataPtrs(srcmf, Efield_fp[lev]); + amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_fp[lev][1], 0, 1); + MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Efield_fp[lev][1], 0, 1); #endif + } if (lev == 0) { varnames.push_back("Ex_fp"); @@ -958,12 +1000,21 @@ WarpX::WritePlotFile () const } dcomp += 3; - PackPlotDataPtrs(srcmf, Bfield_fp[lev]); - amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); + if (do_nodal) + { + amrex::average_node_to_cellcenter(*mf[lev], dcomp , *Bfield_fp[lev][0], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+1, *Bfield_fp[lev][1], 0, 1); + amrex::average_node_to_cellcenter(*mf[lev], dcomp+2, *Bfield_fp[lev][2], 0, 1); + } + else + { + PackPlotDataPtrs(srcmf, Bfield_fp[lev]); + amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); - MultiFab::Copy(*mf[lev], *Bfield_fp[lev][1], 0, dcomp+1, 1, ngrow); + MultiFab::Copy(*mf[lev], *mf[lev], dcomp+1, dcomp+2, 1, ngrow); + MultiFab::Copy(*mf[lev], *Bfield_fp[lev][1], 0, dcomp+1, 1, ngrow); #endif + } if (lev == 0) { varnames.push_back("Bx_fp"); @@ -982,6 +1033,7 @@ WarpX::WritePlotFile () const } else { + if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch"); std::array<std::unique_ptr<MultiFab>, 3> E = getInterpolatedE(lev); PackPlotDataPtrs(srcmf, E); amrex::average_edge_to_cellcenter(*mf[lev], dcomp, srcmf); @@ -1005,6 +1057,7 @@ WarpX::WritePlotFile () const } else { + if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch"); std::array<std::unique_ptr<MultiFab>, 3> B = getInterpolatedB(lev); PackPlotDataPtrs(srcmf, B); amrex::average_face_to_cellcenter(*mf[lev], dcomp, srcmf); @@ -1091,46 +1144,46 @@ WarpX::WritePlotFile () const // Plot fine patch if (plot_finepatch) { - if (plot_raw_fields_guards) { - VisMF::Write(*Efield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); - VisMF::Write(*Efield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); - VisMF::Write(*Efield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); - VisMF::Write(*Bfield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); - VisMF::Write(*Bfield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); - VisMF::Write(*Bfield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); - VisMF::Write(*current_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); - VisMF::Write(*current_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); - VisMF::Write(*current_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); - } else { - const DistributionMapping& dm = DistributionMap(lev); - MultiFab Ex(Efield_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab Ey(Efield_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab Ez(Efield_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab Bx(Bfield_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab By(Bfield_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab Bz(Bfield_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab jx(current_fp[lev][0]->boxArray(), dm, 1, 0); - MultiFab jy(current_fp[lev][1]->boxArray(), dm, 1, 0); - MultiFab jz(current_fp[lev][2]->boxArray(), dm, 1, 0); - MultiFab::Copy(Ex, *Efield_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(Ey, *Efield_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Ez, *Efield_fp[lev][2], 0, 0, 1, 0); - MultiFab::Copy(Bx, *Bfield_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(By, *Bfield_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(Bz, *Bfield_fp[lev][2], 0, 0, 1, 0); - MultiFab::Copy(jx, *current_fp[lev][0], 0, 0, 1, 0); - MultiFab::Copy(jy, *current_fp[lev][1], 0, 0, 1, 0); - MultiFab::Copy(jz, *current_fp[lev][2], 0, 0, 1, 0); - VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); - VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); - VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); - VisMF::Write(Bx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); - VisMF::Write(By, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); - VisMF::Write(Bz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); - VisMF::Write(jx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); - VisMF::Write(jy, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); - VisMF::Write(jz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); - } + if (plot_raw_fields_guards) { + VisMF::Write(*Efield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); + VisMF::Write(*Efield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); + VisMF::Write(*Efield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); + VisMF::Write(*Bfield_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); + VisMF::Write(*Bfield_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); + VisMF::Write(*Bfield_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); + VisMF::Write(*current_fp[lev][0], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); + VisMF::Write(*current_fp[lev][1], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); + VisMF::Write(*current_fp[lev][2], amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); + } else { + const DistributionMapping& dm = DistributionMap(lev); + MultiFab Ex(Efield_fp[lev][0]->boxArray(), dm, 1, 0); + MultiFab Ey(Efield_fp[lev][1]->boxArray(), dm, 1, 0); + MultiFab Ez(Efield_fp[lev][2]->boxArray(), dm, 1, 0); + MultiFab Bx(Bfield_fp[lev][0]->boxArray(), dm, 1, 0); + MultiFab By(Bfield_fp[lev][1]->boxArray(), dm, 1, 0); + MultiFab Bz(Bfield_fp[lev][2]->boxArray(), dm, 1, 0); + MultiFab jx(current_fp[lev][0]->boxArray(), dm, 1, 0); + MultiFab jy(current_fp[lev][1]->boxArray(), dm, 1, 0); + MultiFab jz(current_fp[lev][2]->boxArray(), dm, 1, 0); + MultiFab::Copy(Ex, *Efield_fp[lev][0], 0, 0, 1, 0); + MultiFab::Copy(Ey, *Efield_fp[lev][1], 0, 0, 1, 0); + MultiFab::Copy(Ez, *Efield_fp[lev][2], 0, 0, 1, 0); + MultiFab::Copy(Bx, *Bfield_fp[lev][0], 0, 0, 1, 0); + MultiFab::Copy(By, *Bfield_fp[lev][1], 0, 0, 1, 0); + MultiFab::Copy(Bz, *Bfield_fp[lev][2], 0, 0, 1, 0); + MultiFab::Copy(jx, *current_fp[lev][0], 0, 0, 1, 0); + MultiFab::Copy(jy, *current_fp[lev][1], 0, 0, 1, 0); + MultiFab::Copy(jz, *current_fp[lev][2], 0, 0, 1, 0); + VisMF::Write(Ex, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ex_fp")); + VisMF::Write(Ey, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ey_fp")); + VisMF::Write(Ez, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Ez_fp")); + VisMF::Write(Bx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bx_fp")); + VisMF::Write(By, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "By_fp")); + VisMF::Write(Bz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "Bz_fp")); + VisMF::Write(jx, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jx_fp")); + VisMF::Write(jy, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jy_fp")); + VisMF::Write(jz, amrex::MultiFabFileFullPrefix(lev, raw_plotfilename, level_prefix, "jz_fp")); + } } // Plot coarse patch diff --git a/Source/WarpXParticleContainer.cpp b/Source/WarpXParticleContainer.cpp index c323679a1..2bed6e0d3 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/WarpXParticleContainer.cpp @@ -243,7 +243,6 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, int thread_num, int lev, Real dt ) { Real *jx_ptr, *jy_ptr, *jz_ptr; - const int *jxntot, *jyntot, *jzntot; const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); @@ -293,15 +292,15 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, local_jz_ptr->setVal(0.0, b, 0, 1); }); - jxntot = local_jx[thread_num]->length(); - jyntot = local_jy[thread_num]->length(); - jzntot = local_jz[thread_num]->length(); + auto jxntot = local_jx[thread_num]->length(); + auto jyntot = local_jy[thread_num]->length(); + auto jzntot = local_jz[thread_num]->length(); BL_PROFILE_VAR_START(blp_pxr_cd); warpx_current_deposition( - jx_ptr, &ngJ, jxntot, - jy_ptr, &ngJ, jyntot, - jz_ptr, &ngJ, jzntot, + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), &np_current, m_xp[thread_num].dataPtr(), m_yp[thread_num].dataPtr(), @@ -380,16 +379,16 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, { local_jz_ptr->setVal(0.0, b, 0, 1); }); - jxntot = local_jx[thread_num]->length(); - jyntot = local_jy[thread_num]->length(); - jzntot = local_jz[thread_num]->length(); + auto jxntot = local_jx[thread_num]->length(); + auto jyntot = local_jy[thread_num]->length(); + auto jzntot = local_jz[thread_num]->length(); long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); warpx_current_deposition( - jx_ptr, &ngJ, jxntot, - jy_ptr, &ngJ, jyntot, - jz_ptr, &ngJ, jzntot, + jx_ptr, &ngJ, jxntot.getVect(), + jy_ptr, &ngJ, jyntot.getVect(), + jz_ptr, &ngJ, jzntot.getVect(), &ncrse, m_xp[thread_num].dataPtr() +np_current, m_yp[thread_num].dataPtr() +np_current, @@ -449,7 +448,6 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, long ngRho = rhomf->nGrow(); Real* data_ptr; Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); - const int *rholen; const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); @@ -467,7 +465,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, }); data_ptr = local_rho[thread_num]->dataPtr(); - rholen = local_rho[thread_num]->length(); + auto rholen = local_rho[thread_num]->length(); #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -519,7 +517,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, }); data_ptr = local_rho[thread_num]->dataPtr(); - rholen = local_rho[thread_num]->length(); + auto rholen = local_rho[thread_num]->length(); #if (AMREX_SPACEDIM == 3) const long nx = rholen[0]-1-2*ngRho; const long ny = rholen[1]-1-2*ngRho; @@ -662,7 +660,6 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) // Data on the grid Real* data_ptr; - const int *rholen; FArrayBox& rhofab = (*rho)[pti]; #ifdef _OPENMP Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); @@ -671,11 +668,11 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) local_rho.resize(tile_box); local_rho = 0.0; data_ptr = local_rho.dataPtr(); - rholen = local_rho.length(); + auto rholen = local_rho.length(); #else const std::array<Real, 3>& xyzmin = xyzmin_grid; data_ptr = rhofab.dataPtr(); - rholen = rhofab.length(); + auto rholen = rhofab.length(); #endif #if (AMREX_SPACEDIM == 3) diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H index 578785a14..db7a5f9bd 100644 --- a/Source/WarpX_f.H +++ b/Source/WarpX_f.H @@ -137,7 +137,7 @@ extern "C" const amrex::Real* byg, const int* byg_lo, const int* byg_hi, const amrex::Real* bzg, const int* bzg_lo, const int* bzg_hi, const int* ll4symtry, const int* l_lower_order_in_v, - const long* lvect, + const int* l_nodal, const long* lvect, const long* field_gathe_algo); // Particle pusher (velocity and position) diff --git a/Source/WarpX_picsar.F90 b/Source/WarpX_picsar.F90 index 13168ceca..8e9d2db3a 100644 --- a/Source/WarpX_picsar.F90 +++ b/Source/WarpX_picsar.F90 @@ -75,7 +75,7 @@ contains ex,ey,ez,bx,by,bz,ixyzmin,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & exg,exg_lo,exg_hi,eyg,eyg_lo,eyg_hi,ezg,ezg_lo,ezg_hi, & bxg,bxg_lo,bxg_hi,byg,byg_lo,byg_hi,bzg,bzg_lo,bzg_hi, & - ll4symtry,l_lower_order_in_v, & + ll4symtry,l_lower_order_in_v, l_nodal,& lvect,field_gathe_algo) & bind(C, name="warpx_geteb_energy_conserving") @@ -87,12 +87,12 @@ contains real(amrex_real), intent(in) :: xmin,ymin,zmin,dx,dy,dz integer(c_long), intent(in) :: field_gathe_algo integer(c_long), intent(in) :: np,nox,noy,noz - integer(c_int), intent(in) :: ll4symtry,l_lower_order_in_v + integer(c_int), intent(in) :: ll4symtry,l_lower_order_in_v, l_nodal integer(c_long),intent(in) :: lvect real(amrex_real), intent(in), dimension(np) :: xp,yp,zp real(amrex_real), intent(out), dimension(np) :: ex,ey,ez,bx,by,bz real(amrex_real),intent(in):: exg(*), eyg(*), ezg(*), bxg(*), byg(*), bzg(*) - logical(pxr_logical) :: pxr_ll4symtry, pxr_l_lower_order_in_v + logical(pxr_logical) :: pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal ! Compute the number of valid cells and guard cells integer(c_long) :: exg_nvalid(AMREX_SPACEDIM), eyg_nvalid(AMREX_SPACEDIM), ezg_nvalid(AMREX_SPACEDIM), & @@ -102,7 +102,8 @@ contains pxr_ll4symtry = ll4symtry .eq. 1 pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 - + pxr_l_nodal = l_nodal .eq. 1 + exg_nguards = ixyzmin - exg_lo eyg_nguards = ixyzmin - eyg_lo ezg_nguards = ixyzmin - ezg_lo @@ -124,7 +125,7 @@ contains bxg,bxg_nguards,bxg_nvalid,& byg,byg_nguards,byg_nvalid,& bzg,bzg_nguards,bzg_nvalid,& - pxr_ll4symtry, pxr_l_lower_order_in_v, & + pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal, & lvect, field_gathe_algo ) end subroutine warpx_geteb_energy_conserving |