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. 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..4eda456c23e --- /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.0e-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..be7668b7d58 --- /dev/null +++ b/Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit @@ -0,0 +1,71 @@ +# 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/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 diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index c5946376d52..76298e8479b 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -3996,6 +3996,7 @@ def __init__( "FieldEnergy", "FieldMomentum", "FieldMaximum", + "FieldPoyntingFlux", "RhoMaximum", "ParticleNumber", "LoadBalanceCosts", 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 + } +} + 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/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..ca705bb4b02 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H @@ -0,0 +1,62 @@ +/* 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 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; + + /** + * \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 ComputePoyntingFlux (); + + void WriteCheckpointData (std::string const & dir) final; + + void ReadCheckpointData (std::string const & dir) final; + +private: + + bool use_mid_step_value = false; + +}; + +#endif diff --git a/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp new file mode 100644 index 00000000000..f760516f2b9 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.cpp @@ -0,0 +1,333 @@ +/* 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::literals; + +FieldPoyntingFlux::FieldPoyntingFlux (const std::string& rd_name) + : ReducedDiags{rd_name} +{ + // Resize data array + // lo and hi is 2 + // space dims is AMREX_SPACEDIM + // instantaneous and integrated is 2 + // 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()) + { + 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"}; + +#if defined(WARPX_DIM_3D) + std::vector space_coords = {"x", "y", "z"}; +#elif defined(WARPX_DIM_XZ) + std::vector space_coords = {"x", "z"}; +#elif defined(WARPX_DIM_1D_Z) + std::vector space_coords = {"z"}; +#elif defined(WARPX_DIM_RZ) + 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++) { + 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(); + } + } +} + +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(); + } +} + +void FieldPoyntingFlux::ComputeDiagsMidStep (int /*step*/) +{ + // If this is called, always use the value calculated here. + use_mid_step_value = true; + ComputePoyntingFlux(); +} + +void FieldPoyntingFlux::ComputePoyntingFlux () +{ + using warpx::fields::FieldType; + using ablastr::fields::Direction; + + // Note that this is calculated every step to get the + // full resolution on the integrated data + + int const lev = 0; + + // 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(); + + // Get MultiFab data at given refinement level + 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) + amrex::GpuArray const 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) { + + 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()); + 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; + } + + amrex::Box const boundary = amrex::bdryNode(domain_box, face()); + + // Get cell area + 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; + 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}; + cc[face_dir] = 1; + + // 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()) +#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, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + 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); + + // Find the intersection with the 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 + 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._rt*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::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()); + m_data[ii] = sign*amrex::get<0>(r)/PhysConst::mu0*dA; + + } + + 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; + } + +} + +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"); + } + + chkfile.precision(17); + + 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/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.H b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H index 1a2f51794c6..5a782db7118 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H @@ -49,10 +49,21 @@ 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); + /** \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 5035eac58a8..0ce18174111 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);}}, @@ -124,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) { @@ -142,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 2c942e1df6d..a32de30cc6f 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 * @@ -90,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 a3529cd305d..b0e20584a12 100644 --- a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp @@ -92,6 +92,27 @@ 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::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); diff --git a/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp b/Source/FieldSolver/ImplicitSolvers/SemiImplicitEM.cpp index 41fdf515581..bf8441e1992 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 start_time, // Update WarpX owned Efield_fp to t_{n+1/2} m_WarpX->SetElectricFieldAndApplyBCs( 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/StrangImplicitSpectralEM.cpp b/Source/FieldSolver/ImplicitSolvers/StrangImplicitSpectralEM.cpp index cd672f18f98..b8be6b93c63 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 start_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 aa6ee63f7df..1e6596f5eaa 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; @@ -109,6 +110,7 @@ void ThetaImplicitEM::OneStep ( const amrex::Real start_time, // Update WarpX owned Efield_fp and Bfield_fp to t_{n+theta} UpdateWarpXFields( m_E, start_time ); + m_WarpX->reduced_diags->ComputeDiagsMidStep(a_step); // Advance particles from time n+1/2 to time n+1 m_WarpX->FinishImplicitParticleUpdate();