diff options
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 14 | ||||
| -rw-r--r-- | Source/LaserParticleContainer.cpp | 5 | ||||
| -rw-r--r-- | Source/Make.package | 7 | ||||
| -rw-r--r-- | Source/WarpX.H | 18 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 68 | ||||
| -rw-r--r-- | Source/WarpXBoostedFrameDiagnostic.H | 84 | ||||
| -rw-r--r-- | Source/WarpXBoostedFrameDiagnostic.cpp | 166 | ||||
| -rw-r--r-- | Source/WarpXEvolve.cpp | 8 | ||||
| -rw-r--r-- | Source/WarpXIO.cpp | 46 | ||||
| -rw-r--r-- | Source/WarpXInitData.cpp | 17 | ||||
| -rw-r--r-- | Source/WarpXParticleContainer.cpp | 16 | ||||
| -rw-r--r-- | Source/WarpX_boosted_frame.F90 | 60 | ||||
| -rw-r--r-- | Source/WarpX_f.H | 6 | ||||
| -rw-r--r-- | Tools/read_raw_data.py | 76 |
14 files changed, 549 insertions, 42 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 30359d8a7..7cc4c53c1 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -260,3 +260,17 @@ Diagnostics and output * ``amr.plot_int`` (`integer`) The number of PIC cycles inbetween two consecutive data dumps. Use a negative number to disable data dumping. + +* ``warpx.do_boosted_frame_diagnostic`` (`0 or 1`) + Whether to use the **back-transformed diagnostics** (i.e. diagnostics that + perform on-the-fly conversion to the laboratory frame, when running + boosted-frame simulations) + +* ``warpx.num_snapshots_lab`` (`integer`) + Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1``. + The number of lab-frame snapshots that will be written. + +* ``warpx.dt_snapshots_lab`` (`float`, in seconds) + Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1``. + The time interval inbetween the lab-frame snapshots (where this + time interval is expressed in the laboratory frame). diff --git a/Source/LaserParticleContainer.cpp b/Source/LaserParticleContainer.cpp index b9fbd8d7d..b7d49e402 100644 --- a/Source/LaserParticleContainer.cpp +++ b/Source/LaserParticleContainer.cpp @@ -87,9 +87,8 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies) p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s }; Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); - if (std::abs(dp) > 1.e-14) { - amrex::Abort("Laser plane vector is not perpendicular to the main polarization vector"); - } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, + "Laser plane vector is not perpendicular to the main polarization vector"); p_Y = CrossProduct(nvec, p_X); // The second polarization vector diff --git a/Source/Make.package b/Source/Make.package index c9c8973f4..8255eae63 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -5,14 +5,14 @@ endif CEXE_sources += WarpXWrappers.cpp WarpX_py.cpp CEXE_sources += WarpX.cpp WarpXInitData.cpp WarpXEvolve.cpp WarpXIO.cpp WarpXProb.cpp WarpXRegrid.cpp -CEXE_sources += WarpXTagging.cpp WarpXComm.cpp WarpXMove.cpp +CEXE_sources += WarpXTagging.cpp WarpXComm.cpp WarpXMove.cpp WarpXBoostedFrameDiagnostic.cpp CEXE_sources += ParticleIO.cpp CEXE_sources += ParticleContainer.cpp WarpXParticleContainer.cpp PhysicalParticleContainer.cpp LaserParticleContainer.cpp CEXE_headers += WarpXWrappers.h WarpX_py.H -CEXE_headers += WarpX.H WarpX_f.H WarpXConst.H +CEXE_headers += WarpX.H WarpX_f.H WarpXConst.H WarpXBoostedFrameDiagnostic.H CEXE_headers += ParticleContainer.H WarpXParticleContainer.H PhysicalParticleContainer.H LaserParticleContainer.H @@ -22,7 +22,8 @@ CEXE_sources += PlasmaInjector.cpp CustomDensityProb.cpp CEXE_sources += WarpXPML.cpp CEXE_headers += WarpXPML.H -F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90 WarpX_pml.F90 WarpX_electrostatic.F90 WarpX_filter.F90 +F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90 WarpX_pml.F90 WarpX_electrostatic.F90 +F90EXE_sources += WarpX_boosted_frame.F90 WarpX_filter.F90 ifeq ($(USE_OPENBC_POISSON),TRUE) F90EXE_sources += openbc_f.F90 diff --git a/Source/WarpX.H b/Source/WarpX.H index 6c2b4f402..363d1c4a8 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -16,6 +16,7 @@ #include <ParticleContainer.H> #include <WarpXPML.H> +#include <WarpXBoostedFrameDiagnostic.H> class NoOpPhysBC : public amrex::PhysBCFunctBase @@ -77,6 +78,11 @@ public: static bool use_filter; static bool serialize_ics; + // Back transformation diagnostic + static bool do_boosted_frame_diagnostic; + static int num_snapshots_lab; + static Real dt_snapshots_lab; + // Boosted frame parameters static amrex::Real gamma_boost; static amrex::Real beta_boost; @@ -229,10 +235,15 @@ private: void InitPML (); void ComputePMLFactors (); + void InitDiagnostics (); + void InjectPlasma(int num_shift, int dir); void WriteWarpXHeader(const std::string& name) const; void WriteJobInfo (const std::string& dir) const; + + std::unique_ptr<amrex::MultiFab> GetCellCenteredData(); + static void PackPlotDataPtrs (amrex::Vector<const amrex::MultiFab*>& pmf, const std::array<std::unique_ptr<amrex::MultiFab>,3>& data); @@ -288,6 +299,9 @@ private: // Particle container std::unique_ptr<MultiParticleContainer> mypc; + // Boosted Frame Diagnostics + std::unique_ptr<BoostedFrameDiagnostic> myBFD; + // // Fields: First array for level, second for direction // @@ -323,8 +337,8 @@ private: // Moving window parameters int do_moving_window = 0; int moving_window_dir = -1; - amrex::Real moving_window_x = std::numeric_limits<amrex::Real>::max(); - amrex::Real moving_window_v = std::numeric_limits<amrex::Real>::max(); + amrex::Real moving_window_x = std::numeric_limits<amrex::Real>::max(); + amrex::Real moving_window_v = std::numeric_limits<amrex::Real>::max(); // Plasma injection parameters int do_plasma_injection = 0; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 622b8b0e6..4ea87d6e5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -41,6 +41,10 @@ bool WarpX::use_laser = false; bool WarpX::use_filter = false; bool WarpX::serialize_ics = false; +bool WarpX::do_boosted_frame_diagnostic = false; +int WarpX::num_snapshots_lab = std::numeric_limits<int>::lowest(); +Real WarpX::dt_snapshots_lab = std::numeric_limits<Real>::lowest(); + #if (BL_SPACEDIM == 3) IntVect WarpX::Bx_nodal_flag(1,0,0); IntVect WarpX::By_nodal_flag(0,1,0); @@ -173,21 +177,21 @@ WarpX::ReadParameters () pp.query("verbose", verbose); pp.query("regrid_int", regrid_int); - // Boosted-frame parameters - pp.query("gamma_boost", gamma_boost); - beta_boost = std::sqrt(1.-1./pow(gamma_boost,2)); - if( gamma_boost > 1. ){ - // Read and normalize the boost direction - pp.queryarr("boost_direction", boost_direction); - Real s = 1.0/std::sqrt(boost_direction[0]*boost_direction[0] + - boost_direction[1]*boost_direction[1] + - boost_direction[2]*boost_direction[2]); + // Boosted-frame parameters + pp.query("gamma_boost", gamma_boost); + beta_boost = std::sqrt(1.-1./pow(gamma_boost,2)); + if( gamma_boost > 1. ){ + // Read and normalize the boost direction + pp.queryarr("boost_direction", boost_direction); + Real s = 1.0/std::sqrt(boost_direction[0]*boost_direction[0] + + boost_direction[1]*boost_direction[1] + + boost_direction[2]*boost_direction[2]); boost_direction = { boost_direction[0]*s, - boost_direction[1]*s, - boost_direction[2]*s }; - } - - pp.queryarr("B_external", B_external); + boost_direction[1]*s, + boost_direction[2]*s }; + } + + pp.queryarr("B_external", B_external); pp.query("do_moving_window", do_moving_window); if (do_moving_window) @@ -224,6 +228,31 @@ WarpX::ReadParameters () 0, num_injected_species); } + pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); + if (do_boosted_frame_diagnostic) { + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, + "gamma_boost must be > 1 to use the boosted frame diagnostic."); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (std::abs(boost_direction[0] - 0.0) < 1.0e-12) and + (std::abs(boost_direction[1] - 0.0) < 1.0e-12) and + (std::abs(boost_direction[2] - 1.0) < 1.0e012) , + "The boosted frame diagnostic currently only works if the boost is in the z direction."); + + pp.get("num_snapshots_lab", num_snapshots_lab); + pp.get("dt_snapshots_lab", dt_snapshots_lab); + pp.get("gamma_boost", gamma_boost); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, + "The moving window should be on if using the boosted frame diagnostic."); + + std::string s; + pp.get("moving_window_dir", s); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), + "The boosted frame diagnostic currently only works if the boost is in the z direction."); + } + pp.query("do_electrostatic", do_electrostatic); pp.query("n_buffer", n_buffer); pp.query("const_dt", const_dt); @@ -271,14 +300,11 @@ WarpX::ReadParameters () pp.query("nox", nox); pp.query("noy", noy); pp.query("noz", noz); - if (nox != noy || nox != noz) { - amrex::Abort("warpx.nox, noy and noz must be equal"); - } - if (nox < 1) { - amrex::Abort("warpx.nox must >= 1"); - } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox == noy and nox == noz , + "warpx.nox, noy and noz must be equal"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox >= 1, "warpx.nox must >= 1"); } - + { ParmParse pp("algo"); pp.query("current_deposition", current_deposition_algo); diff --git a/Source/WarpXBoostedFrameDiagnostic.H b/Source/WarpXBoostedFrameDiagnostic.H new file mode 100644 index 000000000..1decde81f --- /dev/null +++ b/Source/WarpXBoostedFrameDiagnostic.H @@ -0,0 +1,84 @@ +#ifndef WARPX_BoostedFrameDiagnostic_H_ +#define WARPX_BoostedFrameDiagnostic_H_ + +#include <vector> + +#include <AMReX_VisMF.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_ParallelDescriptor.H> + +#include "WarpXConst.H" + +/// +/// BoostedFrameDiagnostic is for handling IO when running in a boosted +/// frame of reference. Because of the relativity of simultaneity, events that +/// are synchronized in the simulation frame are not synchronized in the +/// lab frame. Thus, at a given t_boost, we must write slices of data to +/// multiple output files, each one corresponding to a given time in the lab frame. +/// +class BoostedFrameDiagnostic { + + /// + /// LabSnapShot stores metadata corresponding to a single time + /// snapshot in the lab frame. The snapshot is written to disk + /// in the directory "file_name". zmin_lab, zmax_lab, and t_lab + /// are all constant for a given snapshot. current_z_lab and + /// current_z_boost for each snapshot are updated as the + /// simulation time in the boosted frame advances. + /// + struct LabSnapShot { + + std::string file_name; + amrex::Real t_lab; + amrex::Real zmin_lab; + amrex::Real zmax_lab; + amrex::Real current_z_lab; + amrex::Real current_z_boost; + int file_num; + + LabSnapShot(amrex::Real t_lab_in, amrex::Real zmin_lab_in, + amrex::Real zmax_lab_in, int file_num_in); + + /// + /// This snapshot is at time t_lab, and the simulation is at time t_boost. + /// The Lorentz transformation picks out one slice corresponding to both + /// of those times, at position current_z_boost and current_z_lab in the + /// boosted and lab frames, respectively. + /// + void updateCurrentZPositions(amrex::Real t_boost, amrex::Real inv_gamma, + amrex::Real inv_beta); + + /// + /// Write some useful metadata about this snapshot. + /// + void writeSnapShotHeader(); + }; + + amrex::Real gamma_boost_; + amrex::Real inv_gamma_boost_; + amrex::Real beta_boost_; + amrex::Real inv_beta_boost_; + amrex::Real dz_lab_; + amrex::Real inv_dz_lab_; + amrex::Real dt_snapshots_lab_; + amrex::Real dt_boost_; + int N_snapshots_; + int Nz_lab_; + int boost_direction_; + + std::vector<LabSnapShot> snapshots_; + +public: + + BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, + amrex::Real v_window_lab, amrex::Real dt_snapshots_lab, + int N_snapshots, amrex::Real gamma_boost, + amrex::Real dt_boost, int boost_direction); + + void writeLabFrameData(const amrex::MultiFab& cell_centered_data, + const amrex::Geometry& geom, amrex::Real t_boost); + + void writeMetaData(); +}; + +#endif diff --git a/Source/WarpXBoostedFrameDiagnostic.cpp b/Source/WarpXBoostedFrameDiagnostic.cpp new file mode 100644 index 000000000..f900cee8a --- /dev/null +++ b/Source/WarpXBoostedFrameDiagnostic.cpp @@ -0,0 +1,166 @@ +#include "WarpXBoostedFrameDiagnostic.H" +#include <AMReX_MultiFabUtil.H> +#include "WarpX_f.H" + +using namespace amrex; + +BoostedFrameDiagnostic:: +BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, + Real dt_snapshots_lab, int N_snapshots, + Real gamma_boost, Real dt_boost, + int boost_direction) + : gamma_boost_(gamma_boost), + dt_snapshots_lab_(dt_snapshots_lab), + dt_boost_(dt_boost), + N_snapshots_(N_snapshots), + boost_direction_(boost_direction) +{ + inv_gamma_boost_ = 1.0 / gamma_boost_; + beta_boost_ = std::sqrt(1.0 - inv_gamma_boost_*inv_gamma_boost_); + inv_beta_boost_ = 1.0 / beta_boost_; + + dz_lab_ = PhysConst::c * dt_boost_ * inv_beta_boost_ * inv_gamma_boost_; + inv_dz_lab_ = 1.0 / dz_lab_; + Nz_lab_ = static_cast<int>((zmax_lab - zmin_lab) * inv_dz_lab_); + + writeMetaData(); + + for (int i = 0; i < N_snapshots; ++i) { + Real t_lab = i * dt_snapshots_lab_; + LabSnapShot snapshot(t_lab, zmin_lab + v_window_lab * t_lab, + zmax_lab + v_window_lab * t_lab, i); + snapshots_.push_back(snapshot); + } +} + +void +BoostedFrameDiagnostic:: +writeLabFrameData(const MultiFab& cell_centered_data, const Geometry& geom, Real t_boost) { + + const RealBox& domain_z_boost = geom.ProbDomain(); + const Real zlo_boost = domain_z_boost.lo(boost_direction_); + const Real zhi_boost = domain_z_boost.hi(boost_direction_); + + for (int i = 0; i < N_snapshots_; ++i) { + snapshots_[i].updateCurrentZPositions(t_boost, + inv_gamma_boost_, + inv_beta_boost_); + + if ( (snapshots_[i].current_z_boost < zlo_boost) or + (snapshots_[i].current_z_boost > zhi_boost) or + (snapshots_[i].current_z_lab < snapshots_[i].zmin_lab) or + (snapshots_[i].current_z_lab > snapshots_[i].zmax_lab) ) continue; + + // for each z position, fill a slice with the data. + int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; + + const int ncomp = cell_centered_data.nComp(); + const int start_comp = 0; + std::unique_ptr<MultiFab> slice = amrex::get_slice_data(boost_direction_, + snapshots_[i].current_z_boost, + cell_centered_data, geom, + start_comp, ncomp); + + // transform it to the lab frame + for (MFIter mfi(*slice); mfi.isValid(); ++mfi) { + const Box& box = mfi.validbox(); + const Box& tile_box = mfi.tilebox(); + WRPX_LORENTZ_TRANSFORM_Z(BL_TO_FORTRAN_ANYD((*slice)[mfi]), + BL_TO_FORTRAN_BOX(tile_box), + &gamma_boost_, &beta_boost_); + } + + // and write it to disk. + std::stringstream ss; + ss << snapshots_[i].file_name << "/Level_0/" << Concatenate("slice", i_lab, 5); + VisMF::Write(*slice, ss.str()); + } +} + +void +BoostedFrameDiagnostic:: +writeMetaData() +{ + if (ParallelDescriptor::IOProcessor()) { + std::string DiagnosticDirectory = "lab_frame_data"; + + if (!UtilCreateDirectory(DiagnosticDirectory, 0755)) + CreateDirectoryFailed(DiagnosticDirectory); + + std::string HeaderFileName(DiagnosticDirectory + "/Header"); + std::ofstream HeaderFile(HeaderFileName.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + if(!HeaderFile.good()) + FileOpenFailed(HeaderFileName); + + HeaderFile.precision(17); + + VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); + HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); + + HeaderFile << N_snapshots_ << "\n"; + HeaderFile << dt_snapshots_lab_ << "\n"; + HeaderFile << gamma_boost_ << "\n"; + HeaderFile << beta_boost_ << "\n"; + HeaderFile << dz_lab_ << "\n"; + HeaderFile << Nz_lab_ << "\n"; + } +} + +BoostedFrameDiagnostic::LabSnapShot:: +LabSnapShot(Real t_lab_in, Real zmin_lab_in, + Real zmax_lab_in, int file_num_in) + : t_lab(t_lab_in), + zmin_lab(zmin_lab_in), + zmax_lab(zmax_lab_in), + file_num(file_num_in) +{ + current_z_lab = 0.0; + current_z_boost = 0.0; + file_name = Concatenate("lab_frame_data/snapshot", file_num, 5); + + const int nlevels = 1; + const std::string level_prefix = "Level_"; + + if (!UtilCreateDirectory(file_name, 0755)) + CreateDirectoryFailed(file_name); + for(int i(0); i < nlevels; ++i) { + const std::string &fullpath = LevelFullPath(i, file_name); + if (!UtilCreateDirectory(fullpath, 0755)) + CreateDirectoryFailed(fullpath); + } + + ParallelDescriptor::Barrier(); + + writeSnapShotHeader(); +} + +void +BoostedFrameDiagnostic::LabSnapShot:: +updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta) { + current_z_boost = (t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta; + current_z_lab = (t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta; +} + +void +BoostedFrameDiagnostic::LabSnapShot:: +writeSnapShotHeader() { + if (ParallelDescriptor::IOProcessor()) { + std::string HeaderFileName(file_name + "/Header"); + std::ofstream HeaderFile(HeaderFileName.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + if(!HeaderFile.good()) + FileOpenFailed(HeaderFileName); + + HeaderFile.precision(17); + + VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); + HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); + + HeaderFile << t_lab << "\n"; + HeaderFile << zmin_lab << "\n"; + HeaderFile << zmax_lab << "\n"; + } +} diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 86b7bc2ce..df6822962 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -160,8 +160,7 @@ WarpX::EvolveEM (int numsteps) numsteps_max = std::min(istep[0]+numsteps, max_step); } - bool max_time_reached = false; - + bool max_time_reached = false; for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step) { if (warpx_py_print_step) { @@ -255,6 +254,11 @@ WarpX::EvolveEM (int numsteps) t_new[i] = cur_time; } + if (do_boosted_frame_diagnostic) { + std::unique_ptr<MultiFab> cell_centered_data = GetCellCenteredData(); + myBFD->writeLabFrameData(*cell_centered_data, geom[0], cur_time); + } + if (to_make_plot) { FillBoundaryE(); diff --git a/Source/WarpXIO.cpp b/Source/WarpXIO.cpp index adb7b6cf7..0a9efa486 100644 --- a/Source/WarpXIO.cpp +++ b/Source/WarpXIO.cpp @@ -393,6 +393,52 @@ WarpX::InitFromCheckpoint () } +std::unique_ptr<MultiFab> +WarpX::GetCellCenteredData() { + + const int ng = 1; + const int nc = 10; + const int lev = 0; + auto cc = std::unique_ptr<MultiFab>( new MultiFab(boxArray(lev), + DistributionMap(lev), + nc, ng) ); + + Array<const MultiFab*> srcmf(BL_SPACEDIM); + int dcomp = 0; + + // first the electric field + PackPlotDataPtrs(srcmf, Efield_aux[lev]); + amrex::average_edge_to_cellcenter(*cc, dcomp, srcmf); +#if (BL_SPACEDIM == 2) + MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); + amrex::average_node_to_cellcenter(*cc, dcomp+1, *Efield_aux[lev][1], 0, 1); +#endif + dcomp += 3; + + // then the magnetic field + PackPlotDataPtrs(srcmf, Bfield_aux[lev]); + amrex::average_face_to_cellcenter(*cc, dcomp, srcmf); +#if (BL_SPACEDIM == 2) + MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); + MultiFab::Copy(*cc, *Bfield_aux[lev][1], 0, dcomp+1, 1, ng); +#endif + dcomp += 3; + + // then the current density + PackPlotDataPtrs(srcmf, current_fp[lev]); + amrex::average_edge_to_cellcenter(*cc, dcomp, srcmf); +#if (BL_SPACEDIM == 2) + MultiFab::Copy(*cc, *cc, dcomp+1, dcomp+2, 1, ng); + amrex::average_node_to_cellcenter(*cc, dcomp+1, *current_fp[lev][1], 0, 1); +#endif + dcomp += 3; + + const std::unique_ptr<MultiFab>& charge_density = mypc->GetChargeDensity(lev); + amrex::average_node_to_cellcenter(*cc, dcomp, *charge_density, 0, 1); + + return cc; +} + void WarpX::WritePlotFile () const { diff --git a/Source/WarpXInitData.cpp b/Source/WarpXInitData.cpp index a80d3e9dd..873adcbaf 100644 --- a/Source/WarpXInitData.cpp +++ b/Source/WarpXInitData.cpp @@ -29,6 +29,8 @@ WarpX::InitData () ComputePMLFactors(); + InitDiagnostics(); + if (ParallelDescriptor::IOProcessor()) { std::cout << "\nGrids Summary:\n"; printGridSummary(std::cout, 0, finestLevel()); @@ -46,6 +48,21 @@ WarpX::InitData () } void +WarpX::InitDiagnostics () { + if (do_boosted_frame_diagnostic) { + const Real* current_lo = geom[0].ProbLo(); + const Real* current_hi = geom[0].ProbHi(); + Real dt_boost = dt[0]; + + myBFD.reset(new BoostedFrameDiagnostic(current_lo[moving_window_dir], + current_hi[moving_window_dir], + moving_window_v, dt_snapshots_lab, + num_snapshots_lab, gamma_boost, dt_boost, + moving_window_dir)); + } +} + +void WarpX::InitFromScratch () { const Real time = 0.0; diff --git a/Source/WarpXParticleContainer.cpp b/Source/WarpXParticleContainer.cpp index 5313536f2..b3cf42ad4 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/WarpXParticleContainer.cpp @@ -281,15 +281,15 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) data_ptr = rhofab.dataPtr(); rholen = rhofab.length(); #endif - + #if (BL_SPACEDIM == 3) - long nx = box.length(0); - long ny = box.length(1); - long nz = box.length(2); -#elif (BL_SPACEDIM == 2) - long nx = box.length(0); - long ny = 0; - long nz = box.length(1); + const long nx = rholen[0]-1-2*ng; + const long ny = rholen[1]-1-2*ng; + const long nz = rholen[2]-1-2*ng; +#else + const long nx = rholen[0]-1-2*ng; + const long ny = 0; + const long nz = rholen[1]-1-2*ng; #endif long nxg = ng; diff --git a/Source/WarpX_boosted_frame.F90 b/Source/WarpX_boosted_frame.F90 new file mode 100644 index 000000000..e7ad70ef8 --- /dev/null +++ b/Source/WarpX_boosted_frame.F90 @@ -0,0 +1,60 @@ + +module warpx_boosted_frame_module + + use iso_c_binding + use amrex_fort_module, only : amrex_real + use constants + + implicit none + +contains + +! +! Given cell-centered data in the boosted reference frame of the simulation, +! this transforms E and B in place so that the multifab now contains values +! in the lab frame. This routine assumes that the simulation frame is moving +! in the positive z direction with respect to the lab frame. +! + subroutine warpx_lorentz_transform_z(data, dlo, dhi, tlo, thi, gamma_boost, beta_boost) & + bind(C, name="warpx_lorentz_transform_z") + + integer(c_int), intent(in) :: dlo(3), dhi(3) + integer(c_int), intent(in) :: tlo(3), thi(3) + real(amrex_real), intent(inout) :: data(dlo(1):dhi(1),dlo(2):dhi(2),dlo(3):dhi(3), 10) + real(amrex_real), intent(in) :: gamma_boost, beta_boost + + integer i, j, k + real(amrex_real) e_lab, b_lab, j_lab, r_lab + + do k = tlo(3), thi(3) + do j = tlo(2), thi(2) + do i = tlo(1), thi(1) + + ! Transform the transverse E and B fields. Note that ez and bz are not + ! changed by the tranform. + e_lab = gamma_boost * (data(i, j, k, 1) + beta_boost*clight*data(i, j, k, 5)) + b_lab = gamma_boost * (data(i, j, k, 5) + beta_boost*data(i, j, k, 1)/clight) + + data(i, j, k, 1) = e_lab + data(i, j, k, 5) = b_lab + + e_lab = gamma_boost * (data(i, j, k, 2) - beta_boost*clight*data(i, j, k, 4)) + b_lab = gamma_boost * (data(i, j, k, 4) - beta_boost*data(i, j, k, 2)/clight) + + data(i, j, k, 2) = e_lab + data(i, j, k, 4) = b_lab + + ! Transform the charge and current density. Only the z component of j is affected. + j_lab = gamma_boost*(data(i, j, k, 9) + beta_boost*clight*data(i, j, k, 10)) + r_lab = gamma_boost*(data(i, j, k, 10) + beta_boost*data(i, j, k, 9)/clight) + + data(i, j, k, 9) = j_lab + data(i, j, k, 10) = r_lab + + end do + end do + end do + + end subroutine warpx_lorentz_transform_z + +end module warpx_boosted_frame_module diff --git a/Source/WarpX_f.H b/Source/WarpX_f.H index 5f922dab7..62b26e4f6 100644 --- a/Source/WarpX_f.H +++ b/Source/WarpX_f.H @@ -26,6 +26,7 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d +#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z #define WRPX_FILTER warpx_filter_3d #elif (BL_SPACEDIM == 2) @@ -53,6 +54,7 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d +#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z #define WRPX_FILTER warpx_filter_2d #endif @@ -256,6 +258,10 @@ extern "C" const amrex::Real* dt, const amrex::Real* prob_lo, const amrex::Real* prob_hi); + void WRPX_LORENTZ_TRANSFORM_Z(amrex::Real* data, const int* dlo, const int* dhi, + const int* tlo, const int* thi, + const amrex::Real* gamma_boost, const amrex::Real* beta_boost); + // These functions are used to evolve E and B in the PML void WRPX_COMPUTE_DIVB (const int* lo, const int* hi, diff --git a/Tools/read_raw_data.py b/Tools/read_raw_data.py index 1ef08dcbf..7dc36520d 100644 --- a/Tools/read_raw_data.py +++ b/Tools/read_raw_data.py @@ -4,6 +4,7 @@ from collections import namedtuple HeaderInfo = namedtuple('HeaderInfo', ['version', 'how', 'ncomp', 'nghost']) +_component_names = ['Ex', 'Ey', 'Ez', 'Bx', 'By', 'Bz', 'jx', 'jy', 'jz', 'rho'] def read_data(plt_file): ''' @@ -43,6 +44,46 @@ def read_data(plt_file): return all_data +def read_lab_snapshot(snapshot): + ''' + + This reads the data from one of the lab frame snapshots generated when + WarpX is run with boosted frame diagnostics turned on. It returns a + dictionary of numpy arrays, where each key corresponds to one of the + data fields ("Ex", "By,", etc... ). These values are cell-centered. + + ''' + + hdrs = glob(snapshot + "/Level_0/slice?????_H") + hdrs.sort() + + boxes, file_names, offsets, header = _read_header(hdrs[0]) + dom_lo, dom_hi = _combine_boxes(boxes) + domain_size = dom_hi - dom_lo + 1 + + space_dim = len(dom_lo) + + data = {} + for i in range(header.ncomp): + if space_dim == 3: + data[component_names[i]] = np.zeros((domain_size[0], domain_size[1], len(hdrs))) + elif space_dim == 2: + data[component_names[i]] = np.zeros((domain_size[0], len(hdrs))) + + for i, hdr in enumerate(hdrs): + slice_data = _read_slice(snapshot, hdr) + if data is None: + data = slice_data + else: + for k,v in slice_data.items(): + if space_dim == 3: + data[k][:,:,i] = v[:,:,0] + elif space_dim == 2: + data[k][:,i] = v[:,0] + + return data + + def _get_field_names(raw_file): header_files = glob(raw_file + "*_H") return [hf.split("/")[-1][:-2] for hf in header_files] @@ -59,8 +100,7 @@ def _line_to_numpy_arrays(line): return lo_corner, hi_corner, node_type -def _read_header(raw_file, field): - header_file = raw_file + field + "_H" +def _read_header(header_file): with open(header_file, "r") as f: version = int(f.readline()) @@ -105,7 +145,8 @@ def _combine_boxes(boxes): def _read_field(raw_file, field_name): - boxes, file_names, offsets, header = _read_header(raw_file, field_name) + header_file = raw_file + field + "_H" + boxes, file_names, offsets, header = _read_header(header_file) ng = header.nghost lo, hi = _combine_boxes(boxes) @@ -125,6 +166,35 @@ def _read_field(raw_file, field_name): return data + +def _read_slice(snapshot, header_fn): + + boxes, file_names, offsets, header = _read_header(header_fn) + + ng = header.nghost + dom_lo, dom_hi = _combine_boxes(boxes) + + all_data = {} + for i in range(header.ncomp): + all_data[_component_names[i]] = np.zeros(dom_hi - dom_lo + 1) + + for box, fn, offset in zip(boxes, file_names, offsets): + lo = box[0] - dom_lo + hi = box[1] - dom_lo + shape = hi - lo + 1 + size = np.product(shape) + with open(snapshot + "/Level_0/" + fn, "rb") as f: + f.seek(offset) + f.readline() # always skip the first line + arr = np.fromfile(f, 'float64', header.ncomp*size) + for i in range(header.ncomp): + comp_data = arr[i*size:(i+1)*size].reshape(shape, order='F') + data = all_data[_component_names[i]] + data[[slice(l,h+1) for l, h in zip(lo, hi)]] = comp_data + all_data[_component_names[i]] = data + return all_data + + if __name__ == "__main__": import matplotlib matplotlib.use('Agg') |