diff options
Diffstat (limited to 'Source/Initialization/WarpXInitData.cpp')
| -rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 356 |
1 files changed, 349 insertions, 7 deletions
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 29c9a8955..2166adbe2 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -10,6 +10,8 @@ #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> #endif +#include <GpuParser.H> + using namespace amrex; @@ -231,23 +233,228 @@ WarpX::PostRestart () mypc->PostRestart(); } + +namespace { +WarpXParser makeParser (std::string const& parse_function) +{ + WarpXParser parser(parse_function); + parser.registerVariables({"x","y","z"}); + ParmParse pp("my_constants"); + std::set<std::string> symbols = parser.symbols(); + symbols.erase("x"); + symbols.erase("y"); + symbols.erase("z"); + for (auto it = symbols.begin(); it != symbols.end(); ) { + Real v; + if (pp.query(it->c_str(), v)) { + parser.setConstant(*it, v); + it = symbols.erase(it); + } else { + ++it; + } + } + for (auto const& s : symbols) { + amrex::Abort(" ExternalEBFieldOnGrid::makeParser::Unknown symbol "+s); + } + return parser; +} +} + + + + +/* \brief + * This function initializes E, B, rho, and F, at all the levels + * of the multifab. rho and F are initialized with 0. + * The E and B fields are initialized using user-defined inputs. + * The initialization type is set using "B_ext_grid_init_style" + * and "E_ext_grid_init_style". The initialization style is set to "default" + * if not explicitly defined by the user, and the E and B fields are + * initialized with E_external_grid and B_external_grid, respectively, each with + * a default value of 0. + * If the initialization type for the E and B field is "constant", + * then, the E and B fields at all the levels are initialized with + * user-defined values for E_external_grid and B_external_grid. + * If the initialization type for B-field is set to + * "parse_B_ext_grid_function", then, the parser is used to read + * Bx_external_grid_function(x,y,z), By_external_grid_function(x,y,z), + * and Bz_external_grid_function(x,y,z). + * Similarly, if the E-field initialization type is set to + * "parse_E_ext_grid_function", then, the parser is used to read + * Ex_external_grid_function(x,y,z), Ey_external_grid_function(x,y,z), + * and Ex_external_grid_function(x,y,z). The parser for the E and B + * initialization assumes that the function has three independent + * variables, at max, namely, x, y, z. However, any number of constants + * can be used in the function used to define the E and B fields on the grid. + */ + void WarpX::InitLevelData (int lev, Real time) { + + ParmParse pp("warpx"); + + // default values of E_external_grid and B_external_grid + // are used to set the E and B field when "constant" or + // "parser" is not explicitly used in the input. + pp.query("B_ext_grid_init_style", B_ext_grid_s); + std::transform(B_ext_grid_s.begin(), + B_ext_grid_s.end(), + B_ext_grid_s.begin(), + ::tolower); + + pp.query("E_ext_grid_init_style", E_ext_grid_s); + std::transform(E_ext_grid_s.begin(), + E_ext_grid_s.end(), + E_ext_grid_s.begin(), + ::tolower); + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (B_ext_grid_s == "constant") + pp.getarr("B_external_grid", B_external_grid); + + // if the input string is "constant", the values for the + // external grid must be provided in the input. + if (E_ext_grid_s == "constant") + pp.getarr("E_external_grid", E_external_grid); + for (int i = 0; i < 3; ++i) { current_fp[lev][i]->setVal(0.0); - Efield_fp[lev][i]->setVal(E_external_grid[i]); - Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") + Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") + Efield_fp[lev][i]->setVal(E_external_grid[i]); + } + + if (B_ext_grid_s == "parse_b_ext_grid_function") { + + std::vector<std::string> f; + // Parse Bx_external_grid_function + pp.getarr("Bx_external_grid_function(x,y,z)", f); + str_Bx_ext_grid_function.clear(); + for (auto const& s : f) { + str_Bx_ext_grid_function += s; + } + f.clear(); + + // Parse By_external_grid_function + pp.getarr("By_external_grid_function(x,y,z)", f); + str_By_ext_grid_function.clear(); + for (auto const& s : f) { + str_By_ext_grid_function += s; + } + f.clear(); + + // Parse Bz_external_grid_function + pp.getarr("Bz_external_grid_function(x,y,z)", f); + str_Bz_ext_grid_function.clear(); + for (auto const& s : f) { + str_Bz_ext_grid_function += s; + } + f.clear(); + + // Initialize Bfield_fp with external function + MultiFab *Bx, *By, *Bz; + Bx = Bfield_fp[lev][0].get(); + By = Bfield_fp[lev][1].get(); + Bz = Bfield_fp[lev][2].get(); + + bool B_flag = 1; + InitializeExternalFieldsOnGridUsingParser(Bx, By, Bz, lev, B_flag); + } + + if (E_ext_grid_s == "parse_e_ext_grid_function") { + + std::vector<std::string> f; + // Parse Ex_external_grid_function + pp.getarr("Ex_external_grid_function(x,y,z)", f); + str_Ex_ext_grid_function.clear(); + for (auto const& s : f) { + str_Ex_ext_grid_function += s; + } + f.clear(); + + // Parse Ey_external_grid_function + pp.getarr("Ey_external_grid_function(x,y,z)", f); + str_Ey_ext_grid_function.clear(); + for (auto const& s : f) { + str_Ey_ext_grid_function += s; + } + f.clear(); + + // Parse Ez_external_grid_function + pp.getarr("Ez_external_grid_function(x,y,z)", f); + str_Ez_ext_grid_function.clear(); + for (auto const& s : f) { + str_Ez_ext_grid_function += s; + } + f.clear(); + + // Initialize Efield_fp with external function + MultiFab *Ex, *Ey, *Ez; + Ex = Efield_fp[lev][0].get(); + Ey = Efield_fp[lev][1].get(); + Ez = Efield_fp[lev][2].get(); + + bool B_flag = 0; + InitializeExternalFieldsOnGridUsingParser(Ex, Ey, Ez, lev, B_flag); + } if (lev > 0) { for (int i = 0; i < 3; ++i) { - Efield_aux[lev][i]->setVal(E_external_grid[i]); - Bfield_aux[lev][i]->setVal(B_external_grid[i]); - current_cp[lev][i]->setVal(0.0); - Efield_cp[lev][i]->setVal(E_external_grid[i]); - Bfield_cp[lev][i]->setVal(B_external_grid[i]); + if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") { + Bfield_aux[lev][i]->setVal(B_external_grid[i]); + Bfield_cp[lev][i]->setVal(B_external_grid[i]); + } + if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") { + Efield_aux[lev][i]->setVal(E_external_grid[i]); + Efield_cp[lev][i]->setVal(E_external_grid[i]); + } + } + if (B_ext_grid_s == "parse_b_ext_grid_function") { + MultiFab *Bx_aux, *By_aux, *Bz_aux; + Bx_aux = Bfield_aux[lev][0].get(); + By_aux = Bfield_aux[lev][1].get(); + Bz_aux = Bfield_aux[lev][2].get(); + + // Setting b_flag to 1 since we are initializing + // B_external on the grid. + bool B_flag = 1; + InitializeExternalFieldsOnGridUsingParser(Bx_aux, By_aux, + Bz_aux, lev, B_flag); + + MultiFab *Bx_cp, *By_cp, *Bz_cp; + Bx_cp = Bfield_cp[lev][0].get(); + By_cp = Bfield_cp[lev][1].get(); + Bz_cp = Bfield_cp[lev][2].get(); + + InitializeExternalFieldsOnGridUsingParser(Bx_cp, By_cp, + Bz_cp, lev, B_flag); + + } + if (E_ext_grid_s == "parse_e_ext_grid_function") { + + MultiFab *Ex_aux, *Ey_aux, *Ez_aux; + Ex_aux = Efield_aux[lev][0].get(); + Ey_aux = Efield_aux[lev][1].get(); + Ez_aux = Efield_aux[lev][2].get(); + // Setting b_flag to zero since we are initializing + // E_external on the grid here. + bool B_flag = 0; + InitializeExternalFieldsOnGridUsingParser(Ex_aux, Ey_aux, + Ez_aux, lev, B_flag); + + MultiFab *Ex_cp, *Ey_cp, *Ez_cp; + Ex_cp = Efield_cp[lev][0].get(); + Ey_cp = Efield_cp[lev][1].get(); + Ez_cp = Efield_cp[lev][2].get(); + + InitializeExternalFieldsOnGridUsingParser(Ex_cp, Ey_cp, + Ez_cp, lev, B_flag); + } } @@ -306,3 +513,138 @@ WarpX::InitLevelDataFFT (int lev, Real time) } #endif + +/* \brief + * This function initializes the E and B fields on each level + * using the parser and the user-defined function for the external fields. + * Depending on the bool value of the B_flag, the subroutine will parse + * the Bx_/By_/Bz_external_grid_function (if B_flag==1) + * or parse the Ex_/Ey_/Ez_external_grid_function (if B_flag==0). + * 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 +WarpX::InitializeExternalFieldsOnGridUsingParser ( + MultiFab *mfx, MultiFab *mfy, MultiFab *mfz, + const int lev, const bool B_flag) +{ + std::unique_ptr<ParserWrapper> xfield_parsewrap; + std::unique_ptr<ParserWrapper> yfield_parsewrap; + std::unique_ptr<ParserWrapper> zfield_parsewrap; + + if (B_flag == 1) { + xfield_parsewrap.reset(new ParserWrapper + (makeParser(str_Bx_ext_grid_function))); + yfield_parsewrap.reset(new ParserWrapper + (makeParser(str_By_ext_grid_function))); + zfield_parsewrap.reset(new ParserWrapper + (makeParser(str_Bz_ext_grid_function))); + } else { + xfield_parsewrap.reset(new ParserWrapper + (makeParser(str_Ex_ext_grid_function))); + yfield_parsewrap.reset(new ParserWrapper + (makeParser(str_Ey_ext_grid_function))); + zfield_parsewrap.reset(new ParserWrapper + (makeParser(str_Ez_ext_grid_function))); + } + + ParserWrapper *xfield_wrap = xfield_parsewrap.get(); + ParserWrapper *yfield_wrap = yfield_parsewrap.get(); + ParserWrapper *zfield_wrap = zfield_parsewrap.get(); + + const auto dx_lev = geom[lev].CellSizeArray(); + const RealBox& real_box = geom[lev].ProbDomain(); + for ( MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + IntVect x_nodal_flag, y_nodal_flag, z_nodal_flag; + if (B_flag == 1) { + x_nodal_flag = Bx_nodal_flag; + y_nodal_flag = By_nodal_flag; + z_nodal_flag = Bz_nodal_flag; + } else { + x_nodal_flag = Ex_nodal_flag; + y_nodal_flag = Ey_nodal_flag; + z_nodal_flag = Ez_nodal_flag; + } + const Box& tbx = mfi.tilebox(x_nodal_flag); + const Box& tby = mfi.tilebox(y_nodal_flag); + const Box& tbz = mfi.tilebox(z_nodal_flag); + + auto const& mfxfab = mfx->array(mfi); + auto const& mfyfab = mfy->array(mfi); + auto const& mfzfab = mfz->array(mfi); + + auto const& mfx_IndexType = (*mfx).ixType(); + auto const& mfy_IndexType = (*mfy).ixType(); + auto const& mfz_IndexType = (*mfz).ixType(); + + // Initialize IntVect based on the index type of multiFab + // 0 if cell-centered, 1 if node-centered. + IntVect mfx_type(AMREX_D_DECL(0,0,0)); + IntVect mfy_type(AMREX_D_DECL(0,0,0)); + IntVect mfz_type(AMREX_D_DECL(0,0,0)); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + mfx_type[idim] = mfx_IndexType.nodeCentered(idim); + mfy_type[idim] = mfy_IndexType.nodeCentered(idim); + mfz_type[idim] = mfz_IndexType.nodeCentered(idim); + } + + amrex::ParallelFor (tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // Shift required in the x-, y-, or z- position + // depending on the index type of the multifab + Real fac_x = (1.0 - mfx_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#else + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the x-component of the field. + mfxfab(i,j,k) = xfield_wrap->getField(x,y,z); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real fac_x = (1.0 - mfy_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==3) + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfx_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the y-component of the field. + mfyfab(i,j,k) = yfield_wrap->getField(x,y,z); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real fac_x = (1.0 - mfz_type[0]) * dx_lev[0]*0.5; + Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; +#if (AMREX_SPACEDIM==2) + Real y = 0.0; + Real fac_z = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==3) + Real fac_y = (1.0 - mfx_type[1]) * dx_lev[1]*0.5; + Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; + Real fac_z = (1.0 - mfz_type[2]) * dx_lev[2]*0.5; + Real z = k*dx_lev[2] + real_box.lo(2) + fac_z; +#endif + // Initialize the z-component of the field. + mfzfab(i,j,k) = zfield_wrap->getField(x,y,z); + }, + amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3 + ); + + } + +} |