From 021f36f6d208fe77f8f63563aa0548d10737297b Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 19 Nov 2024 18:20:59 -0800 Subject: [PATCH 01/23] Add FieldPoyntingFlux reduced diagnostic --- Python/pywarpx/picmi.py | 1 + .../Diagnostics/ReducedDiags/CMakeLists.txt | 1 + .../ReducedDiags/FieldPoyntingFlux.H | 39 ++++ .../ReducedDiags/FieldPoyntingFlux.cpp | 201 ++++++++++++++++++ Source/Diagnostics/ReducedDiags/Make.package | 1 + .../ReducedDiags/MultiReducedDiags.cpp | 2 + 6 files changed, 245 insertions(+) create mode 100644 Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H create mode 100644 Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index afd28851f70..ac4232b25af 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -3968,6 +3968,7 @@ def __init__( "FieldEnergy", "FieldMomentum", "FieldMaximum", + "FieldPoyntingFlux", "RhoMaximum", "ParticleNumber", "LoadBalanceCosts", diff --git a/Source/Diagnostics/ReducedDiags/CMakeLists.txt b/Source/Diagnostics/ReducedDiags/CMakeLists.txt index bbf1b6b65b0..4fbfc489aba 100644 --- a/Source/Diagnostics/ReducedDiags/CMakeLists.txt +++ b/Source/Diagnostics/ReducedDiags/CMakeLists.txt @@ -9,6 +9,7 @@ foreach(D IN LISTS WarpX_DIMS) FieldEnergy.cpp FieldMaximum.cpp FieldMomentum.cpp + FieldPoyntingFlux.cpp FieldProbe.cpp FieldProbeParticleContainer.cpp FieldReduction.cpp diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H new file mode 100644 index 00000000000..caeefec677e --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -0,0 +1,39 @@ +/* Copyright 2019-2020 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDPOYTINGFLUX_H_ +#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDPOYTINGFLUX_H_ + +#include "ReducedDiags.H" + +#include + +/** + * \brief This class mainly contains a function that computes the field Poynting flux. + */ +class FieldPoyntingFlux : public ReducedDiags +{ +public: + + /** + * \brief Constructor + * + * \param[in] rd_name reduced diags names + */ + FieldPoyntingFlux(const std::string& rd_name); + + /** + * \brief This function computes the electromagnetic Poynting flux, + * obtained by integrating the electromagnetic Poynting flux density g = eps0 * (E x B) + * on the surface of the domain. + * + * \param[in] step current time step + */ + void ComputeDiags(int step) final; +}; + +#endif diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp new file mode 100644 index 00000000000..b08e5689a20 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -0,0 +1,201 @@ +/* Copyright 2019-2020 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "FieldPoyntingFlux.H" + +#include "Fields.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXConst.H" +#include "WarpX.H" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace amrex; +using warpx::fields::FieldType; + +FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) + : ReducedDiags{rd_name} +{ + // RZ coordinate is not working +#if (defined WARPX_DIM_RZ) + WARPX_ABORT_WITH_MESSAGE( + "FieldPoyntingFlux reduced diagnostics not implemented in RZ geometry"); +#endif + + // Resize data array + // lo and hi is 2 + // vector components is 3 + // space dims is AMREX_SPACEDIM + // The order will be (Sx, Sy, Sz) for low faces, then high faces + m_data.resize(2*3*AMREX_SPACEDIM, 0.0_rt); + + if (amrex::ParallelDescriptor::IOProcessor()) + { + if (m_write_header) + { + // Open file + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, std::ofstream::out}; + + int c = 0; + + // Write header row + ofs << "#"; + ofs << "[" << c++ << "]step()"; + ofs << m_sep; + ofs << "[" << c++ << "]time(s)"; + + std::vector sides = {"lo", "hi"}; + std::vector coords = {"x", "y", "z"}; + + // Only on level 0 + for (int iside = 0; iside < 2; iside++) { + for (int ic = 0; ic < AMREX_SPACEDIM; ic++) { + for (int iv = 0; iv < 3; iv++) { + ofs << m_sep; + ofs << "[" << c++ << "]flux_" + coords[iv] + "_" + sides[iside] + "_" + coords[ic] +"(W)"; + }}} + + ofs << "\n"; + ofs.close(); + } + } +} + +void FieldPoyntingFlux::ComputeDiags (int step) +{ + using ablastr::fields::Direction; + + // Check if the diags should be done + if (!m_intervals.contains(step+1)) { return; } + + int const lev = 0; + + // Get a reference to WarpX instance + auto & warpx = WarpX::GetInstance(); + amrex::Box domain_box = warpx.Geom(lev).Domain(); + domain_box.surroundingNodes(); + + // Get MultiFab data at given refinement level + const amrex::MultiFab & Ex = *warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, lev); + const amrex::MultiFab & Ey = *warpx.m_fields.get(FieldType::Efield_aux, Direction{1}, lev); + const amrex::MultiFab & Ez = *warpx.m_fields.get(FieldType::Efield_aux, Direction{2}, lev); + const amrex::MultiFab & Bx = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{0}, lev); + const amrex::MultiFab & By = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{1}, lev); + const amrex::MultiFab & Bz = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{2}, lev); + + // Coarsening ratio (no coarsening) + const amrex::GpuArray cr{1,1,1}; + // Reduction component (fourth component in Array4) + constexpr int comp = 0; + + // Index type (staggering) of each MultiFab + // (with third component set to zero in 2D) + amrex::GpuArray Ex_stag{0,0,0}; + amrex::GpuArray Ey_stag{0,0,0}; + amrex::GpuArray Ez_stag{0,0,0}; + amrex::GpuArray Bx_stag{0,0,0}; + amrex::GpuArray By_stag{0,0,0}; + amrex::GpuArray Bz_stag{0,0,0}; + for (int i = 0; i < AMREX_SPACEDIM; ++i) + { + Ex_stag[i] = Ex.ixType()[i]; + Ey_stag[i] = Ey.ixType()[i]; + Ez_stag[i] = Ez.ixType()[i]; + Bx_stag[i] = Bx.ixType()[i]; + By_stag[i] = By.ixType()[i]; + Bz_stag[i] = Bz.ixType()[i]; + } + + for (amrex::OrientationIter face; face; ++face) { + amrex::Box boundary = amrex::bdryNode(domain_box, face()); + + // Get cell area + const amrex::Real *dx = warpx.Geom(lev).CellSize(); + std::array dxtemp = {AMREX_D_DECL(dx[0], dx[1], dx[2])}; + dxtemp[face().coordDir()] = 1._rt; + const amrex::Real dA = AMREX_D_TERM(dxtemp[0], *dxtemp[1], *dxtemp[2]); + + // Node-centered in the face direction, Cell-centered in other directions + amrex::GpuArray cc{0,0,0}; + cc[face().coordDir()] = 1; + + amrex::ReduceOps reduce_ops; + amrex::ReduceData reduce_data(reduce_ops); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + // Loop over boxes, interpolate E,B data to cell face centers + // and compute sum over cells of (E x B) components + for (amrex::MFIter mfi(Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Array4 & Ex_arr = Ex[mfi].array(); + const amrex::Array4 & Ey_arr = Ey[mfi].array(); + const amrex::Array4 & Ez_arr = Ez[mfi].array(); + const amrex::Array4 & Bx_arr = Bx[mfi].array(); + const amrex::Array4 & By_arr = By[mfi].array(); + const amrex::Array4 & Bz_arr = Bz[mfi].array(); + + amrex::Box box = enclosedCells(mfi.nodaltilebox()); + box.surroundingNodes(face().coordDir()); + + // Find the intersection with the boundary + box &= boundary; + + // Compute E x B + reduce_ops.eval(box, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple + { + const amrex::Real Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); + const amrex::Real Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); + const amrex::Real Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); + + const amrex::Real Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); + const amrex::Real By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); + const amrex::Real Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); + + return {Ey_cc * Bz_cc - Ez_cc * By_cc, + Ez_cc * Bx_cc - Ex_cc * Bz_cc, + Ex_cc * By_cc - Ey_cc * Bx_cc}; + }); + } + + auto r = reduce_data.value(); + int ii = int(face())*3; + m_data[ii+0] = amrex::get<0>(r)/PhysConst::mu0*dA; + m_data[ii+1] = amrex::get<1>(r)/PhysConst::mu0*dA; + m_data[ii+2] = amrex::get<2>(r)/PhysConst::mu0*dA; + + } + + amrex::ParallelDescriptor::ReduceRealSum(m_data.data(), static_cast(m_data.size())); + +} diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package index 2611831a3dd..4d2e4d7def9 100644 --- a/Source/Diagnostics/ReducedDiags/Make.package +++ b/Source/Diagnostics/ReducedDiags/Make.package @@ -7,6 +7,7 @@ CEXE_sources += DifferentialLuminosity.cpp CEXE_sources += FieldEnergy.cpp CEXE_sources += FieldMaximum.cpp CEXE_sources += FieldMomentum.cpp +CEXE_sources += FieldPoyntingFlux.cpp CEXE_sources += FieldProbe.cpp CEXE_sources += FieldProbeParticleContainer.cpp CEXE_sources += FieldReduction.cpp diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 5035eac58a8..5ddc6bc4fe7 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -13,6 +13,7 @@ #include "FieldEnergy.H" #include "FieldMaximum.H" #include "FieldMomentum.H" +#include "FieldPoyntingFlux.H" #include "FieldProbe.H" #include "FieldReduction.H" #include "LoadBalanceCosts.H" @@ -66,6 +67,7 @@ MultiReducedDiags::MultiReducedDiags () {"FieldEnergy", [](CS s){return std::make_unique(s);}}, {"FieldMaximum", [](CS s){return std::make_unique(s);}}, {"FieldMomentum", [](CS s){return std::make_unique(s);}}, + {"FieldPoyntingFlux", [](CS s){return std::make_unique(s);}}, {"FieldProbe", [](CS s){return std::make_unique(s);}}, {"FieldReduction", [](CS s){return std::make_unique(s);}}, {"LoadBalanceCosts", [](CS s){return std::make_unique(s);}}, From 86864c0fff4cb2a0f3a24ec16a77f5846cea6fa7 Mon Sep 17 00:00:00 2001 From: David Grote Date: Wed, 20 Nov 2024 09:09:19 -0800 Subject: [PATCH 02/23] Fix consts --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index b08e5689a20..0e57088d407 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -135,7 +135,7 @@ void FieldPoyntingFlux::ComputeDiags (int step) } for (amrex::OrientationIter face; face; ++face) { - amrex::Box boundary = amrex::bdryNode(domain_box, face()); + amrex::Box const boundary = amrex::bdryNode(domain_box, face()); // Get cell area const amrex::Real *dx = warpx.Geom(lev).CellSize(); @@ -189,7 +189,7 @@ void FieldPoyntingFlux::ComputeDiags (int step) } auto r = reduce_data.value(); - int ii = int(face())*3; + int const ii = int(face())*3; m_data[ii+0] = amrex::get<0>(r)/PhysConst::mu0*dA; m_data[ii+1] = amrex::get<1>(r)/PhysConst::mu0*dA; m_data[ii+2] = amrex::get<2>(r)/PhysConst::mu0*dA; From 9ab05a07c4d3822f47710242abf4df8f6ad68f00 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 26 Nov 2024 09:28:13 -0800 Subject: [PATCH 03/23] Use correct suffices for data file heading --- .../ReducedDiags/FieldPoyntingFlux.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 0e57088d407..73ccd271393 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -73,14 +73,27 @@ FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) ofs << "[" << c++ << "]time(s)"; std::vector sides = {"lo", "hi"}; - std::vector coords = {"x", "y", "z"}; + +#if defined(WARPX_DIM_3D) + std::vector vector_coords = {"x", "y", "z"}; + std::vector space_coords = {"x", "y", "z"}; +#elif defined(WARPX_DIM_XZ) + std::vector vector_coords = {"x", "y", "z"}; + std::vector space_coords = {"x", "z"}; +#elif defined(WARPX_DIM_1D_Z) + std::vector vector_coords = {"x", "y", "z"}; + std::vector space_coords = {"z"}; +#elif defined(WARPX_DIM_RZ) + std::vector vector_coords = {"r", "t", "z"}; + std::vector space_coords = {"r", "z"}; +#endif // Only on level 0 for (int iside = 0; iside < 2; iside++) { for (int ic = 0; ic < AMREX_SPACEDIM; ic++) { for (int iv = 0; iv < 3; iv++) { ofs << m_sep; - ofs << "[" << c++ << "]flux_" + coords[iv] + "_" + sides[iside] + "_" + coords[ic] +"(W)"; + ofs << "[" << c++ << "]flux_" + vector_coords[iv] + "_" + sides[iside] + "_" + space_coords[ic] +"(W)"; }}} ofs << "\n"; From 86daa4a61bf512abce829286d4af3be6861bf5b6 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 26 Nov 2024 09:31:47 -0800 Subject: [PATCH 04/23] Add ComputeDiagsMidStep --- .../ReducedDiags/FieldPoyntingFlux.H | 23 +++++++++++++++++-- .../ReducedDiags/FieldPoyntingFlux.cpp | 19 +++++++++++++++ .../ReducedDiags/MultiReducedDiags.H | 4 ++++ .../ReducedDiags/MultiReducedDiags.cpp | 14 +++++++++++ .../Diagnostics/ReducedDiags/ReducedDiags.H | 7 ++++++ .../Diagnostics/ReducedDiags/ReducedDiags.cpp | 7 ++++++ .../ImplicitSolvers/SemiImplicitEM.cpp | 2 ++ .../StrangImplicitSpectralEM.cpp | 2 ++ .../ImplicitSolvers/ThetaImplicitEM.cpp | 2 ++ 9 files changed, 78 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H index caeefec677e..487ff563848 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -24,7 +24,23 @@ public: * * \param[in] rd_name reduced diags names */ - FieldPoyntingFlux(const std::string& rd_name); + FieldPoyntingFlux (const std::string& rd_name); + + /** + * \brief Call the routine to compute the Poynting flux if needed + * + * \param[in] step current time step + */ + void ComputeDiags (int step) final; + + /** + * \brief Call the routine to compute the Poynting flux at the mid step time level + * + * \param[in] step current time step + */ + void ComputeDiagsMidStep (int step) final; + +private: /** * \brief This function computes the electromagnetic Poynting flux, @@ -33,7 +49,10 @@ public: * * \param[in] step current time step */ - void ComputeDiags(int step) final; + void ComputePoyntingFlux (int step); + + bool use_mid_step_value = false; + }; #endif diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 73ccd271393..4c8abda1a58 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -103,6 +103,25 @@ FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) } void FieldPoyntingFlux::ComputeDiags (int step) +{ + // This will be called at the end of the time step. Only calculate the + // flux if it had not already been calculated mid step. + amrex::Print() << "ComputeDiags\n"; + if (!use_mid_step_value) { + ComputePoyntingFlux(step); + amrex::Print() << " calling ComputePoyntingFlux\n"; + } +} + +void FieldPoyntingFlux::ComputeDiagsMidStep (int step) +{ + // If this is called, always use the value calculated here. + amrex::Print() << "ComputeDiagsMidStep\n"; + use_mid_step_value = true; + ComputePoyntingFlux(step); +} + +void FieldPoyntingFlux::ComputePoyntingFlux (int step) { using ablastr::fields::Direction; diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H index 1a2f51794c6..64157e14791 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H @@ -49,6 +49,10 @@ public: * @param[in] step current iteration time */ void ComputeDiags (int step); + /** Loop over all ReducedDiags and call their ComputeDiagsMidStep + * @param[in] step current iteration time */ + void ComputeDiagsMidStep (int step); + /** Loop over all ReducedDiags and call their WriteToFile * @param[in] step current iteration time */ void WriteToFile (int step); diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 5ddc6bc4fe7..27faaf9f3ca 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -126,6 +126,20 @@ void MultiReducedDiags::ComputeDiags (int step) } // end void MultiReducedDiags::ComputeDiags +// call functions to compute diags at the mid step time level +void MultiReducedDiags::ComputeDiagsMidStep (int step) +{ + WARPX_PROFILE("MultiReducedDiags::ComputeDiagsMidStep()"); + + // loop over all reduced diags + for (int i_rd = 0; i_rd < static_cast(m_rd_names.size()); ++i_rd) + { + m_multi_rd[i_rd] -> ComputeDiagsMidStep(step); + } + // end loop over all reduced diags +} +// end void MultiReducedDiags::ComputeDiagsMidStep + // function to write data void MultiReducedDiags::WriteToFile (int step) { diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.H b/Source/Diagnostics/ReducedDiags/ReducedDiags.H index 2c942e1df6d..a5c787e636a 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.H +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.H @@ -83,6 +83,13 @@ public: */ virtual void ComputeDiags (int step) = 0; + /** + * function to compute diags at the mid step time level + * + * @param[in] step current time step + */ + virtual void ComputeDiagsMidStep (int step); + /** * write to file function * diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp index ec31d9de81c..c48f1511d84 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp @@ -86,6 +86,13 @@ void ReducedDiags::LoadBalance () // load balancing operations } +void ReducedDiags::ComputeDiagsMidStep (int /*step*/) +{ + // Defines an empty function ComputeDiagsMidStep() to be overwritten if needed. + // Function used to calculate the diagnostic at the mid step time leve + // (instead of at the end of the step). +} + void ReducedDiags::BackwardCompatibility () const { const amrex::ParmParse pp_rd_name(m_rd_name); diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp index 117c3baecaa..536e2a1b749 100644 --- a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp @@ -5,6 +5,7 @@ * License: BSD-3-Clause-LBNL */ #include "SemiImplicitEM.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "WarpX.H" using warpx::fields::FieldType; @@ -83,6 +84,7 @@ void SemiImplicitEM::OneStep ( amrex::Real a_time, // Update WarpX owned Efield_fp to t_{n+1/2} m_WarpX->SetElectricFieldAndApplyBCs( m_E ); + m_WarpX->reduced_diags->ComputeDiagsMidStep(a_step); // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate(); diff --git a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp b/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp index 501cbed10eb..d37f58f6ec3 100644 --- a/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp @@ -6,6 +6,7 @@ */ #include "Fields.H" #include "StrangImplicitSpectralEM.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "WarpX.H" using namespace warpx::fields; @@ -84,6 +85,7 @@ void StrangImplicitSpectralEM::OneStep ( amrex::Real a_time, // Update WarpX owned Efield_fp and Bfield_fp to t_{n+1/2} UpdateWarpXFields( m_E, half_time ); + m_WarpX->reduced_diags->ComputeDiagsMidStep(a_step); // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate(); diff --git a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp index 8ca592517ac..b9771393c00 100644 --- a/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp +++ b/Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp @@ -6,6 +6,7 @@ */ #include "Fields.H" #include "ThetaImplicitEM.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "WarpX.H" using warpx::fields::FieldType; @@ -111,6 +112,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real a_time, // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta} UpdateWarpXFields( m_E, theta_time ); + m_WarpX->reduced_diags->ComputeDiagsMidStep(a_step); // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate(); From 46d81e74d0bcb143b31c54a41dbcdd7e107dd510 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 2 Dec 2024 09:35:30 -0800 Subject: [PATCH 05/23] Add CI test in reduced_diags/inputs_test_3d_reduced_diags --- Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags b/Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags index dc0c57264ba..cc1b658c27f 100644 --- a/Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags +++ b/Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags @@ -68,7 +68,7 @@ photons.uz_th = 0.2 ################################# ###### REDUCED DIAGS ############ ################################# -warpx.reduced_diags_names = EP NP EF PP PF MF MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral Edotj +warpx.reduced_diags_names = EP NP EF PP PF MF PX MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral Edotj EP.type = ParticleEnergy EP.intervals = 200 EF.type = FieldEnergy @@ -79,6 +79,8 @@ PF.type = FieldMomentum PF.intervals = 200 MF.type = FieldMaximum MF.intervals = 200 +PX.type = FieldPoyntingFlux +PX.intervals = 200 FP.type = FieldProbe FP.intervals = 200 #The probe is placed at a cell center to match the value in the plotfile From bbbc5e1201944194f46148ff8cdff8cf089bf017 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 2 Dec 2024 09:46:28 -0800 Subject: [PATCH 06/23] Add temporary call to FillBoundaryB --- Source/Evolve/WarpXEvolve.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index e9540be3da7..3bbd137b1ba 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -436,6 +436,13 @@ WarpX::OneStep_nosub (Real cur_time) EvolveG(0.5_rt * dt[0], DtType::SecondHalf); EvolveB(0.5_rt * dt[0], DtType::SecondHalf); // We now have B^{n+1} + // Temporary hack so that the B guard cells are set as needed for + // the Poynting flux diagnostic. + // Maybe this should be done here anyway to make sure that everything + // is fully up to date at the end of the time step? + // This could be useful for the Python interface + FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); + if (do_pml) { DampPML(); FillBoundaryE(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points); From 77899d5b7749c08e85b467f16ecd6c7e9861e53c Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 2 Dec 2024 10:19:47 -0800 Subject: [PATCH 07/23] Some clean up --- .../Diagnostics/ReducedDiags/FieldPoyntingFlux.H | 4 ++-- .../Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H index 487ff563848..390190fa19b 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -40,8 +40,6 @@ public: */ void ComputeDiagsMidStep (int step) final; -private: - /** * \brief This function computes the electromagnetic Poynting flux, * obtained by integrating the electromagnetic Poynting flux density g = eps0 * (E x B) @@ -51,6 +49,8 @@ private: */ void ComputePoyntingFlux (int step); +private: + bool use_mid_step_value = false; }; diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 4c8abda1a58..0f000d4a80a 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -38,8 +38,7 @@ #include #include -using namespace amrex; -using warpx::fields::FieldType; +using namespace amrex::literals; FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) : ReducedDiags{rd_name} @@ -123,6 +122,7 @@ void FieldPoyntingFlux::ComputeDiagsMidStep (int step) void FieldPoyntingFlux::ComputePoyntingFlux (int step) { + using warpx::fields::FieldType; using ablastr::fields::Direction; // Check if the diags should be done @@ -171,7 +171,7 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // Get cell area const amrex::Real *dx = warpx.Geom(lev).CellSize(); - std::array dxtemp = {AMREX_D_DECL(dx[0], dx[1], dx[2])}; + std::array dxtemp = {AMREX_D_DECL(dx[0], dx[1], dx[2])}; dxtemp[face().coordDir()] = 1._rt; const amrex::Real dA = AMREX_D_TERM(dxtemp[0], *dxtemp[1], *dxtemp[2]); @@ -179,15 +179,15 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) amrex::GpuArray cc{0,0,0}; cc[face().coordDir()] = 1; - amrex::ReduceOps reduce_ops; - amrex::ReduceData reduce_data(reduce_ops); + amrex::ReduceOps reduce_ops; + amrex::ReduceData reduce_data(reduce_ops); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif // Loop over boxes, interpolate E,B data to cell face centers // and compute sum over cells of (E x B) components - for (amrex::MFIter mfi(Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi) + for (amrex::MFIter mfi(Ex, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { const amrex::Array4 & Ex_arr = Ex[mfi].array(); const amrex::Array4 & Ey_arr = Ey[mfi].array(); @@ -204,7 +204,7 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // Compute E x B reduce_ops.eval(box, reduce_data, - [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple { const amrex::Real Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); const amrex::Real Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); From fe6613453c1ed6d677fb3931e80883869c4558ab Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 5 Dec 2024 17:21:50 -0800 Subject: [PATCH 08/23] Remove the temporary FillBoundaryB as it is not needed The periodic boundaries are handled explicitly in the calculation --- Source/Evolve/WarpXEvolve.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 3bbd137b1ba..e9540be3da7 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -436,13 +436,6 @@ WarpX::OneStep_nosub (Real cur_time) EvolveG(0.5_rt * dt[0], DtType::SecondHalf); EvolveB(0.5_rt * dt[0], DtType::SecondHalf); // We now have B^{n+1} - // Temporary hack so that the B guard cells are set as needed for - // the Poynting flux diagnostic. - // Maybe this should be done here anyway to make sure that everything - // is fully up to date at the end of the time step? - // This could be useful for the Python interface - FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points); - if (do_pml) { DampPML(); FillBoundaryE(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points); From 2476d42bb462140f5286ed4edc901d7d1721d52f Mon Sep 17 00:00:00 2001 From: David Grote Date: Thu, 5 Dec 2024 17:59:34 -0800 Subject: [PATCH 09/23] Explicitly handle periodic boundary conditions --- .../ReducedDiags/FieldPoyntingFlux.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 0f000d4a80a..a4d4b68d60b 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -167,17 +167,30 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) } for (amrex::OrientationIter face; face; ++face) { + + int const face_dir = face().coordDir(); + + if (face().isHigh() && WarpX::field_boundary_hi[face_dir] == FieldBoundaryType::Periodic) { + // For upper periodic boundaries, copy the lower value instead of regenerating it. + int const iu = int(face())*3; + int const il = int(face().flip())*3; + m_data[iu+0] = m_data[il+0]; + m_data[iu+1] = m_data[il+1]; + m_data[iu+2] = m_data[il+2]; + continue; + } + amrex::Box const boundary = amrex::bdryNode(domain_box, face()); // Get cell area const amrex::Real *dx = warpx.Geom(lev).CellSize(); std::array dxtemp = {AMREX_D_DECL(dx[0], dx[1], dx[2])}; - dxtemp[face().coordDir()] = 1._rt; + dxtemp[face_dir] = 1._rt; const amrex::Real dA = AMREX_D_TERM(dxtemp[0], *dxtemp[1], *dxtemp[2]); // Node-centered in the face direction, Cell-centered in other directions amrex::GpuArray cc{0,0,0}; - cc[face().coordDir()] = 1; + cc[face_dir] = 1; amrex::ReduceOps reduce_ops; amrex::ReduceData reduce_data(reduce_ops); @@ -197,7 +210,7 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) const amrex::Array4 & Bz_arr = Bz[mfi].array(); amrex::Box box = enclosedCells(mfi.nodaltilebox()); - box.surroundingNodes(face().coordDir()); + box.surroundingNodes(face_dir); // Find the intersection with the boundary box &= boundary; From d76f23c7e70ca707fbf6f3d98d46e8ecfb2930be Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 6 Dec 2024 10:45:31 -0800 Subject: [PATCH 10/23] Code clean up --- .../ReducedDiags/FieldPoyntingFlux.cpp | 45 +++++++++---------- 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index a4d4b68d60b..485c03bb19a 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -105,17 +105,14 @@ void FieldPoyntingFlux::ComputeDiags (int step) { // This will be called at the end of the time step. Only calculate the // flux if it had not already been calculated mid step. - amrex::Print() << "ComputeDiags\n"; if (!use_mid_step_value) { ComputePoyntingFlux(step); - amrex::Print() << " calling ComputePoyntingFlux\n"; } } void FieldPoyntingFlux::ComputeDiagsMidStep (int step) { // If this is called, always use the value calculated here. - amrex::Print() << "ComputeDiagsMidStep\n"; use_mid_step_value = true; ComputePoyntingFlux(step); } @@ -136,15 +133,15 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) domain_box.surroundingNodes(); // Get MultiFab data at given refinement level - const amrex::MultiFab & Ex = *warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, lev); - const amrex::MultiFab & Ey = *warpx.m_fields.get(FieldType::Efield_aux, Direction{1}, lev); - const amrex::MultiFab & Ez = *warpx.m_fields.get(FieldType::Efield_aux, Direction{2}, lev); - const amrex::MultiFab & Bx = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{0}, lev); - const amrex::MultiFab & By = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{1}, lev); - const amrex::MultiFab & Bz = *warpx.m_fields.get(FieldType::Bfield_aux, Direction{2}, lev); + amrex::MultiFab const & Ex = *warpx.m_fields.get(FieldType::Efield_fp, Direction{0}, lev); + amrex::MultiFab const & Ey = *warpx.m_fields.get(FieldType::Efield_fp, Direction{1}, lev); + amrex::MultiFab const & Ez = *warpx.m_fields.get(FieldType::Efield_fp, Direction{2}, lev); + amrex::MultiFab const & Bx = *warpx.m_fields.get(FieldType::Bfield_fp, Direction{0}, lev); + amrex::MultiFab const & By = *warpx.m_fields.get(FieldType::Bfield_fp, Direction{1}, lev); + amrex::MultiFab const & Bz = *warpx.m_fields.get(FieldType::Bfield_fp, Direction{2}, lev); // Coarsening ratio (no coarsening) - const amrex::GpuArray cr{1,1,1}; + amrex::GpuArray const cr{1,1,1}; // Reduction component (fourth component in Array4) constexpr int comp = 0; @@ -183,10 +180,10 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) amrex::Box const boundary = amrex::bdryNode(domain_box, face()); // Get cell area - const amrex::Real *dx = warpx.Geom(lev).CellSize(); + amrex::Real const *dx = warpx.Geom(lev).CellSize(); std::array dxtemp = {AMREX_D_DECL(dx[0], dx[1], dx[2])}; dxtemp[face_dir] = 1._rt; - const amrex::Real dA = AMREX_D_TERM(dxtemp[0], *dxtemp[1], *dxtemp[2]); + amrex::Real const dA = AMREX_D_TERM(dxtemp[0], *dxtemp[1], *dxtemp[2]); // Node-centered in the face direction, Cell-centered in other directions amrex::GpuArray cc{0,0,0}; @@ -202,12 +199,12 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // and compute sum over cells of (E x B) components for (amrex::MFIter mfi(Ex, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const amrex::Array4 & Ex_arr = Ex[mfi].array(); - const amrex::Array4 & Ey_arr = Ey[mfi].array(); - const amrex::Array4 & Ez_arr = Ez[mfi].array(); - const amrex::Array4 & Bx_arr = Bx[mfi].array(); - const amrex::Array4 & By_arr = By[mfi].array(); - const amrex::Array4 & Bz_arr = Bz[mfi].array(); + amrex::Array4 const & Ex_arr = Ex[mfi].array(); + amrex::Array4 const & Ey_arr = Ey[mfi].array(); + amrex::Array4 const & Ez_arr = Ez[mfi].array(); + amrex::Array4 const & Bx_arr = Bx[mfi].array(); + amrex::Array4 const & By_arr = By[mfi].array(); + amrex::Array4 const & Bz_arr = Bz[mfi].array(); amrex::Box box = enclosedCells(mfi.nodaltilebox()); box.surroundingNodes(face_dir); @@ -219,13 +216,13 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) reduce_ops.eval(box, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple { - const amrex::Real Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); - const amrex::Real Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); - const amrex::Real Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); + amrex::Real const Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); + amrex::Real const Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); + amrex::Real const Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); - const amrex::Real Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); - const amrex::Real By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); - const amrex::Real Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); + amrex::Real const Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); + amrex::Real const By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); + amrex::Real const Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); return {Ey_cc * Bz_cc - Ez_cc * By_cc, Ez_cc * Bx_cc - Ex_cc * Bz_cc, From a96502ab2c7e2304d1ae9f56433e89a5e17fa7b6 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 6 Dec 2024 15:20:22 -0800 Subject: [PATCH 11/23] Update diagnostic Now only calculates component normal to the surface. Now works in RZ (for only one mode) --- .../ReducedDiags/FieldPoyntingFlux.cpp | 128 ++++++++++++------ 1 file changed, 89 insertions(+), 39 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 485c03bb19a..580220a6eda 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -43,18 +43,12 @@ using namespace amrex::literals; FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) : ReducedDiags{rd_name} { - // RZ coordinate is not working -#if (defined WARPX_DIM_RZ) - WARPX_ABORT_WITH_MESSAGE( - "FieldPoyntingFlux reduced diagnostics not implemented in RZ geometry"); -#endif - // Resize data array // lo and hi is 2 - // vector components is 3 // space dims is AMREX_SPACEDIM + // instantaneous and integrated is 2 // The order will be (Sx, Sy, Sz) for low faces, then high faces - m_data.resize(2*3*AMREX_SPACEDIM, 0.0_rt); + m_data.resize(2*AMREX_SPACEDIM*2, 0.0_rt); if (amrex::ParallelDescriptor::IOProcessor()) { @@ -74,26 +68,26 @@ FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) std::vector sides = {"lo", "hi"}; #if defined(WARPX_DIM_3D) - std::vector vector_coords = {"x", "y", "z"}; std::vector space_coords = {"x", "y", "z"}; #elif defined(WARPX_DIM_XZ) - std::vector vector_coords = {"x", "y", "z"}; std::vector space_coords = {"x", "z"}; #elif defined(WARPX_DIM_1D_Z) - std::vector vector_coords = {"x", "y", "z"}; std::vector space_coords = {"z"}; #elif defined(WARPX_DIM_RZ) - std::vector vector_coords = {"r", "t", "z"}; std::vector space_coords = {"r", "z"}; #endif // Only on level 0 for (int iside = 0; iside < 2; iside++) { for (int ic = 0; ic < AMREX_SPACEDIM; ic++) { - for (int iv = 0; iv < 3; iv++) { - ofs << m_sep; - ofs << "[" << c++ << "]flux_" + vector_coords[iv] + "_" + sides[iside] + "_" + space_coords[ic] +"(W)"; - }}} + ofs << m_sep; + ofs << "[" << c++ << "]outward_power_" + sides[iside] + "_" + space_coords[ic] +"(W)"; + }} + for (int iside = 0; iside < 2; iside++) { + for (int ic = 0; ic < AMREX_SPACEDIM; ic++) { + ofs << m_sep; + ofs << "[" << c++ << "]integrated_energy_loss_" + sides[iside] + "_" + space_coords[ic] +"(J)"; + }} ofs << "\n"; ofs.close(); @@ -129,6 +123,13 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // Get a reference to WarpX instance auto & warpx = WarpX::GetInstance(); + + // RZ coordinate only working with one mode +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(warpx.n_rz_azimuthal_modes == 1, + "FieldPoyntingFlux reduced diagnostics only implemented in RZ geometry for one mode"); +#endif + amrex::Box domain_box = warpx.Geom(lev).Domain(); domain_box.surroundingNodes(); @@ -142,6 +143,7 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // Coarsening ratio (no coarsening) amrex::GpuArray const cr{1,1,1}; + // Reduction component (fourth component in Array4) constexpr int comp = 0; @@ -169,11 +171,10 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) if (face().isHigh() && WarpX::field_boundary_hi[face_dir] == FieldBoundaryType::Periodic) { // For upper periodic boundaries, copy the lower value instead of regenerating it. - int const iu = int(face())*3; - int const il = int(face().flip())*3; - m_data[iu+0] = m_data[il+0]; - m_data[iu+1] = m_data[il+1]; - m_data[iu+2] = m_data[il+2]; + int const iu = int(face()); + int const il = int(face().flip()); + m_data[iu] = -m_data[il]; + m_data[iu + 2*AMREX_SPACEDIM] = -m_data[il + 2*AMREX_SPACEDIM]; continue; } @@ -189,8 +190,21 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) amrex::GpuArray cc{0,0,0}; cc[face_dir] = 1; - amrex::ReduceOps reduce_ops; - amrex::ReduceData reduce_data(reduce_ops); + // Only calculate the ExB term that is normal to the surface. + // normal_dir is the normal direction relative to the WarpX coordinates +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) + // For 2D : it is either 0, or 2 + int const normal_dir = 2*face_dir; +#elif (defined WARPX_DIM_1D_Z) + // For 1D : it is always 2 + int const normal_dir = 2; +#else + // For 3D : it is the same as the face direction + int const normal_dir = face_dir; +#endif + + amrex::ReduceOps reduce_ops; + amrex::ReduceData reduce_data(reduce_ops); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -212,29 +226,65 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) // Find the intersection with the boundary box &= boundary; +#if defined(WARPX_DIM_RZ) + // Lower corner of box physical domain + amrex::XDim3 const xyzmin = WarpX::LowerCorner(box, lev, 0._rt); + amrex::Dim3 const lo = amrex::lbound(box); + amrex::Real const dr = warpx.Geom(lev).CellSize(lev); + amrex::Real const rmin = xyzmin.x; + int const irmin = lo.x; +#endif + + auto area_factor = [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::ignore_unused(i,j,k); +#if defined WARPX_DIM_RZ + amrex::Real r; + if (normal_dir == 0) { + r = rmin + (i - irmin)*dr; + } else { + r = rmin + (i + 0.5_rt - irmin)*dr; + } + return 2.*MathConst::pi*r; +#else + return 1._rt; +#endif + }; + // Compute E x B reduce_ops.eval(box, reduce_data, - [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::GpuTuple { - amrex::Real const Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); - amrex::Real const Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); - amrex::Real const Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); - - amrex::Real const Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); - amrex::Real const By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); - amrex::Real const Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); - - return {Ey_cc * Bz_cc - Ez_cc * By_cc, - Ez_cc * Bx_cc - Ex_cc * Bz_cc, - Ex_cc * By_cc - Ey_cc * Bx_cc}; + amrex::Real Ex_cc = 0._rt, Ey_cc = 0._rt, Ez_cc = 0._rt; + amrex::Real Bx_cc = 0._rt, By_cc = 0._rt, Bz_cc = 0._rt; + + if (normal_dir == 1 || normal_dir == 2) { + Ex_cc = ablastr::coarsen::sample::Interp(Ex_arr, Ex_stag, cc, cr, i, j, k, comp); + Bx_cc = ablastr::coarsen::sample::Interp(Bx_arr, Bx_stag, cc, cr, i, j, k, comp); + } + + if (normal_dir == 0 || normal_dir == 2) { + Ey_cc = ablastr::coarsen::sample::Interp(Ey_arr, Ey_stag, cc, cr, i, j, k, comp); + By_cc = ablastr::coarsen::sample::Interp(By_arr, By_stag, cc, cr, i, j, k, comp); + } + if (normal_dir == 0 || normal_dir == 1) { + Ez_cc = ablastr::coarsen::sample::Interp(Ez_arr, Ez_stag, cc, cr, i, j, k, comp); + Bz_cc = ablastr::coarsen::sample::Interp(Bz_arr, Bz_stag, cc, cr, i, j, k, comp); + } + + amrex::Real const af = area_factor(i,j,k); + if (normal_dir == 0) { return af*(Ey_cc * Bz_cc - Ez_cc * By_cc); } + else if (normal_dir == 1) { return af*(Ez_cc * Bx_cc - Ex_cc * Bz_cc); } + else { return af*(Ex_cc * By_cc - Ey_cc * Bx_cc); } }); } + int const sign = (face().isLow() ? -1 : 1); auto r = reduce_data.value(); - int const ii = int(face())*3; - m_data[ii+0] = amrex::get<0>(r)/PhysConst::mu0*dA; - m_data[ii+1] = amrex::get<1>(r)/PhysConst::mu0*dA; - m_data[ii+2] = amrex::get<2>(r)/PhysConst::mu0*dA; + int const ii = int(face()); + m_data[ii] = sign*amrex::get<0>(r)/PhysConst::mu0*dA; + + amrex::Real const dt = warpx.getdt(lev); + m_data[ii + 2*AMREX_SPACEDIM] += m_data[ii]*dt; } From 86d72ad3a1ab139cfe3de45ab2e14194b7480cd9 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 9 Dec 2024 11:18:25 -0800 Subject: [PATCH 12/23] Add checkpoint/restart capability for time integrated diagnostic --- .../FlushFormats/FlushFormatCheckpoint.H | 2 ++ .../FlushFormats/FlushFormatCheckpoint.cpp | 12 +++++++ .../ReducedDiags/FieldPoyntingFlux.H | 4 +++ .../ReducedDiags/FieldPoyntingFlux.cpp | 33 +++++++++++++++++++ .../ReducedDiags/MultiReducedDiags.H | 7 ++++ .../ReducedDiags/MultiReducedDiags.cpp | 23 +++++++++++++ .../Diagnostics/ReducedDiags/ReducedDiags.H | 14 ++++++++ .../Diagnostics/ReducedDiags/ReducedDiags.cpp | 14 ++++++++ Source/Diagnostics/WarpXIO.cpp | 3 ++ 9 files changed, 112 insertions(+) diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H index 5c26ac97f61..429855a8196 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H @@ -34,6 +34,8 @@ class FlushFormatCheckpoint final : public FlushFormatPlotfile const amrex::Vector& particle_diags) const; void WriteDMaps (const std::string& dir, int nlev) const; + + void WriteReducedDiagsData (std::string const & dir) const; }; #endif // WARPX_FLUSHFORMATCHECKPOINT_H_ diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 4d721dd6abe..2d046bd7f45 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -5,6 +5,7 @@ # include "BoundaryConditions/PML_RZ.H" #endif #include "Diagnostics/ParticleDiag/ParticleDiag.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include "Fields.H" #include "Particles/WarpXParticleContainer.H" #include "Utils/TextMsg.H" @@ -171,6 +172,8 @@ FlushFormatCheckpoint::WriteToFile ( WriteDMaps(checkpointname, nlev); + WriteReducedDiagsData(checkpointname); + VisMF::SetHeaderVersion(current_version); } @@ -242,3 +245,12 @@ FlushFormatCheckpoint::WriteDMaps (const std::string& dir, int nlev) const } } } + +void +FlushFormatCheckpoint::WriteReducedDiagsData (std::string const & dir) const +{ + if (ParallelDescriptor::IOProcessor()) { + auto & warpx = WarpX::GetInstance(); + warpx.reduced_diags->WriteCheckpointData(dir); + } +} diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H index 390190fa19b..2fcc898dd65 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -49,6 +49,10 @@ public: */ void ComputePoyntingFlux (int step); + void WriteCheckpointData (std::string const & dir); + + void ReadCheckpointData (std::string const & dir); + private: bool use_mid_step_value = false; diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 580220a6eda..69b4c95e149 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -291,3 +291,36 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) amrex::ParallelDescriptor::ReduceRealSum(m_data.data(), static_cast(m_data.size())); } + +void +FieldPoyntingFlux::WriteCheckpointData (std::string const & dir) +{ + // Write out the current values of the time integrated data + std::ofstream chkfile{dir + "/FieldPoyntingFlux_data.txt", std::ofstream::out}; + if (!chkfile.good()) { + WARPX_ABORT_WITH_MESSAGE("FieldPoyntingFlux::WriteCheckpointData: could not open file for writing checkpoint data"); + } + + for (int i=0; i < 2*AMREX_SPACEDIM; i++) { + chkfile << m_data[2*AMREX_SPACEDIM + i] << "\n"; + } +} + +void +FieldPoyntingFlux::ReadCheckpointData (std::string const & dir) +{ + // Read in the current values of the time integrated data + std::ifstream chkfile{dir + "/FieldPoyntingFlux_data.txt", std::ifstream::in}; + if (!chkfile.good()) { + WARPX_ABORT_WITH_MESSAGE("FieldPoyntingFlux::ReadCheckpointData: could not open file for reading checkpoint data"); + } + + for (int i=0; i < 2*AMREX_SPACEDIM; i++) { + amrex::Real data; + if (chkfile >> data) { + m_data[2*AMREX_SPACEDIM + i] = data; + } else { + WARPX_ABORT_WITH_MESSAGE("FieldPoyntingFlux::ReadCheckpointData: could not read in time integrated data"); + } + } +} diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H index 64157e14791..5a782db7118 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H @@ -57,6 +57,13 @@ public: * @param[in] step current iteration time */ void WriteToFile (int step); + /** \brief Loop over all ReducedDiags and call their WriteCheckpointData + * @param[in] dir checkpoint directory */ + void WriteCheckpointData (std::string const & dir); + + /** \brief Loop over all ReducedDiags and call their ReadCheckpointData + * @param[in] dir checkpoint directory */ + void ReadCheckpointData (std::string const & dir); }; #endif diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 27faaf9f3ca..0ce18174111 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -158,3 +158,26 @@ void MultiReducedDiags::WriteToFile (int step) // end loop over all reduced diags } // end void MultiReducedDiags::WriteToFile + +void MultiReducedDiags::WriteCheckpointData (std::string const & dir) +{ + // Only the I/O rank does + if ( !ParallelDescriptor::IOProcessor() ) { return; } + + // loop over all reduced diags + for (int i_rd = 0; i_rd < static_cast(m_rd_names.size()); ++i_rd) + { + m_multi_rd[i_rd]->WriteCheckpointData(dir); + } + // end loop over all reduced diags +} + +void MultiReducedDiags::ReadCheckpointData (std::string const & dir) +{ + // loop over all reduced diags + for (int i_rd = 0; i_rd < static_cast(m_rd_names.size()); ++i_rd) + { + m_multi_rd[i_rd]->ReadCheckpointData(dir); + } + // end loop over all reduced diags +} diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.H b/Source/Diagnostics/ReducedDiags/ReducedDiags.H index a5c787e636a..a32de30cc6f 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.H +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.H @@ -97,6 +97,20 @@ public: */ virtual void WriteToFile (int step) const; + /** + * \brief Write out checkpoint related data + * + * \param[in] dir Directory where checkpoint data is written + */ + virtual void WriteCheckpointData (std::string const & dir); + + /** + * \brief Read in checkpoint related data + * + * \param[in] dir Directory where checkpoint data was written + */ + virtual void ReadCheckpointData (std::string const & dir); + /** * This function queries deprecated input parameters and aborts * the run if one of them is specified. diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp index 7abf91295db..b0e20584a12 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp @@ -99,6 +99,20 @@ void ReducedDiags::ComputeDiagsMidStep (int /*step*/) // (instead of at the end of the step). } +void ReducedDiags::WriteCheckpointData (std::string const & /*dir*/) +{ + // Defines an empty function WriteCheckpointData() to be overwritten if needed. + // Function used to write out and data needed by the diagnostic in + // the checkpoint. +} + +void ReducedDiags::ReadCheckpointData (std::string const & /*dir*/) +{ + // Defines an empty function ReadCheckpointData() to be overwritten if needed. + // Function used to read in any data that was written out in the checkpoint + // when doing a restart. +} + void ReducedDiags::BackwardCompatibility () const { const amrex::ParmParse pp_rd_name(m_rd_name); diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 43415daf151..acc5e7d8d15 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -19,6 +19,7 @@ #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" #include "Diagnostics/MultiDiagnostics.H" +#include "Diagnostics/ReducedDiags/MultiReducedDiags.H" #include #include @@ -400,6 +401,8 @@ WarpX::InitFromCheckpoint () if (EB::enabled()) { InitializeEBGridData(maxLevel()); } + reduced_diags->ReadCheckpointData(restart_chkfile); + // Initialize particles mypc->AllocData(); mypc->Restart(restart_chkfile); From 2fc79fe87c5035499e950af17413d1a9c538f19a Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 9 Dec 2024 11:41:24 -0800 Subject: [PATCH 13/23] Set precision on checkpoint data --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 69b4c95e149..9416633c56d 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -301,6 +301,8 @@ FieldPoyntingFlux::WriteCheckpointData (std::string const & dir) WARPX_ABORT_WITH_MESSAGE("FieldPoyntingFlux::WriteCheckpointData: could not open file for writing checkpoint data"); } + chkfile.precision(17); + for (int i=0; i < 2*AMREX_SPACEDIM; i++) { chkfile << m_data[2*AMREX_SPACEDIM + i] << "\n"; } From bbf8a557580079da47ab4b83944619b5637158ea Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 9 Dec 2024 11:50:46 -0800 Subject: [PATCH 14/23] Add documentation --- Docs/source/usage/parameters.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index c3756d44d96..c0a1d16c0a7 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -3150,6 +3150,12 @@ This shifts analysis from post-processing to runtime calculation of reduction op Note that the fields are averaged on the cell centers before their maximum values are computed. + * ``FieldPoyntingFlux`` + Integrates the normal Poynting flux over each domain boundary surface and also integrates the flux over time. + This provides the power and total energy loss into or out of the simulation domain. + The output columns are the flux for each dimension on the lower boundaries, then the higher boundaries, + then the integrated energy loss for each dimension on the the lower and higher boundaries. + * ``FieldProbe`` This type computes the value of each component of the electric and magnetic fields and of the Poynting vector (a measure of electromagnetic flux) at points in the domain. From f4e5563ccd72e1965e33f5958cf9666f924cdaca Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 9 Dec 2024 12:01:17 -0800 Subject: [PATCH 15/23] Calculate every step independently of the intervals --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 9416633c56d..adea8b88a0b 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -116,8 +116,8 @@ void FieldPoyntingFlux::ComputePoyntingFlux (int step) using warpx::fields::FieldType; using ablastr::fields::Direction; - // Check if the diags should be done - if (!m_intervals.contains(step+1)) { return; } + // Note that this is calculated every step to get the + // full resolution on the integrated data int const lev = 0; From 761e77d2a0b61ea94cc9940df07ec3a76b616b07 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 9 Dec 2024 13:24:24 -0800 Subject: [PATCH 16/23] Clean up unused parameter --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H | 2 +- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H index 2fcc898dd65..cf9a0feaeaf 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -47,7 +47,7 @@ public: * * \param[in] step current time step */ - void ComputePoyntingFlux (int step); + void ComputePoyntingFlux (); void WriteCheckpointData (std::string const & dir); diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index adea8b88a0b..fd0974a1ab3 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -100,7 +100,7 @@ void FieldPoyntingFlux::ComputeDiags (int step) // This will be called at the end of the time step. Only calculate the // flux if it had not already been calculated mid step. if (!use_mid_step_value) { - ComputePoyntingFlux(step); + ComputePoyntingFlux(); } } @@ -108,10 +108,10 @@ void FieldPoyntingFlux::ComputeDiagsMidStep (int step) { // If this is called, always use the value calculated here. use_mid_step_value = true; - ComputePoyntingFlux(step); + ComputePoyntingFlux(); } -void FieldPoyntingFlux::ComputePoyntingFlux (int step) +void FieldPoyntingFlux::ComputePoyntingFlux () { using warpx::fields::FieldType; using ablastr::fields::Direction; From e84d1e2fb6f99bc22d855b60cd62d74e432af3a9 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 10 Dec 2024 08:52:48 -0800 Subject: [PATCH 17/23] Remove unused argument --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index fd0974a1ab3..365fb6b9a9e 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -95,7 +95,7 @@ FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) } } -void FieldPoyntingFlux::ComputeDiags (int step) +void FieldPoyntingFlux::ComputeDiags (int /*step*/) { // This will be called at the end of the time step. Only calculate the // flux if it had not already been calculated mid step. @@ -104,7 +104,7 @@ void FieldPoyntingFlux::ComputeDiags (int step) } } -void FieldPoyntingFlux::ComputeDiagsMidStep (int step) +void FieldPoyntingFlux::ComputeDiagsMidStep (int /*step*/) { // If this is called, always use the value calculated here. use_mid_step_value = true; From a7e9ddcaaf65b6fe37325dc6a0c73b7916b07154 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 10 Dec 2024 09:32:26 -0800 Subject: [PATCH 18/23] Made methods final --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H index cf9a0feaeaf..ca705bb4b02 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -49,9 +49,9 @@ public: */ void ComputePoyntingFlux (); - void WriteCheckpointData (std::string const & dir); + void WriteCheckpointData (std::string const & dir) final; - void ReadCheckpointData (std::string const & dir); + void ReadCheckpointData (std::string const & dir) final; private: From 763d4949e4ce0fa354efe1015de3e3e090dd3be1 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 10 Dec 2024 12:57:01 -0800 Subject: [PATCH 19/23] Fix literal suffix and comment --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 365fb6b9a9e..3157d5d0490 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -47,7 +47,8 @@ FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) // lo and hi is 2 // space dims is AMREX_SPACEDIM // instantaneous and integrated is 2 - // The order will be (Sx, Sy, Sz) for low faces, then high faces + // The order will be outward flux for low faces, then high faces, + // energy loss for low faces, then high faces m_data.resize(2*AMREX_SPACEDIM*2, 0.0_rt); if (amrex::ParallelDescriptor::IOProcessor()) @@ -244,7 +245,7 @@ void FieldPoyntingFlux::ComputePoyntingFlux () } else { r = rmin + (i + 0.5_rt - irmin)*dr; } - return 2.*MathConst::pi*r; + return 2._rt*MathConst::pi*r; #else return 1._rt; #endif From 7d487c975bd4d9cebd9112452957aea6a8a88727 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 13 Dec 2024 16:14:56 -0800 Subject: [PATCH 20/23] Fix the type of the boundary box --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index 3157d5d0490..b610b66f760 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -225,7 +225,9 @@ void FieldPoyntingFlux::ComputePoyntingFlux () box.surroundingNodes(face_dir); // Find the intersection with the boundary - box &= boundary; + // boundary needs to have the same type as box + amrex::Box const boundary_matched = amrex::convert(boundary, box.ixType()); + box &= boundary_matched; #if defined(WARPX_DIM_RZ) // Lower corner of box physical domain From 2203ed9b1054bb1eb09fb389f95843b112ebfab8 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 13 Dec 2024 16:27:30 -0800 Subject: [PATCH 21/23] Fix integrated term for parallel --- Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp index b610b66f760..f760516f2b9 100644 --- a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -286,12 +286,14 @@ void FieldPoyntingFlux::ComputePoyntingFlux () int const ii = int(face()); m_data[ii] = sign*amrex::get<0>(r)/PhysConst::mu0*dA; - amrex::Real const dt = warpx.getdt(lev); - m_data[ii + 2*AMREX_SPACEDIM] += m_data[ii]*dt; - } - amrex::ParallelDescriptor::ReduceRealSum(m_data.data(), static_cast(m_data.size())); + amrex::ParallelDescriptor::ReduceRealSum(m_data.data(), 2*AMREX_SPACEDIM); + + amrex::Real const dt = warpx.getdt(lev); + for (int ii=0 ; ii < 2*AMREX_SPACEDIM ; ii++) { + m_data[ii + 2*AMREX_SPACEDIM] += m_data[ii]*dt; + } } From 0e4a4175105f5d466cf3fccbf09b636819502cfe Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 13 Dec 2024 17:06:19 -0800 Subject: [PATCH 22/23] Add stringent CI test that checks for machine level energy accouting --- Examples/Tests/pec/CMakeLists.txt | 20 ++++++ .../pec/analysis_pec_insulator_implicit.py | 67 +++++++++++++++++ ...nputs_test_2d_pec_field_insulator_implicit | 72 +++++++++++++++++++ ...st_2d_pec_field_insulator_implicit_restart | 5 ++ .../test_2d_pec_field_insulator_implicit.json | 14 ++++ ..._pec_field_insulator_implicit_restart.json | 14 ++++ 6 files changed, 192 insertions(+) create mode 100755 Examples/Tests/pec/analysis_pec_insulator_implicit.py create mode 100644 Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit create mode 100644 Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart create mode 100644 Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit.json create mode 100644 Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit_restart.json diff --git a/Examples/Tests/pec/CMakeLists.txt b/Examples/Tests/pec/CMakeLists.txt index e0bab40d058..bc5807b2069 100644 --- a/Examples/Tests/pec/CMakeLists.txt +++ b/Examples/Tests/pec/CMakeLists.txt @@ -40,3 +40,23 @@ add_warpx_test( diags/diag1000010 # output OFF # dependency ) + +add_warpx_test( + test_2d_pec_field_insulator_implicit # name + 2 # dims + 2 # nprocs + inputs_test_2d_pec_field_insulator_implicit # inputs + analysis_pec_insulator_implicit.py # analysis + diags/diag1000020 # output + OFF # dependency +) + +add_warpx_test( + test_2d_pec_field_insulator_implicit_restart # name + 2 # dims + 2 # nprocs + inputs_test_2d_pec_field_insulator_implicit_restart # inputs + analysis_pec_insulator_implicit.py # analysis + diags/diag1000020 # output + OFF # dependency +) diff --git a/Examples/Tests/pec/analysis_pec_insulator_implicit.py b/Examples/Tests/pec/analysis_pec_insulator_implicit.py new file mode 100755 index 00000000000..bdd5e291133 --- /dev/null +++ b/Examples/Tests/pec/analysis_pec_insulator_implicit.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 + +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL +# +# This is a script that analyses the simulation results from +# the scripts `inputs_test_2d_pec_field_insulator_implicit` and +# `inputs_test_2d_pec_field_insulator_implicit_restart`. +# The scripts model an insulator boundary condition on part of the +# upper x boundary that pushes B field into the domain. The implicit +# solver is used, converging to machine tolerance. The energy accounting +# should be exact to machine precision, so that the energy is the system +# should be the same as the amount of energy pushed in from the boundary. +# This is checked using the FieldEnergy and FieldPoyntingFlux reduced +# diagnostics. +import os +import sys + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +sys.path.insert(1, "../../../../warpx/Regression/Checksum/") +from checksumAPI import evaluate_checksum + +# this will be the name of the plot file +fn = sys.argv[1] + +EE = np.loadtxt(f"{fn}/../reducedfiles/fieldenergy.txt", skiprows=1) +SS = np.loadtxt(f"{fn}/../reducedfiles/poyntingflux.txt", skiprows=1) +SSsum = SS[:,2:6].sum(1) +EEloss = SS[:,7:].sum(1) + +dt = EE[1,1] + +fig, ax = plt.subplots() +ax.plot(EE[:,0], EE[:,2], label="field energy") +ax.plot(SS[:,0], -EEloss, label="-flux*dt") +ax.legend() +ax.set_xlabel("time (s)") +ax.set_ylabel("energy (J)") +fig.savefig("energy_history.png") + +fig, ax = plt.subplots() +ax.plot(EE[:,0], (EE[:,2] + EEloss)/EE[:,2].max()) +ax.set_xlabel("time (s)") +ax.set_ylabel("energy difference/max energy (1)") +fig.savefig("energy_difference.png") + +tolerance_rel = 1.e-13 + +energy_difference_fraction = np.abs((EE[:,2] + EEloss)/EE[:,2].max()).max() +print(f"energy accounting error = {energy_difference_fraction}") +print(f"tolerance_rel = {tolerance_rel}") + +assert energy_difference_fraction < tolerance_rel + +# compare checksums +evaluate_checksum( + test_name=os.path.split(os.getcwd())[1], + output_file=sys.argv[1], +) diff --git a/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit new file mode 100644 index 00000000000..773de2bd1f4 --- /dev/null +++ b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit @@ -0,0 +1,72 @@ +# Maximum number of time steps +max_step = 20 + +# number of grid points +amr.n_cell = 32 32 +amr.blocking_factor = 16 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +# Geometry +geometry.dims = 2 +geometry.prob_lo = 0. 2.e-2 # physical domain +geometry.prob_hi = 1.e-2 3.e-2 + +# Boundary condition +boundary.field_lo = neumann periodic +boundary.field_hi = pec_insulator periodic + +insulator.area_x_hi(y,z) = (2.25e-2 <= z and z <= 2.75e-2) +insulator.By_x_hi(y,z,t) = min(t/1.0e-12,1)*1.e1*3.3e-4 + +warpx.serialize_initial_conditions = 1 + +# Implicit setup +warpx.const_dt = 7.37079480234276e-13/2. + +algo.maxwell_solver = Yee +algo.evolve_scheme = "theta_implicit_em" +#algo.evolve_scheme = "semi_implicit_em" + +implicit_evolve.theta = 0.5 +#implicit_evolve.max_particle_iterations = 21 +#implicit_evolve.particle_tolerance = 1.0e-12 + +implicit_evolve.nonlinear_solver = "picard" +picard.verbose = true +picard.max_iterations = 25 +picard.relative_tolerance = 0.0 +picard.absolute_tolerance = 0.0 +picard.require_convergence = false + +#implicit_evolve.nonlinear_solver = "newton" +#newton.verbose = true +#newton.max_iterations = 20 +#newton.relative_tolerance = 1.0e-20 +#newton.absolute_tolerance = 0.0 +#newton.require_convergence = false + +#gmres.verbose_int = 2 +#gmres.max_iterations = 1000 +#gmres.relative_tolerance = 1.0e-20 +#gmres.absolute_tolerance = 0.0 + +# Verbosity +warpx.verbose = 1 + +# Diagnostics +diagnostics.diags_names = diag1 chk +diag1.intervals = 20 +diag1.diag_type = Full + +chk.intervals = 10 +chk.diag_type = Full +chk.format = checkpoint + +warpx.reduced_diags_names = fieldenergy poyntingflux +poyntingflux.type = FieldPoyntingFlux +poyntingflux.intervals = 1 +fieldenergy.type = FieldEnergy +fieldenergy.intervals = 1 + diff --git a/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart new file mode 100644 index 00000000000..35b78d01acd --- /dev/null +++ b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart @@ -0,0 +1,5 @@ +# base input parameters +FILE = inputs_test_2d_pec_field_insulator_implicit + +# test input parameters +amr.restart = "../test_2d_pec_field_insulator_implicit/diags/chk000010" diff --git a/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit.json b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit.json new file mode 100644 index 00000000000..fcb3081f6ae --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit.json @@ -0,0 +1,14 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.35907571934346943, + "Bz": 0.0, + "Ex": 36840284.366667606, + "Ey": 0.0, + "Ez": 107777138.0847348, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + } +} + diff --git a/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit_restart.json b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit_restart.json new file mode 100644 index 00000000000..fcb3081f6ae --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit_restart.json @@ -0,0 +1,14 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.35907571934346943, + "Bz": 0.0, + "Ex": 36840284.366667606, + "Ey": 0.0, + "Ez": 107777138.0847348, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + } +} + From cece997ae0d77d9ce0e4eb3fd27cbe02f5fddb0d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 14 Dec 2024 01:07:05 +0000 Subject: [PATCH 23/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../Tests/pec/analysis_pec_insulator_implicit.py | 16 ++++++++-------- .../inputs_test_2d_pec_field_insulator_implicit | 1 - 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/Examples/Tests/pec/analysis_pec_insulator_implicit.py b/Examples/Tests/pec/analysis_pec_insulator_implicit.py index bdd5e291133..4eda456c23e 100755 --- a/Examples/Tests/pec/analysis_pec_insulator_implicit.py +++ b/Examples/Tests/pec/analysis_pec_insulator_implicit.py @@ -33,28 +33,28 @@ EE = np.loadtxt(f"{fn}/../reducedfiles/fieldenergy.txt", skiprows=1) SS = np.loadtxt(f"{fn}/../reducedfiles/poyntingflux.txt", skiprows=1) -SSsum = SS[:,2:6].sum(1) -EEloss = SS[:,7:].sum(1) +SSsum = SS[:, 2:6].sum(1) +EEloss = SS[:, 7:].sum(1) -dt = EE[1,1] +dt = EE[1, 1] fig, ax = plt.subplots() -ax.plot(EE[:,0], EE[:,2], label="field energy") -ax.plot(SS[:,0], -EEloss, label="-flux*dt") +ax.plot(EE[:, 0], EE[:, 2], label="field energy") +ax.plot(SS[:, 0], -EEloss, label="-flux*dt") ax.legend() ax.set_xlabel("time (s)") ax.set_ylabel("energy (J)") fig.savefig("energy_history.png") fig, ax = plt.subplots() -ax.plot(EE[:,0], (EE[:,2] + EEloss)/EE[:,2].max()) +ax.plot(EE[:, 0], (EE[:, 2] + EEloss) / EE[:, 2].max()) ax.set_xlabel("time (s)") ax.set_ylabel("energy difference/max energy (1)") fig.savefig("energy_difference.png") -tolerance_rel = 1.e-13 +tolerance_rel = 1.0e-13 -energy_difference_fraction = np.abs((EE[:,2] + EEloss)/EE[:,2].max()).max() +energy_difference_fraction = np.abs((EE[:, 2] + EEloss) / EE[:, 2].max()).max() print(f"energy accounting error = {energy_difference_fraction}") print(f"tolerance_rel = {tolerance_rel}") diff --git a/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit index 773de2bd1f4..be7668b7d58 100644 --- a/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit +++ b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit @@ -69,4 +69,3 @@ poyntingflux.type = FieldPoyntingFlux poyntingflux.intervals = 1 fieldenergy.type = FieldEnergy fieldenergy.intervals = 1 -