diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Diagnostics/FieldIO.cpp | 200 | ||||
| -rw-r--r-- | Source/Utils/write_atomic_data_cpp.py | 1 | ||||
| -rw-r--r-- | Source/WarpX.H | 31 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 20 |
4 files changed, 209 insertions, 43 deletions
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index e1bb8cb54..43f204da0 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -352,6 +352,31 @@ AverageAndPackVectorField( MultiFab& mf_avg, } } +/** \brief Takes the specified component of the three fields and + * averages and packs them into the MultiFab mf_avg. + */ +void +AverageAndPackVectorFieldComponents (MultiFab& mf_avg, + const std::array< std::unique_ptr<MultiFab>, 3 >& vector_field, + const DistributionMapping& dm, + int& dcomp, const int ngrow ) +{ + if (vector_field[0]->nComp() > 1) { + std::array<std::unique_ptr<MultiFab>,3> vector_field_component; + vector_field_component[0].reset(new MultiFab(vector_field[0]->boxArray(), dm, 1, vector_field[0]->nGrowVect())); + vector_field_component[1].reset(new MultiFab(vector_field[1]->boxArray(), dm, 1, vector_field[1]->nGrowVect())); + vector_field_component[2].reset(new MultiFab(vector_field[2]->boxArray(), dm, 1, vector_field[2]->nGrowVect())); + for (int icomp=0 ; icomp < vector_field[0]->nComp() ; icomp++) { + // This is a bit inefficient since it copies the data. Can this be done with pointers? + MultiFab::Copy(*vector_field_component[0], *vector_field[0], icomp, 0, 1, vector_field[0]->nGrowVect()); + MultiFab::Copy(*vector_field_component[1], *vector_field[1], icomp, 0, 1, vector_field[1]->nGrowVect()); + MultiFab::Copy(*vector_field_component[2], *vector_field[2], icomp, 0, 1, vector_field[2]->nGrowVect()); + AverageAndPackVectorField(mf_avg, vector_field_component, dm, dcomp, ngrow); + dcomp += 3; + } + } +} + /** \brief Take a MultiFab `scalar_field` * averages it to the cell center, and stores the * resulting MultiFab in mf_avg (in the components dcomp) @@ -378,6 +403,74 @@ AverageAndPackScalarField( MultiFab& mf_avg, } } +/** \brief Takes the specified component of the scalar and + * averages and packs it into the MultiFab mf_avg. + */ +void +AverageAndPackScalarFieldComponent (MultiFab& mf_avg, + const MultiFab& scalar_field, + const int icomp, + const int dcomp, const int ngrow ) +{ + MultiFab scalar_field_component(scalar_field, amrex::make_alias, icomp, 1); + AverageAndPackScalarField(mf_avg, scalar_field_component, dcomp, ngrow); +} + +/** \brief Generate mode variable name + */ +std::string +ComponentName(std::string fieldname, int mode, std::string type) +{ + if (type == "real") { + return fieldname + "_" + std::to_string(mode) + "_" + "real"; + } else if (type == "imag") { + return fieldname + "_" + std::to_string(mode) + "_" + "imag"; + } else { + AMREX_ALWAYS_ASSERT( false ); + } + // This should never be done + return ""; +} + +/* \brief Copy vector field component data into the MultiFab that will be written out + */ +void +CopyVectorFieldComponentsToMultiFab (int lev, amrex::Vector<MultiFab>& mf_avg, MultiFab& mf_tmp, + int icomp, int& dcomp, int ngrow, + std::string fieldname, Vector<std::string>& varnames) +{ + if (mf_tmp.nComp() > 3) { + if (lev==0) varnames.push_back(ComponentName(fieldname, 0, "real")); + MultiFab::Copy( mf_avg[lev], mf_tmp, 3+icomp, dcomp++, 1, ngrow); + int const nmodes = mf_tmp.nComp()/6; + for (int mode=1 ; mode < nmodes ; mode++) { + if (lev==0) varnames.push_back(ComponentName(fieldname, mode, "real")); + MultiFab::Copy( mf_avg[lev], mf_tmp, 3*2*mode+icomp, dcomp++, 1, ngrow); + if (lev==0) varnames.push_back(ComponentName(fieldname, mode, "imag")); + MultiFab::Copy( mf_avg[lev], mf_tmp, 3*2*mode+3+icomp, dcomp++, 1, ngrow); + } + } +} + +/* \brief Copy scalar field component data into the MultiFab that will be written out + */ +void +CopyScalarFieldComponentsToMultiFab (int lev, amrex::Vector<MultiFab>& mf_avg, MultiFab& mf_tmp, + int& dcomp, int ngrow, int n_rz_azimuthal_modes, + std::string fieldname, Vector<std::string>& varnames) +{ + if (n_rz_azimuthal_modes > 1) { + if (lev==0) varnames.push_back(ComponentName(fieldname, 0, "real")); + AverageAndPackScalarFieldComponent(mf_avg[lev], mf_tmp, 0, dcomp++, ngrow); + for (int mode=1 ; mode < n_rz_azimuthal_modes ; mode++) { + if (lev==0) varnames.push_back(ComponentName(fieldname, mode, "real")); + AverageAndPackScalarFieldComponent(mf_avg[lev], mf_tmp, 2*mode-1, dcomp++, ngrow); + if (lev==0) varnames.push_back(ComponentName(fieldname, mode, "imag")); + AverageAndPackScalarFieldComponent(mf_avg[lev], mf_tmp, 2*mode , dcomp++, ngrow); + } + } +} + /** \brief Add variable names to the list. */ void @@ -396,11 +489,28 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, amrex::Vector<MultiFab>& mf_avg, const int ngrow) const { // Count how many different fields should be written (ncomp) - const int ncomp = fields_to_plot.size() + int ncomp = fields_to_plot.size() + static_cast<int>(plot_finepatch)*6 + static_cast<int>(plot_crsepatch)*6 + static_cast<int>(costs[0] != nullptr and plot_costs); + // Add in the RZ modes + if (n_rz_azimuthal_modes > 1) { + for (std::string field : fields_to_plot) { + if (is_in_vector({"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz", "rho", "F"}, {field})) { + ncomp += 2*n_rz_azimuthal_modes - 1; + } + } + if (plot_finepatch) { + ncomp += 6*(2*n_rz_azimuthal_modes - 1); + } + } + + int nvecs = 3; + if (n_rz_azimuthal_modes > 1) { + nvecs += 3*(2*n_rz_azimuthal_modes - 1); + } + // Loop over levels of refinement for (int lev = 0; lev <= finest_level; ++lev) { @@ -413,19 +523,25 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, // Build mf_tmp_E is at least one component of E is requested if (is_in_vector(fields_to_plot, {"Ex", "Ey", "Ez"} )){ // Allocate temp MultiFab with 3 components - mf_tmp_E = MultiFab(grids[lev], dmap[lev], 3, ngrow); + mf_tmp_E = MultiFab(grids[lev], dmap[lev], nvecs, ngrow); // Fill MultiFab mf_tmp_E with averaged E AverageAndPackVectorField(mf_tmp_E, Efield_aux[lev], dmap[lev], 0, ngrow); + int dcomp = 3; + AverageAndPackVectorFieldComponents(mf_tmp_E, Efield_aux[lev], dmap[lev], dcomp, ngrow); } // Same for B if (is_in_vector(fields_to_plot, {"Bx", "By", "Bz"} )){ - mf_tmp_B = MultiFab(grids[lev], dmap[lev], 3, ngrow); + mf_tmp_B = MultiFab(grids[lev], dmap[lev], nvecs, ngrow); AverageAndPackVectorField(mf_tmp_B, Bfield_aux[lev], dmap[lev], 0, ngrow); + int dcomp = 3; + AverageAndPackVectorFieldComponents(mf_tmp_B, Bfield_aux[lev], dmap[lev], dcomp, ngrow); } // Same for J if (is_in_vector(fields_to_plot, {"jx", "jy", "jz"} )){ - mf_tmp_J = MultiFab(grids[lev], dmap[lev], 3, ngrow); + mf_tmp_J = MultiFab(grids[lev], dmap[lev], nvecs, ngrow); AverageAndPackVectorField(mf_tmp_J, current_fp[lev], dmap[lev], 0, ngrow); + int dcomp = 3; + AverageAndPackVectorFieldComponents(mf_tmp_J, current_fp[lev], dmap[lev], dcomp, ngrow); } int dcomp; @@ -433,37 +549,51 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, // mf_avg[lev] add the corresponding names to `varnames`. // plot_fine_patch and plot_coarse_patch are treated separately // (after this for loop). - for (dcomp=0; dcomp<fields_to_plot.size(); dcomp++){ - std::string fieldname = fields_to_plot[dcomp]; + dcomp = 0; + for (int ifield=0; ifield<fields_to_plot.size(); ifield++){ + std::string fieldname = fields_to_plot[ifield]; if(lev==0) varnames.push_back(fieldname); if (fieldname == "Ex"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_E, 0, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_E, 0, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_E, 0, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "Ey"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_E, 1, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_E, 1, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_E, 1, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "Ez"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_E, 2, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_E, 2, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_E, 2, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "Bx"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_B, 0, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_B, 0, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_B, 0, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "By"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_B, 1, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_B, 1, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_B, 1, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "Bz"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_B, 2, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_B, 2, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_B, 2, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "jx"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_J, 0, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_J, 0, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_J, 0, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "jy"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_J, 1, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_J, 1, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_J, 1, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "jz"){ - MultiFab::Copy( mf_avg[lev], mf_tmp_J, 2, dcomp, 1, ngrow); + MultiFab::Copy( mf_avg[lev], mf_tmp_J, 2, dcomp++, 1, ngrow); + CopyVectorFieldComponentsToMultiFab(lev, mf_avg, mf_tmp_J, 2, dcomp, ngrow, fieldname, varnames); } else if (fieldname == "rho"){ - AverageAndPackScalarField( mf_avg[lev], *rho_fp[lev], dcomp, ngrow ); + AverageAndPackScalarField( mf_avg[lev], *rho_fp[lev], dcomp++, ngrow ); + CopyScalarFieldComponentsToMultiFab(lev, mf_avg, *rho_fp[lev], dcomp, ngrow, n_rz_azimuthal_modes, + fieldname, varnames); } else if (fieldname == "F"){ - AverageAndPackScalarField( mf_avg[lev], *F_fp[lev], dcomp, ngrow); + AverageAndPackScalarField( mf_avg[lev], *F_fp[lev], dcomp++, ngrow ); + CopyScalarFieldComponentsToMultiFab(lev, mf_avg, *F_fp[lev], dcomp, ngrow, n_rz_azimuthal_modes, + fieldname, varnames); } else if (fieldname == "part_per_cell") { MultiFab temp_dat(grids[lev],mf_avg[lev].DistributionMap(),1,0); temp_dat.setVal(0); // MultiFab containing number of particles in each cell mypc->Increment(temp_dat, lev); - AverageAndPackScalarField( mf_avg[lev], temp_dat, dcomp, ngrow ); + AverageAndPackScalarField( mf_avg[lev], temp_dat, dcomp++, ngrow ); } else if (fieldname == "part_per_grid"){ const Vector<long>& npart_in_grid = mypc->NumberOfParticlesInGrid(lev); // MultiFab containing number of particles per grid @@ -473,7 +603,7 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, #endif for (MFIter mfi(mf_avg[lev]); mfi.isValid(); ++mfi) { (mf_avg[lev])[mfi].setVal(static_cast<Real>(npart_in_grid[mfi.index()]), - dcomp); + dcomp++); } } else if (fieldname == "part_per_proc"){ const Vector<long>& npart_in_grid = mypc->NumberOfParticlesInGrid(lev); @@ -486,7 +616,7 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, for (MFIter mfi(mf_avg[lev]); mfi.isValid(); ++mfi) { n_per_proc += npart_in_grid[mfi.index()]; } - mf_avg[lev].setVal(static_cast<Real>(n_per_proc), dcomp,1); + mf_avg[lev].setVal(static_cast<Real>(n_per_proc), dcomp++,1); } else if (fieldname == "proc_number"){ // MultiFab containing the Processor ID #ifdef _OPENMP @@ -494,11 +624,11 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, #endif for (MFIter mfi(mf_avg[lev]); mfi.isValid(); ++mfi) { (mf_avg[lev])[mfi].setVal(static_cast<Real>(ParallelDescriptor::MyProc()), - dcomp); + dcomp++); } } else if (fieldname == "divB"){ if (do_nodal) amrex::Abort("TODO: do_nodal && plot divb"); - ComputeDivB(mf_avg[lev], dcomp, + ComputeDivB(mf_avg[lev], dcomp++, {Bfield_aux[lev][0].get(), Bfield_aux[lev][1].get(), Bfield_aux[lev][2].get()}, @@ -512,7 +642,7 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, Efield_aux[lev][1].get(), Efield_aux[lev][2].get()}, WarpX::CellSize(lev) ); - AverageAndPackScalarField( mf_avg[lev], dive, dcomp, ngrow ); + AverageAndPackScalarField( mf_avg[lev], dive, dcomp++, ngrow ); } else { amrex::Abort("unknown field in fields_to_plot: " + fieldname); } @@ -520,11 +650,31 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, if (plot_finepatch) { AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow ); - if (lev == 0) AddToVarNames(varnames, "E", "_fp"); dcomp += 3; + AverageAndPackVectorFieldComponents(mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow); + if (lev == 0) { + AddToVarNames(varnames, "E", "_fp"); + if (n_rz_azimuthal_modes > 1) { + AddToVarNames(varnames, "E", ComponentName("_fp", 0, "real")); + for (int mode=1 ; mode < n_rz_azimuthal_modes ; mode++) { + AddToVarNames(varnames, "E", ComponentName("_fp", mode, "real")); + AddToVarNames(varnames, "E", ComponentName("_fp", mode, "imag")); + } + } + } AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow ); - if (lev == 0) AddToVarNames(varnames, "B", "_fp"); dcomp += 3; + AverageAndPackVectorFieldComponents(mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow); + if (lev == 0) { + AddToVarNames(varnames, "B", "_fp"); + if (n_rz_azimuthal_modes > 1) { + AddToVarNames(varnames, "B", ComponentName("_fp", 0, "real")); + for (int mode=1 ; mode < n_rz_azimuthal_modes ; mode++) { + AddToVarNames(varnames, "B", ComponentName("_fp", mode, "real")); + AddToVarNames(varnames, "B", ComponentName("_fp", mode, "imag")); + } + } + } } if (plot_crsepatch) diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py index 35baf2cdf..09b7b8300 100644 --- a/Source/Utils/write_atomic_data_cpp.py +++ b/Source/Utils/write_atomic_data_cpp.py @@ -6,7 +6,6 @@ IonizationEnergiesTable.H, which contains tables + metadata. import re, os import numpy as np -from scipy.constants import e filename = os.path.join( '.', 'atomic_data.txt' ) with open(filename) as f: diff --git a/Source/WarpX.H b/Source/WarpX.H index bafba1b6c..b9f969abb 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -354,6 +354,21 @@ public: const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi, std::array<amrex::Real, 3> const beta = {{0,0,0}} ) const; + /** + * \brief + * This function initializes the E and B fields on each level + * using the parser and the user-defined function for the external fields. + * The subroutine will parse the x_/y_z_external_grid_function and + * then, the B or E multifab is initialized based on the (x,y,z) position + * on the staggered yee-grid or cell-centered grid. + */ + void InitializeExternalFieldsOnGridUsingParser ( + amrex::MultiFab *mfx, amrex::MultiFab *mfy, amrex::MultiFab *mfz, + ParserWrapper *xfield_parser, ParserWrapper *yfield_parser, + ParserWrapper *zfield_parser, amrex::IntVect x_nodal_flag, + amrex::IntVect y_nodal_flag, amrex::IntVect z_nodal_flag, + const int lev); + protected: /** @@ -383,22 +398,6 @@ protected: */ void InitLevelData (int lev, amrex::Real time); - /** - * \brief - * This function initializes the E and B fields on each level - * using the parser and the user-defined function for the external fields. - * The subroutine will parse the x_/y_z_external_grid_function and - * then, the B or E multifab is initialized based on the (x,y,z) position - * on the staggered yee-grid or cell-centered grid. - */ - void InitializeExternalFieldsOnGridUsingParser ( - amrex::MultiFab *mfx, amrex::MultiFab *mfy, amrex::MultiFab *mfz, - ParserWrapper *xfield_parser, ParserWrapper *yfield_parser, - ParserWrapper *zfield_parser, amrex::IntVect x_nodal_flag, - amrex::IntVect y_nodal_flag, amrex::IntVect z_nodal_flag, - const int lev); - - //! Tagging cells for refinement virtual void ErrorEst (int lev, amrex::TagBoxArray& tags, amrex::Real time, int /*ngrow*/) final; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e2be8849f..7c5fe1834 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -260,6 +260,25 @@ WarpX::WarpX () // at different levels (the stencil depends on c*dt/dz) nci_godfrey_filter_exeybz.resize(nlevs_max); nci_godfrey_filter_bxbyez.resize(nlevs_max); + + // Sanity checks. Must be done after calling the MultiParticleContainer + // constructor, as it reads additional parameters + // (e.g., use_fdtd_nci_corr) + +#ifndef WARPX_USE_PSATD + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + not ( do_pml && do_nodal ), + "PML + do_nodal for finite-difference not implemented" + ); +#endif + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + not ( do_dive_cleaning && do_nodal ), + "divE cleaning + do_nodal not implemented" + ); +#ifdef WARPX_USE_PSATD + AMREX_ALWAYS_ASSERT(use_fdtd_nci_corr == 0); + AMREX_ALWAYS_ASSERT(do_subcycling == 0); +#endif } WarpX::~WarpX () @@ -644,7 +663,6 @@ WarpX::ReadParameters () } } - } // This is a virtual function. |