diff --git a/.readthedocs.yml b/.readthedocs.yml index 3da9bc77140..95f86fe4ff2 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -9,14 +9,17 @@ version: 2 build: os: ubuntu-22.04 tools: - python: "3.11" + python: "mambaforge-latest" + # python: "3.11" sphinx: - configuration: Docs/source/conf.py + configuration: Docs/source/conf.py -python: - install: - - requirements: Docs/requirements.txt +conda: + environment: Docs/conda.yml +# python: +# install: +# - requirements: Docs/requirements.txt formats: - htmlzip diff --git a/Docs/Doxyfile b/Docs/Doxyfile index 5fbb7651b18..f7740bc0328 100644 --- a/Docs/Doxyfile +++ b/Docs/Doxyfile @@ -2245,7 +2245,7 @@ ENABLE_PREPROCESSING = YES # The default value is: NO. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -MACRO_EXPANSION = NO +MACRO_EXPANSION = YES # If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES then # the macro expansion is limited to the macros specified with the PREDEFINED and @@ -2253,7 +2253,7 @@ MACRO_EXPANSION = NO # The default value is: NO. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -EXPAND_ONLY_PREDEF = NO +EXPAND_ONLY_PREDEF = YES # If the SEARCH_INCLUDES tag is set to YES, the include files in the # INCLUDE_PATH will be searched if a #include is found. @@ -2305,6 +2305,8 @@ PREDEFINED = AMREX_Linux=1 \ WARPX_QED=1 \ WARPX_QED_TABLE_GEN=1 +PREDEFINED += "AMREX_ENUM(CLASS,...)=\"enum class CLASS : int { __VA_ARGS__ };\"" + # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this # tag can be used to specify a list of macro names that should be expanded. The # macro definition that is found in the sources will be used. Use the PREDEFINED @@ -2312,7 +2314,7 @@ PREDEFINED = AMREX_Linux=1 \ # definition found in the source code. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -EXPAND_AS_DEFINED = +EXPAND_AS_DEFINED = AMREX_ENUM # If the SKIP_FUNCTION_MACROS tag is set to YES then doxygen's preprocessor will # remove all references to function-like macros that are alone on a line, have diff --git a/Docs/README.md b/Docs/README.md index e6fac921b04..6d3903ab327 100644 --- a/Docs/README.md +++ b/Docs/README.md @@ -9,12 +9,13 @@ More information can be found in Docs/source/developers/documentation.rst. Install the Python requirements for compiling the documentation: ``` -python3 -m pip install -r Docs/requirements.txt +cd Docs/ +python3 -m pip install -r requirements.txt ``` ### Compiling the documentation -`cd` into the `Docs/` directory and type +Still in the `Docs/` directory, type ``` make html ``` diff --git a/Docs/conda.yml b/Docs/conda.yml new file mode 100644 index 00000000000..1e23c203b2b --- /dev/null +++ b/Docs/conda.yml @@ -0,0 +1,12 @@ +name: readthedocs + +channels: + - conda-forge + - nodefaults + +dependencies: + - python + - doxygen + - pip + - pip: + - -r requirements.txt diff --git a/Docs/requirements.txt b/Docs/requirements.txt index a8c2af0e474..bc34e69cd65 100644 --- a/Docs/requirements.txt +++ b/Docs/requirements.txt @@ -5,7 +5,7 @@ # License: BSD-3-Clause-LBNL # WarpX PICMI bindings w/o C++ component (used for autoclass docs) --e Python +-e ../Python breathe docutils>=0.17.1 diff --git a/Docs/source/developers/documentation.rst b/Docs/source/developers/documentation.rst index a5013299336..5d604bcf9b3 100644 --- a/Docs/source/developers/documentation.rst +++ b/Docs/source/developers/documentation.rst @@ -56,16 +56,15 @@ First, make sure you are in the root directory of WarpX's source and install the .. code-block:: sh - python3 -m pip install -r Docs/requirements.txt + cd Docs/ + python3 -m pip install -r requirements.txt You will also need Doxygen (macOS: ``brew install doxygen``; Ubuntu: ``sudo apt install doxygen``). -Then, to compile the documentation, use +Still in the ``Docs/`` directory, compile the documentation via .. code-block:: sh - cd Docs/ - make html # This will first compile the Doxygen documentation (execute doxygen) # and then build html pages from rst files using sphinx and breathe. diff --git a/Docs/source/developers/fields.rst b/Docs/source/developers/fields.rst index d0af160afef..9d980119814 100644 --- a/Docs/source/developers/fields.rst +++ b/Docs/source/developers/fields.rst @@ -37,6 +37,13 @@ The ``MultiFab`` constructor (for, e.g., ``Ex`` on level ``lev``) is called in ` By default, the ``MultiFab`` are set to ``0`` at initialization. They can be assigned a different value in ``WarpX::InitLevelData``. +Field Names +----------- + +The commonly used WarpX field names are defined in: + +.. doxygenenum:: warpx::fields::FieldType + Field solver ------------ diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst index 47610e0d7ba..275a4dd134d 100644 --- a/Docs/source/usage/workflows/python_extend.rst +++ b/Docs/source/usage/workflows/python_extend.rst @@ -134,9 +134,12 @@ This example accesses the :math:`E_x(x,y,z)` field at level 0 after every time s warpx = sim.extension.warpx # data access - E_x_mf = warpx.multifab(f"Efield_fp[x][level=0]") + # vector field E, component x, on the fine patch of MR level 0 + E_x_mf = warpx.multifab("Efield_fp", dir=0, level=0) + # scalar field rho, on the fine patch of MR level 0 + rho_mf = warpx.multifab("rho_fp", level=0) - # compute + # compute on E_x_mf # iterate over mesh-refinement levels for lev in range(warpx.finest_level + 1): # grow (aka guard/ghost/halo) regions diff --git a/Examples/Physics_applications/spacecraft_charging/inputs_test_rz_spacecraft_charging_picmi.py b/Examples/Physics_applications/spacecraft_charging/inputs_test_rz_spacecraft_charging_picmi.py index e3bc888f600..9ce8bb8433c 100644 --- a/Examples/Physics_applications/spacecraft_charging/inputs_test_rz_spacecraft_charging_picmi.py +++ b/Examples/Physics_applications/spacecraft_charging/inputs_test_rz_spacecraft_charging_picmi.py @@ -121,13 +121,13 @@ def compute_virtual_charge_on_spacecraft(): # Compute integral of rho over volume of the domain # (i.e. total charge of the plasma particles) rho_integral = ( - (rho[1 : nr - 1, 1 : nz - 1] * r[1 : nr - 1, np.newaxis]).sum() * dr * dz + (rho[1 : nr - 1, 1 : nz - 1] * r[1 : nr - 1, np.newaxis]).sum() + * 2 + * np.pi + * dr + * dz ) - # Due to an oddity in WarpX (which will probably be solved later) - # we need to multiply `rho` by `-epsilon_0` to get the correct charge - rho_integral *= 2 * np.pi * -scc.epsilon_0 # does this oddity still exist? - # Compute charge of the spacecraft, based on Gauss theorem q_spacecraft = -rho_integral - scc.epsilon_0 * grad_phi_integral print("Virtual charge on the spacecraft: %e" % q_spacecraft) diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 91d821d6646..f45ca222e69 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -1234,7 +1234,7 @@ PML::CheckPoint ( { using ablastr::fields::Direction; - if (fields.has(FieldType::pml_E_fp, Direction{0}, 0)) + if (fields.has_vector(FieldType::pml_E_fp, 0)) { ablastr::fields::VectorField pml_E_fp = fields.get_alldirs(FieldType::pml_E_fp, 0); ablastr::fields::VectorField pml_B_fp = fields.get_alldirs(FieldType::pml_B_fp, 0); @@ -1246,7 +1246,7 @@ PML::CheckPoint ( VisMF::AsyncWrite(*pml_B_fp[2], dir+"_Bz_fp"); } - if (fields.has(FieldType::pml_E_cp, Direction{0}, 0)) + if (fields.has_vector(FieldType::pml_E_cp, 0)) { ablastr::fields::VectorField pml_E_cp = fields.get_alldirs(FieldType::pml_E_cp, 0); ablastr::fields::VectorField pml_B_cp = fields.get_alldirs(FieldType::pml_B_cp, 0); @@ -1267,7 +1267,7 @@ PML::Restart ( { using ablastr::fields::Direction; - if (fields.has(FieldType::pml_E_fp, Direction{0}, 0)) + if (fields.has_vector(FieldType::pml_E_fp, 0)) { ablastr::fields::VectorField pml_E_fp = fields.get_alldirs(FieldType::pml_E_fp, 0); ablastr::fields::VectorField pml_B_fp = fields.get_alldirs(FieldType::pml_B_fp, 0); @@ -1279,7 +1279,7 @@ PML::Restart ( VisMF::Read(*pml_B_fp[2], dir+"_Bz_fp"); } - if (fields.has(FieldType::pml_E_cp, Direction{0}, 0)) + if (fields.has_vector(FieldType::pml_E_cp, 0)) { ablastr::fields::VectorField pml_E_cp = fields.get_alldirs(FieldType::pml_E_cp, 0); ablastr::fields::VectorField pml_B_cp = fields.get_alldirs(FieldType::pml_B_cp, 0); diff --git a/Source/Diagnostics/BTDiagnostics.H b/Source/Diagnostics/BTDiagnostics.H index d5dd67226b7..d11db98276b 100644 --- a/Source/Diagnostics/BTDiagnostics.H +++ b/Source/Diagnostics/BTDiagnostics.H @@ -241,7 +241,7 @@ private: * will be used by all snapshots to obtain lab-frame data at the respective * z slice location. */ - amrex::Vector > m_cell_centered_data; + std::string const m_cell_centered_data_name; /** Vector of pointers to compute cell-centered data, per level, per component * using the coarsening-ratio provided by the user. */ @@ -346,7 +346,7 @@ private: * \param[in] i_buffer snapshot index */ void SetSnapshotFullStatus (int i_buffer); - /** Vector of field-data stored in the cell-centered multifab, m_cell_centered_data. + /** Vector of field-data stored in the cell-centered MultiFab. * All the fields are stored regardless of the specific fields to plot selected * by the user. */ diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index e00c30aa78e..631de298861 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -56,7 +56,8 @@ namespace } BTDiagnostics::BTDiagnostics (int i, const std::string& name) - : Diagnostics{i, name} + : Diagnostics{i, name}, + m_cell_centered_data_name("BTD_cell_centered_data_" + name) { ReadParameters(); } @@ -83,7 +84,6 @@ void BTDiagnostics::DerivedInitData () m_old_z_boost.resize(m_num_buffers); m_buffer_counter.resize(m_num_buffers); m_snapshot_ncells_lab.resize(m_num_buffers); - m_cell_centered_data.resize(nmax_lev); m_cell_center_functors.resize(nmax_lev); m_max_buffer_multifabs.resize(m_num_buffers); m_buffer_flush_counter.resize(m_num_buffers); @@ -519,7 +519,10 @@ BTDiagnostics::DefineCellCenteredMultiFab(int lev) #else const int ncomps = static_cast(m_cellcenter_varnames.size()); #endif - WarpX::AllocInitMultiFab(m_cell_centered_data[lev], ba, dmap, ncomps, amrex::IntVect(ngrow), lev, "cellcentered_BTD", 0._rt); + bool const remake = false; + bool const redistribute_on_remake = false; + warpx.m_fields.alloc_init(m_cell_centered_data_name, lev, ba, dmap, ncomps, amrex::IntVect(ngrow), 0.0_rt, + remake, redistribute_on_remake); } @@ -540,12 +543,14 @@ BTDiagnostics::InitializeFieldFunctors (int lev) #else auto & warpx = WarpX::GetInstance(); + auto & fields = warpx.m_fields; + // Clear any pre-existing vector to release stored data // This ensures that when domain is load-balanced, the functors point // to the correct field-data pointers m_all_field_functors[lev].clear(); // For back-transformed data, all the components are cell-centered and stored - // in a single multifab, m_cell_centered_data. + // in a single multifab. // Therefore, size of functors at all levels is 1. const int num_BT_functors = 1; m_all_field_functors[lev].resize(num_BT_functors); @@ -554,11 +559,11 @@ BTDiagnostics::InitializeFieldFunctors (int lev) // Create an object of class BackTransformFunctor for (int i = 0; i < num_BT_functors; ++i) { - // coarsening ratio is not provided since the source MultiFab, m_cell_centered_data + // coarsening ratio is not provided since the source MultiFab // is coarsened based on the user-defined m_crse_ratio const int nvars = static_cast(m_varnames.size()); m_all_field_functors[lev][i] = std::make_unique( - m_cell_centered_data[lev].get(), lev, + fields.get(m_cell_centered_data_name, lev), lev, nvars, m_num_buffers, m_varnames, m_varnames_fields); } @@ -570,23 +575,23 @@ BTDiagnostics::InitializeFieldFunctors (int lev) m_cell_center_functors.at(lev).size()); for (int comp=0; comp(warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{0}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "Ey" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Efield_aux, Direction{1}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{1}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "Ez" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Efield_aux, Direction{2}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{2}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "Bx" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{0}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{0}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "By" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{1}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{1}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "Bz" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{2}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{2}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "jx" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp,Direction{0}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp,Direction{0}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "jy" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp,Direction{1}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp,Direction{1}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "jz" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp,Direction{2}, lev), lev, m_crse_ratio); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp,Direction{2}, lev), lev, m_crse_ratio); } else if ( m_cellcenter_varnames[comp] == "rho" ){ m_cell_center_functors[lev][comp] = std::make_unique(lev, m_crse_ratio); } @@ -601,8 +606,9 @@ BTDiagnostics::UpdateVarnamesForRZopenPMD () { #ifdef WARPX_DIM_RZ auto & warpx = WarpX::GetInstance(); + auto & fields = warpx.m_fields; using ablastr::fields::Direction; - const int ncomp_multimodefab = warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, 0)->nComp(); + const int ncomp_multimodefab = fields.get(FieldType::Efield_aux, Direction{0}, 0)->nComp(); const int ncomp = ncomp_multimodefab; @@ -663,21 +669,22 @@ BTDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev) using ablastr::fields::Direction; auto & warpx = WarpX::GetInstance(); - const int ncomp_multimodefab = warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, 0)->nComp(); + auto & fields = warpx.m_fields; + const int ncomp_multimodefab = fields.get(FieldType::Efield_aux, Direction{0}, 0)->nComp(); const int ncomp = ncomp_multimodefab; // Clear any pre-existing vector to release stored data // This ensures that when domain is load-balanced, the functors point // to the correct field-data pointers m_all_field_functors[lev].clear(); // For back-transformed data, all the components are cell-centered and stored - // in a single multifab, m_cell_centered_data. + // in a single MultiFab. // Therefore, size of functors at all levels is 1 const int num_BT_functors = 1; m_all_field_functors[lev].resize(num_BT_functors); for (int i = 0; i < num_BT_functors; ++i) { const int nvars = static_cast(m_varnames.size()); m_all_field_functors[lev][i] = std::make_unique( - m_cell_centered_data[lev].get(), lev, + fields.get(m_cell_centered_data_name, lev), lev, nvars, m_num_buffers, m_varnames, m_varnames_fields); } @@ -689,23 +696,23 @@ BTDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev) const auto m_cell_center_functors_at_lev_size = static_cast(m_cell_center_functors.at(lev).size()); for (int comp=0; comp(warpx.m_fields.get(FieldType::Efield_aux, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "Et" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Efield_aux, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "Ez" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Efield_aux, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Efield_aux, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "Br" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "Bt" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "Bz" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::Bfield_aux, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::Bfield_aux, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "jr" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp, Direction{0}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "jt" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp, Direction{1}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "jz" ){ - m_cell_center_functors[lev][comp] = std::make_unique(warpx.m_fields.get(FieldType::current_fp, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); + m_cell_center_functors[lev][comp] = std::make_unique(fields.get(FieldType::current_fp, Direction{2}, lev), lev, m_crse_ratio, false, ncomp); } else if ( m_cellcenter_varnames_fields[comp] == "rho" ){ m_cell_center_functors[lev][comp] = std::make_unique(lev, m_crse_ratio, false, -1, false, ncomp); } @@ -795,6 +802,8 @@ BTDiagnostics::PrepareFieldDataForOutput () if (!m_do_back_transformed_fields) { return; } auto & warpx = WarpX::GetInstance(); + auto & fields = warpx.m_fields; + // In this function, we will get cell-centered data for every level, lev, // using the cell-center functors and their respective operators() // Call m_cell_center_functors->operator @@ -804,21 +813,23 @@ BTDiagnostics::PrepareFieldDataForOutput () for (int icomp = 0; icompoperator()(*m_cell_centered_data[lev], icomp_dst); + // stores it in cell-centered MultiFab. + m_cell_center_functors[lev][icomp]->operator()(*fields.get(m_cell_centered_data_name, lev), icomp_dst); icomp_dst += m_cell_center_functors[lev][icomp]->nComp(); } // Check that the proper number of user-requested components are cell-centered AMREX_ALWAYS_ASSERT( icomp_dst == m_cellcenter_varnames.size() ); // fill boundary call is required to average_down (flatten) data to // the coarsest level. - ablastr::utils::communication::FillBoundary(*m_cell_centered_data[lev], WarpX::do_single_precision_comms, + ablastr::utils::communication::FillBoundary(*fields.get(m_cell_centered_data_name, lev), + WarpX::do_single_precision_comms, warpx.Geom(lev).periodicity()); } // Flattening out MF over levels for (int lev = warpx.finestLevel(); lev > 0; --lev) { - ablastr::coarsen::sample::Coarsen(*m_cell_centered_data[lev - 1], *m_cell_centered_data[lev], 0, 0, + ablastr::coarsen::sample::Coarsen(*fields.get(m_cell_centered_data_name, lev - 1), + *fields.get(m_cell_centered_data_name, lev), 0, 0, static_cast(m_cellcenter_varnames.size()), 0, WarpX::RefRatio(lev-1) ); } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 93d265d598f..a685afd28e7 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -1147,7 +1147,7 @@ WarpX::PushParticlesandDeposit (int lev, amrex::Real cur_time, DtType a_dt_type, m_fields.get(FieldType::current_fp, Direction{1}, lev), m_fields.get(FieldType::current_fp, Direction{2}, lev), lev); - if (m_fields.has(FieldType::current_buf, Direction{0}, lev)) { + if (m_fields.has_vector(FieldType::current_buf, lev)) { ApplyInverseVolumeScalingToCurrentDensity( m_fields.get(FieldType::current_buf, Direction{0}, lev), m_fields.get(FieldType::current_buf, Direction{1}, lev), diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp index 63b51cb8416..c6a1e206200 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -79,19 +79,19 @@ void FiniteDifferenceSolver::EvolveB ( fields.get(FieldType::G_fp, lev) : fields.get(FieldType::G_cp, lev); } ablastr::fields::VectorField face_areas; - if (fields.has(FieldType::face_areas, Direction{0}, lev)) { + if (fields.has_vector(FieldType::face_areas, lev)) { face_areas = fields.get_alldirs(FieldType::face_areas, lev); } ablastr::fields::VectorField area_mod; - if (fields.has(FieldType::area_mod, Direction{0}, lev)) { + if (fields.has_vector(FieldType::area_mod, lev)) { area_mod = fields.get_alldirs(FieldType::area_mod, lev); } ablastr::fields::VectorField ECTRhofield; - if (fields.has(FieldType::ECTRhofield, Direction{0}, lev)) { + if (fields.has_vector(FieldType::ECTRhofield, lev)) { ECTRhofield = fields.get_alldirs(FieldType::ECTRhofield, lev); } ablastr::fields::VectorField Venl; - if (fields.has(FieldType::Venl, Direction{0}, lev)) { + if (fields.has_vector(FieldType::Venl, lev)) { Venl = fields.get_alldirs(FieldType::Venl, lev); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index db8e80cc972..03a9866fb98 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -73,19 +73,19 @@ void FiniteDifferenceSolver::EvolveE ( } ablastr::fields::VectorField edge_lengths; - if (fields.has(FieldType::edge_lengths, Direction{0}, lev)) { + if (fields.has_vector(FieldType::edge_lengths, lev)) { edge_lengths = fields.get_alldirs(FieldType::edge_lengths, lev); } ablastr::fields::VectorField face_areas; - if (fields.has(FieldType::face_areas, Direction{0}, lev)) { + if (fields.has_vector(FieldType::face_areas, lev)) { face_areas = fields.get_alldirs(FieldType::face_areas, lev); } ablastr::fields::VectorField area_mod; - if (fields.has(FieldType::area_mod, Direction{0}, lev)) { + if (fields.has_vector(FieldType::area_mod, lev)) { area_mod = fields.get_alldirs(FieldType::area_mod, lev); } ablastr::fields::VectorField ECTRhofield; - if (fields.has(FieldType::ECTRhofield, Direction{0}, lev)) { + if (fields.has_vector(FieldType::ECTRhofield, lev)) { ECTRhofield = fields.get_alldirs(FieldType::ECTRhofield, lev); } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp index 9ecae05516d..7a1a05d560d 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp @@ -69,7 +69,7 @@ void FiniteDifferenceSolver::EvolveEPML ( const ablastr::fields::VectorField Jfield = (patch_type == PatchType::fine) ? fields.get_alldirs(FieldType::pml_j_fp, level) : fields.get_alldirs(FieldType::pml_j_cp, level); ablastr::fields::VectorField edge_lengths; - if (fields.has(FieldType::pml_edge_lengths, Direction{0}, level)) { + if (fields.has_vector(FieldType::pml_edge_lengths, level)) { edge_lengths = fields.get_alldirs(FieldType::pml_edge_lengths, level); } amrex::MultiFab * Ffield = nullptr; diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp index 84efe8bf45a..5c28ff1f3c7 100644 --- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp +++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp @@ -170,15 +170,14 @@ WarpX::computeVectorPotential (ablastr::fields::MultiLevelVectorField const& cur amrex::Vector> sorted_curr; amrex::Vector> sorted_A; for (int lev = 0; lev <= finest_level; ++lev) { - sorted_curr.emplace_back(amrex::Array ({curr[lev][Direction{0}], - curr[lev][Direction{1}], - curr[lev][Direction{2}]})); - sorted_A.emplace_back(amrex::Array ({A[lev][Direction{0}], - A[lev][Direction{1}], - A[lev][Direction{2}]})); + sorted_curr.emplace_back(amrex::Array ({curr[lev][0], + curr[lev][1], + curr[lev][2]})); + sorted_A.emplace_back(amrex::Array ({A[lev][0], + A[lev][1], + A[lev][2]})); } -#if defined(AMREX_USE_EB) const ablastr::fields::MultiLevelVectorField Bfield_fp = m_fields.get_mr_levels_alldirs(FieldType::Bfield_fp, finest_level); const std::optional post_A_calculation( { @@ -187,13 +186,13 @@ WarpX::computeVectorPotential (ablastr::fields::MultiLevelVectorField const& cur m_fields.get_mr_levels_alldirs(FieldType::vector_potential_grad_buf_b_stag, finest_level) }); +#if defined(AMREX_USE_EB) amrex::Vector factories; for (int lev = 0; lev <= finest_level; ++lev) { factories.push_back(&WarpX::fieldEBFactory(lev)); } const std::optional > eb_farray_box_factory({factories}); #else - const std::optional post_A_calculation; const std::optional > eb_farray_box_factory; #endif diff --git a/Source/Fields.H b/Source/Fields.H index 8b753bc1008..d470ac33e43 100644 --- a/Source/Fields.H +++ b/Source/Fields.H @@ -18,25 +18,30 @@ namespace warpx::fields { + /** Unique identifiers for WarpX scalar and vector fields. + * + * These are implemented as amrex::MultiFab (one or one per component "direction", + * respectively) and stored in the ablastr::fields::MultiFabRegister . + */ AMREX_ENUM(FieldType, None, - Efield_aux, //!< Field that the particles gather from. Obtained from Efield_fp (and Efield_cp when using MR); see UpdateAuxilaryData - Bfield_aux, //!< Field that the particles gather from. Obtained from Bfield_fp (and Bfield_cp when using MR); see UpdateAuxilaryData - Efield_fp, //!< The field that is updated by the field solver at each timestep - Bfield_fp, //!< The field that is updated by the field solver at each timestep - Efield_fp_external, //!< Stores grid particle fields provided by the user as through an openPMD file - Bfield_fp_external, //!< Stores grid particle fields provided by the user as through an openPMD file - current_fp, //!< The current that is used as a source for the field solver - current_fp_nodal, //!< Only used when using nodal current deposition - current_fp_vay, //!< Only used when using Vay current deposition - current_buf, //!< Particles that are close to the edge of the MR patch (i.e. in the deposition buffer) deposit to this field. - current_store, //!< Only used when doing subcycling with mesh refinement, for book-keeping of currents - rho_buf, //!< Particles that are close to the edge of the MR patch (i.e. in the deposition buffer) deposit to this field. - rho_fp, //!< The charge density that is used as a source for the field solver (mostly for labframe electrostatic and PSATD) - F_fp, //!< Used for divE cleaning - G_fp, //!< Used for divB cleaning - phi_fp, //!< Obtained by the Poisson solver, for labframe electrostatic - vector_potential_fp, //!< Obtained by the magnetostatic solver + Efield_aux, /**< Field that the particles gather from. Obtained from Efield_fp (and Efield_cp when using MR); see UpdateAuxilaryData */ + Bfield_aux, /**< Field that the particles gather from. Obtained from Bfield_fp (and Bfield_cp when using MR); see UpdateAuxilaryData */ + Efield_fp, /**< The field that is updated by the field solver at each timestep */ + Bfield_fp, /**< The field that is updated by the field solver at each timestep */ + Efield_fp_external, /**< Stores grid particle fields provided by the user as through an openPMD file */ + Bfield_fp_external, /**< Stores grid particle fields provided by the user as through an openPMD file */ + current_fp, /**< The current that is used as a source for the field solver */ + current_fp_nodal, /**< Only used when using nodal current deposition */ + current_fp_vay, /**< Only used when using Vay current deposition */ + current_buf, /**< Particles that are close to the edge of the MR patch (i.e. in the deposition buffer) deposit to this field. */ + current_store, /**< Only used when doing subcycling with mesh refinement, for book-keeping of currents */ + rho_buf, /**< Particles that are close to the edge of the MR patch (i.e. in the deposition buffer) deposit to this field. */ + rho_fp, /**< The charge density that is used as a source for the field solver (mostly for labframe electrostatic and PSATD) */ + F_fp, /**< Used for divE cleaning */ + G_fp, /**< Used for divB cleaning */ + phi_fp, /**< Obtained by the Poisson solver, for labframe electrostatic */ + vector_potential_fp, /**< Obtained by the magnetostatic solver */ vector_potential_fp_nodal, vector_potential_grad_buf_e_stag, vector_potential_grad_buf_b_stag, @@ -76,7 +81,7 @@ namespace warpx::fields Bfield_avg_fp, Efield_avg_cp, Bfield_avg_cp, - B_old, //!< Stores the value of B at the beginning of the timestep, for the implicit solver + B_old, /**< Stores the value of B at the beginning of the timestep, for the implicit solver */ ECTRhofield, Venl ); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index ac797d1e706..d64632d964a 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -196,7 +196,7 @@ WarpX::UpdateAuxilaryDataStagToNodal () { if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { Array,3> Btmp; - if (m_fields.has(FieldType::Bfield_cax, Direction{0}, lev)) { + if (m_fields.has_vector(FieldType::Bfield_cax, lev)) { for (int i = 0; i < 3; ++i) { Btmp[i] = std::make_unique( *m_fields.get(FieldType::Bfield_cax, Direction{i}, lev), amrex::make_alias, 0, 1); @@ -290,7 +290,7 @@ WarpX::UpdateAuxilaryDataStagToNodal () { if (electromagnetic_solver_id != ElectromagneticSolverAlgo::None) { Array,3> Etmp; - if (m_fields.has(FieldType::Efield_cax, Direction{0}, lev)) { + if (m_fields.has_vector(FieldType::Efield_cax, lev)) { for (int i = 0; i < 3; ++i) { Etmp[i] = std::make_unique( *m_fields.get(FieldType::Efield_cax, Direction{i}, lev), amrex::make_alias, 0, 1); @@ -450,7 +450,7 @@ WarpX::UpdateAuxilaryDataSameType () Bfield_aux[lev - 1][2]->nComp(), ng_src, ng, WarpX::do_single_precision_comms, crse_period); - if (m_fields.has(FieldType::Bfield_cax, Direction{0}, lev)) + if (m_fields.has_vector(FieldType::Bfield_cax, lev)) { MultiFab::Copy(*m_fields.get(FieldType::Bfield_cax, Direction{0}, lev), dBx, 0, 0, m_fields.get(FieldType::Bfield_cax, Direction{0}, lev)->nComp(), ng); MultiFab::Copy(*m_fields.get(FieldType::Bfield_cax, Direction{1}, lev), dBy, 0, 0, m_fields.get(FieldType::Bfield_cax, Direction{1}, lev)->nComp(), ng); @@ -535,7 +535,7 @@ WarpX::UpdateAuxilaryDataSameType () WarpX::do_single_precision_comms, crse_period); - if (m_fields.has(FieldType::Efield_cax, Direction{0}, lev)) + if (m_fields.has_vector(FieldType::Efield_cax, lev)) { MultiFab::Copy(*m_fields.get(FieldType::Efield_cax, Direction{0}, lev), dEx, 0, 0, m_fields.get(FieldType::Efield_cax, Direction{0}, lev)->nComp(), ng); MultiFab::Copy(*m_fields.get(FieldType::Efield_cax, Direction{1}, lev), dEy, 0, 0, m_fields.get(FieldType::Efield_cax, Direction{1}, lev)->nComp(), ng); diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 10849a0c0d5..c804bb12797 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -586,7 +586,7 @@ LaserParticleContainer::Evolve (ablastr::fields::MultiFabRegister& fields, amrex::LayoutData* cost = WarpX::getCosts(lev); const bool has_rho = fields.has(FieldType::rho_fp, lev); - const bool has_buffer = fields.has(FieldType::current_buf, lev); + const bool has_buffer = fields.has_vector(FieldType::current_buf, lev); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 07997a61f0c..26f9fee38d3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1753,9 +1753,9 @@ PhysicalParticleContainer::Evolve (ablastr::fields::MultiFabRegister& fields, const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); const bool has_rho = fields.has(FieldType::rho_fp, lev); - const bool has_cjx = fields.has(FieldType::current_buf, Direction{0}, lev); - const bool has_cEx = fields.has(FieldType::Efield_cax, Direction{0}, lev); - const bool has_buffer = has_cEx || has_cjx; + const bool has_J_buf = fields.has_vector(FieldType::current_buf, lev); + const bool has_E_cax = fields.has_vector(FieldType::Efield_cax, lev); + const bool has_buffer = has_E_cax || has_J_buf; amrex::MultiFab & Ex = *fields.get(FieldType::Efield_aux, Direction{0}, lev); amrex::MultiFab & Ey = *fields.get(FieldType::Efield_aux, Direction{1}, lev); @@ -1850,7 +1850,7 @@ PhysicalParticleContainer::Evolve (ablastr::fields::MultiFabRegister& fields, pti, lev, current_masks, gather_masks ); } - const long np_current = has_cjx ? nfine_current : np; + const long np_current = has_J_buf ? nfine_current : np; if (has_rho && ! skip_deposition && ! do_not_deposit) { // Deposit charge before particle push, in component 0 of MultiFab rho. @@ -1870,7 +1870,7 @@ PhysicalParticleContainer::Evolve (ablastr::fields::MultiFabRegister& fields, if (! do_not_push) { - const long np_gather = has_cEx ? nfine_gather : np; + const long np_gather = has_E_cax ? nfine_gather : np; int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp index 857d23dc588..0aab95f78f8 100644 --- a/Source/Python/WarpX.cpp +++ b/Source/Python/WarpX.cpp @@ -120,17 +120,35 @@ void init_WarpX (py::module& m) R"doc(Registry to all WarpX MultiFab (fields).)doc" ) .def("multifab", - [](WarpX & wx, std::string multifab_name, int level) { - if (wx.m_fields.has(multifab_name, level)) { - return wx.m_fields.get(multifab_name, level); + [](WarpX & wx, std::string internal_name) { + if (wx.m_fields.internal_has(internal_name)) { + return wx.m_fields.internal_get(internal_name); + } else { + throw std::runtime_error("MultiFab '" + internal_name + "' is unknown or is not allocated!"); + } + }, + py::arg("internal_name"), + py::return_value_policy::reference_internal, + R"doc(Return a MultiFab by its internal name (deprecated). + +The multifab('internal_name') signature is deprecated. +Please use: +- multifab('prefix', level=...) for scalar fields +- multifab('prefix', dir=..., level=...) for vector field components +where 'prefix' is the part of 'internal_name';' before the [])doc" + ) + .def("multifab", + [](WarpX & wx, std::string scalar_name, int level) { + if (wx.m_fields.has(scalar_name, level)) { + return wx.m_fields.get(scalar_name, level); } else { - throw std::runtime_error("The MultiFab '" + multifab_name + "' is unknown or is not allocated!"); + throw std::runtime_error("The scalar field '" + scalar_name + "' is unknown or is not allocated!"); } }, - py::arg("multifab_name"), + py::arg("scalar_name"), py::arg("level"), py::return_value_policy::reference_internal, - R"doc(Return MultiFabs by name and level, e.g., ``\"Efield_aux\"``, ``\"Efield_fp"``, ... + R"doc(Return scalar fields (MultiFabs) by name and level, e.g., ``\"rho_fp\"``, ``\"phi_fp"``, ... The physical fields in WarpX have the following naming: @@ -141,18 +159,18 @@ The physical fields in WarpX have the following naming: (only for level 1 and higher).)doc" ) .def("multifab", - [](WarpX & wx, std::string multifab_name, Direction dir, int level) { - if (wx.m_fields.has(multifab_name, dir, level)) { - return wx.m_fields.get(multifab_name, dir, level); + [](WarpX & wx, std::string vector_name, Direction dir, int level) { + if (wx.m_fields.has(vector_name, dir, level)) { + return wx.m_fields.get(vector_name, dir, level); } else { - throw std::runtime_error("The MultiFab '" + multifab_name + "' is unknown or is not allocated!"); + throw std::runtime_error("The vector field '" + vector_name + "' is unknown or is not allocated!"); } }, - py::arg("multifab_name"), + py::arg("vector_name"), py::arg("dir"), py::arg("level"), py::return_value_policy::reference_internal, - R"doc(Return MultiFabs by name, direction, and level, e.g., ``\"Efield_aux\"``, ``\"Efield_fp"``, ... + R"doc(Return the component of a vector field (MultiFab) by name, direction, and level, e.g., ``\"Efield_aux\"``, ``\"Efield_fp"``, ... The physical fields in WarpX have the following naming: diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 88b6e1d28c1..91153f4ed8e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1388,8 +1388,9 @@ WarpX::ReadParameters () // Instead, if warpx.grid_type=collocated, the momentum-conserving and // energy conserving field gathering algorithms are equivalent (forces // gathered from the collocated grid) and no fields centering occurs. - if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving && - WarpX::grid_type != GridType::Collocated) + if ((WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving + && WarpX::grid_type != GridType::Collocated) + || WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) { utils::parser::queryWithParser( pp_warpx, "field_centering_nox", field_centering_nox); diff --git a/Source/ablastr/fields/MultiFabRegister.H b/Source/ablastr/fields/MultiFabRegister.H index f33eed1c5a6..21df20c1678 100644 --- a/Source/ablastr/fields/MultiFabRegister.H +++ b/Source/ablastr/fields/MultiFabRegister.H @@ -661,14 +661,21 @@ namespace ablastr::fields int level ) const; - private: + /** Temporary test function for legacy Python bindings */ + [[nodiscard]] bool + internal_has ( + std::string const & internal_name + ); [[nodiscard]] amrex::MultiFab * internal_get ( - std::string const & key + std::string const & internal_name ); + + private: + [[nodiscard]] amrex::MultiFab const * internal_get ( - std::string const & key + std::string const & internal_name ) const; amrex::MultiFab* diff --git a/Source/ablastr/fields/MultiFabRegister.cpp b/Source/ablastr/fields/MultiFabRegister.cpp index 106a3aede79..2c384a90089 100644 --- a/Source/ablastr/fields/MultiFabRegister.cpp +++ b/Source/ablastr/fields/MultiFabRegister.cpp @@ -336,32 +336,40 @@ namespace ablastr::fields return count == 3; } + bool + MultiFabRegister::internal_has ( + std::string const & internal_name + ) + { + return m_mf_register.count(internal_name) > 0; + } + amrex::MultiFab* MultiFabRegister::internal_get ( - std::string const & key + std::string const & internal_name ) { - if (m_mf_register.count(key) == 0) { + if (m_mf_register.count(internal_name) == 0) { // FIXME: temporary, throw a std::runtime_error // throw std::runtime_error("MultiFabRegister::get name does not exist in register: " + key); return nullptr; } - amrex::MultiFab & mf = m_mf_register.at(key).m_mf; + amrex::MultiFab & mf = m_mf_register.at(internal_name).m_mf; return &mf; } amrex::MultiFab const * MultiFabRegister::internal_get ( - std::string const & key + std::string const & internal_name ) const { - if (m_mf_register.count(key) == 0) { + if (m_mf_register.count(internal_name) == 0) { // FIXME: temporary, throw a std::runtime_error - // throw std::runtime_error("MultiFabRegister::get name does not exist in register: " + key); + // throw std::runtime_error("MultiFabRegister::get name does not exist in register: " + internal_name); return nullptr; } - amrex::MultiFab const & mf = m_mf_register.at(key).m_mf; + amrex::MultiFab const & mf = m_mf_register.at(internal_name).m_mf; return &mf; } diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H index 727280d630b..8b4f9cea9a1 100644 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -14,6 +14,7 @@ #include #include #include +#include #include #if defined(ABLASTR_USE_FFT) && defined(WARPX_DIM_3D) @@ -53,6 +54,7 @@ #include #include +#include namespace ablastr::fields { @@ -275,7 +277,7 @@ computePhi ( #endif // Use the Multigrid (MLMG) solver if selected or on refined patches // but first scale rho appropriately - rho[lev]->mult(-1._rt / ablastr::constant::SI::ep0); // TODO: when do we "un-multiply" this? We need to document this side-effect! + rho[lev]->mult(-1._rt / ablastr::constant::SI::ep0); #ifdef WARPX_DIM_RZ constexpr bool is_rz = true; @@ -310,13 +312,17 @@ computePhi ( auto linop_nodelap = std::make_unique(); if (eb_enabled) { #if defined(AMREX_USE_EB) - linop_nodelap->define( - amrex::Vector{geom[lev]}, - amrex::Vector{grids[lev]}, - amrex::Vector{dmap[lev]}, - info, - amrex::Vector{eb_farray_box_factory.value()[lev]} - ); + if constexpr(std::is_same_v) { + throw std::runtime_error("EB requested by eb_farray_box_factory not provided!"); + } else { + linop_nodelap->define( + amrex::Vector{geom[lev]}, + amrex::Vector{grids[lev]}, + amrex::Vector{dmap[lev]}, + info, + amrex::Vector{eb_farray_box_factory.value()[lev]} + ); + } #endif } else { @@ -409,6 +415,8 @@ computePhi ( post_phi_calculation.value()(mlmg, lev); } } + rho[lev]->mult(-ablastr::constant::SI::ep0); // Multiply rho by epsilon again + } // loop over lev(els) } // computePhi } // namespace ablastr::fields