diff options
Diffstat (limited to 'Source/Diagnostics/FieldIO.cpp')
| -rw-r--r-- | Source/Diagnostics/FieldIO.cpp | 184 |
1 files changed, 157 insertions, 27 deletions
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 5c523804e..a1fd61b25 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -204,6 +204,20 @@ WriteOpenPMDFields( const std::string& filename, } #endif // WARPX_USE_OPENPMD +void +ConstructTotalRZField(std::array< std::unique_ptr<MultiFab>, 3 >& mf_total, + const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field) +{ + // Sum over the real components, giving quantity at theta=0 + MultiFab::Copy(*mf_total[0], *vector_field[0], 0, 0, 1, vector_field[0]->nGrowVect()); + MultiFab::Copy(*mf_total[1], *vector_field[1], 0, 0, 1, vector_field[1]->nGrowVect()); + MultiFab::Copy(*mf_total[2], *vector_field[2], 0, 0, 1, vector_field[2]->nGrowVect()); + for (int ic=2 ; ic < vector_field[0]->nComp() ; ic += 2) { + MultiFab::Add(*mf_total[0], *vector_field[0], ic, 0, 1, vector_field[0]->nGrowVect()); + MultiFab::Add(*mf_total[1], *vector_field[1], ic, 0, 1, vector_field[1]->nGrowVect()); + MultiFab::Add(*mf_total[2], *vector_field[2], ic, 0, 1, vector_field[2]->nGrowVect()); + } +} void PackPlotDataPtrs (Vector<const MultiFab*>& pmf, @@ -228,6 +242,7 @@ PackPlotDataPtrs (Vector<const MultiFab*>& pmf, void AverageAndPackVectorField( MultiFab& mf_avg, const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field, + const DistributionMapping& dm, const int dcomp, const int ngrow ) { // The object below is temporary, and is needed because @@ -256,23 +271,101 @@ AverageAndPackVectorField( MultiFab& mf_avg, // - Face centered, in the same way as B on a Yee grid } else if ( vector_field[0]->is_nodal(0) ){ - PackPlotDataPtrs(srcmf, vector_field); - amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); + // Note that average_face_to_cellcenter operates only on the number of + // arrays equal to the number of dimensions. So, for 2D, PackPlotDataPtrs + // packs in the x and z (or r and z) arrays, which are then cell averaged. + // The Copy code then copies the z from the 2nd to the 3rd field, + // and copies over directly the y (or theta) component (which is + // already cell centered). + if (vector_field[0]->nComp() > 1) { +#ifdef WARPX_RZ + // When there are more than one components, the total + // fields needs to be constructed in temporary MultiFabs + // Note that mf_total is declared in the same way as + // vector_field so that it can be passed into PackPlotDataPtrs. + std::array<std::unique_ptr<MultiFab>,3> mf_total; + mf_total[0].reset(new MultiFab(vector_field[0]->boxArray(), dm, 1, vector_field[0]->nGrowVect())); + mf_total[1].reset(new MultiFab(vector_field[1]->boxArray(), dm, 1, vector_field[1]->nGrowVect())); + mf_total[2].reset(new MultiFab(vector_field[2]->boxArray(), dm, 1, vector_field[2]->nGrowVect())); + ConstructTotalRZField(mf_total, vector_field); + PackPlotDataPtrs(srcmf, mf_total); + amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); + MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); + MultiFab::Copy( mf_avg, *mf_total[1], 0, dcomp+1, 1, ngrow); + // Also add the real and imaginary parts of each mode. + /* Not supported yet + for (int i=0 ; i < vector_field[0]->nComp() ; i++) { + MultiFab v_comp0(*vector_field[0], amrex::make_alias, i, 1); + MultiFab v_comp1(*vector_field[1], amrex::make_alias, i, 1); + MultiFab v_comp2(*vector_field[2], amrex::make_alias, i, 1); + srcmf[0] = &v_comp0; + srcmf[1] = &v_comp2; + int id = dcomp + 3*(i + 1); + amrex::average_face_to_cellcenter( mf_avg, id, srcmf, ngrow); + MultiFab::Copy( mf_avg, mf_avg, id+1, id+2, 1, ngrow); + MultiFab::Copy( mf_avg, v_comp1, 0, id+1, 1, ngrow); + } + */ +#else + amrex::Abort("AverageAndPackVectorField not implemented for ncomp > 1"); +#endif + } else { + PackPlotDataPtrs(srcmf, vector_field); + amrex::average_face_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); - MultiFab::Copy( mf_avg, *vector_field[1], 0, dcomp+1, 1, ngrow); + MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); + MultiFab::Copy( mf_avg, *vector_field[1], 0, dcomp+1, 1, ngrow); #endif + } // - Edge centered, in the same way as E on a Yee grid } else if ( !vector_field[0]->is_nodal(0) ){ - PackPlotDataPtrs(srcmf, vector_field); - amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); + // See comment above, though here, the y (or theta) component + // has node centering. + if (vector_field[0]->nComp() > 1) { +#ifdef WARPX_RZ + // When there are more than one components, the total + // fields needs to be constructed in temporary MultiFabs + // Note that mf_total is declared in the same way as + // vector_field so that it can be passed into PackPlotDataPtrs. + std::array<std::unique_ptr<MultiFab>,3> mf_total; + mf_total[0].reset(new MultiFab(vector_field[0]->boxArray(), dm, 1, vector_field[0]->nGrowVect())); + mf_total[1].reset(new MultiFab(vector_field[1]->boxArray(), dm, 1, vector_field[1]->nGrowVect())); + mf_total[2].reset(new MultiFab(vector_field[2]->boxArray(), dm, 1, vector_field[2]->nGrowVect())); + ConstructTotalRZField(mf_total, vector_field); + PackPlotDataPtrs(srcmf, mf_total); + amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); + MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); + amrex::average_node_to_cellcenter( mf_avg, dcomp+1, + *mf_total[1], 0, 1, ngrow); + // Also add the real and imaginary parts of each mode. + /* Not supported yet + for (int i=0 ; i < vector_field[0]->nComp() ; i++) { + MultiFab v_comp0(*vector_field[0], amrex::make_alias, i, 1); + MultiFab v_comp1(*vector_field[1], amrex::make_alias, i, 1); + MultiFab v_comp2(*vector_field[2], amrex::make_alias, i, 1); + srcmf[0] = &v_comp0; + srcmf[1] = &v_comp2; + int id = dcomp + 3*(i + 1); + amrex::average_edge_to_cellcenter( mf_avg, id, srcmf, ngrow); + MultiFab::Copy( mf_avg, mf_avg, id+1, id+2, 1, ngrow); + amrex::average_node_to_cellcenter( mf_avg, id+1, + v_comp1, 0, 1, ngrow); + } + */ +#else + amrex::Abort("AverageAndPackVectorField not implemented for ncomp > 1"); +#endif + } else { + PackPlotDataPtrs(srcmf, vector_field); + amrex::average_edge_to_cellcenter( mf_avg, dcomp, srcmf, ngrow); #if (AMREX_SPACEDIM == 2) - MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); - amrex::average_node_to_cellcenter( mf_avg, dcomp+1, - *vector_field[1], 0, 1, ngrow); + MultiFab::Copy( mf_avg, mf_avg, dcomp+1, dcomp+2, 1, ngrow); + amrex::average_node_to_cellcenter( mf_avg, dcomp+1, + *vector_field[1], 0, 1, ngrow); #endif + } } else { amrex::Abort("Unknown staggering."); @@ -305,6 +398,35 @@ AverageAndPackScalarField( MultiFab& mf_avg, } } +/** \brief Add variable names to the list. + * If there are more that one mode, add the + * name of the total field and then the + * names of the real and imaginary parts of each + * mode. + */ +void +AddToVarNames( Vector<std::string>& varnames, + std::string name, std::string suffix, + int n_rz_azimuthal_modes ) { + auto coords = {"x", "y", "z"}; + for(auto coord:coords) varnames.push_back(name+coord+suffix); + /* Not supported yet + auto coordsRZ = {"r", "theta", "z"}; + if (n_rz_azimuthal_modes > 1) { + // Note that the names are added in the same order as the fields + // are packed in AverageAndPackVectorField + for (int i = 0 ; i < n_rz_azimuthal_modes ; i++) { + for(auto coord:coordsRZ) { + varnames.push_back(name + coord + suffix + std::to_string(i) + "_real"); + } + for(auto coord:coordsRZ) { + varnames.push_back(name + coord + suffix + std::to_string(i) + "_imag"); + } + } + } + */ +} + /** \brief Write the different fields that are meant for output, * into the vector of MultiFab `mf_avg` (one MultiFab per level) * after averaging them to the cell centers. @@ -313,6 +435,14 @@ void WarpX::AverageAndPackFields ( Vector<std::string>& varnames, amrex::Vector<MultiFab>& mf_avg, const int ngrow) const { + // Factor to account for quantities that have multiple components. + // If n_rz_azimuthal_modes > 1, allow space for total field and the real and + // imaginary part of each node. For now, also include the + // imaginary part of mode 0 for code symmetry, even though + // it is always zero. + int modes_factor = 1; + if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1; + // Count how many different fields should be written (ncomp) const int ncomp = fields_to_plot.size() + static_cast<int>(plot_finepatch)*6 @@ -333,17 +463,17 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, // Allocate temp MultiFab with 3 components mf_tmp_E = MultiFab(grids[lev], dmap[lev], 3, ngrow); // Fill MultiFab mf_tmp_E with averaged E - AverageAndPackVectorField(mf_tmp_E, Efield_aux[lev], 0, ngrow); + AverageAndPackVectorField(mf_tmp_E, Efield_aux[lev], dmap[lev], 0, ngrow); } // Same for B if (is_in_vector(fields_to_plot, {"Bx", "By", "Bz"} )){ mf_tmp_B = MultiFab(grids[lev], dmap[lev], 3, ngrow); - AverageAndPackVectorField(mf_tmp_B, Bfield_aux[lev], 0, ngrow); + AverageAndPackVectorField(mf_tmp_B, Bfield_aux[lev], dmap[lev], 0, ngrow); } // Same for J if (is_in_vector(fields_to_plot, {"jx", "jy", "jz"} )){ mf_tmp_J = MultiFab(grids[lev], dmap[lev], 3, ngrow); - AverageAndPackVectorField(mf_tmp_J, current_fp[lev], 0, ngrow); + AverageAndPackVectorField(mf_tmp_J, current_fp[lev], dmap[lev], 0, ngrow); } int dcomp; @@ -437,12 +567,12 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, } if (plot_finepatch) { - AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dcomp, ngrow ); - if(lev==0) for(auto name:{"Ex_fp","Ey_fp","Ez_fp"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dcomp, ngrow ); - if(lev==0) for(auto name:{"Bx_fp","By_fp","Bz_fp"}) varnames.push_back(name); - dcomp += 3; + AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow ); + if (lev == 0) AddToVarNames(varnames, "E", "_fp", n_rz_azimuthal_modes); + dcomp += 3*modes_factor; + AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow ); + if (lev == 0) AddToVarNames(varnames, "B", "_fp", n_rz_azimuthal_modes); + dcomp += 3*modes_factor; } if (plot_crsepatch) @@ -455,11 +585,11 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, { if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch"); std::array<std::unique_ptr<MultiFab>, 3> E = getInterpolatedE(lev); - AverageAndPackVectorField( mf_avg[lev], E, dcomp, ngrow ); + AverageAndPackVectorField( mf_avg[lev], E, dmap[lev], dcomp, ngrow ); } - if(lev==0) for(auto name:{"Ex_cp","Ey_cp","Ez_cp"}) varnames.push_back(name); - dcomp += 3; + if (lev == 0) AddToVarNames(varnames, "E", "_cp", n_rz_azimuthal_modes); + dcomp += 3*modes_factor; // now the magnetic field if (lev == 0) @@ -470,10 +600,10 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, { if (do_nodal) amrex::Abort("TODO: do_nodal && plot_crsepatch"); std::array<std::unique_ptr<MultiFab>, 3> B = getInterpolatedB(lev); - AverageAndPackVectorField( mf_avg[lev], B, dcomp, ngrow ); + AverageAndPackVectorField( mf_avg[lev], B, dmap[lev], dcomp, ngrow ); } - if(lev==0) for(auto name:{"Bx_cp","By_cp","Bz_cp"}) varnames.push_back(name); - dcomp += 3; + if (lev == 0) AddToVarNames(varnames, "B", "_cp", n_rz_azimuthal_modes); + dcomp += 3*modes_factor; } if (costs[0] != nullptr) @@ -535,8 +665,8 @@ WriteRawField( const MultiFab& F, const DistributionMapping& dm, VisMF::Write(F, prefix); } else { // Copy original MultiFab into one that does not have guard cells - MultiFab tmpF( F.boxArray(), dm, 1, 0); - MultiFab::Copy(tmpF, F, 0, 0, 1, 0); + MultiFab tmpF( F.boxArray(), dm, F.nComp(), 0); + MultiFab::Copy(tmpF, F, 0, 0, F.nComp(), 0); VisMF::Write(tmpF, prefix); } @@ -558,7 +688,7 @@ WriteZeroRawField( const MultiFab& F, const DistributionMapping& dm, std::string prefix = amrex::MultiFabFileFullPrefix(lev, filename, level_prefix, field_name); - MultiFab tmpF(F.boxArray(), dm, 1, ng); + MultiFab tmpF(F.boxArray(), dm, F.nComp(), ng); tmpF.setVal(0.); VisMF::Write(tmpF, prefix); } |