diff options
Diffstat (limited to 'Source/FieldSolver/ElectrostaticSolver.cpp')
| -rw-r--r-- | Source/FieldSolver/ElectrostaticSolver.cpp | 193 |
1 files changed, 179 insertions, 14 deletions
diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index ab2fc2b25..7aea1d4e4 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -266,15 +266,35 @@ WarpX::computePhiRZ (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho Array<LinOpBCType,AMREX_SPACEDIM> lobc, hibc; lobc[0] = LinOpBCType::Neumann; hibc[0] = LinOpBCType::Dirichlet; - if ( Geom(0).isPeriodic(1) ) { + std::array<bool,AMREX_SPACEDIM> dirichlet_flag; + dirichlet_flag[0] = false; + Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo, phi_bc_values_hi; + if ( WarpX::field_boundary_lo[1] == FieldBoundaryType::Periodic + && WarpX::field_boundary_hi[1] == FieldBoundaryType::Periodic ) { lobc[1] = LinOpBCType::Periodic; hibc[1] = LinOpBCType::Periodic; - } else { + dirichlet_flag[1] = false; + } else if ( WarpX::field_boundary_lo[1] == FieldBoundaryType::PEC + && WarpX::field_boundary_hi[1] == FieldBoundaryType::PEC ) { // Use Dirichlet boundary condition by default. // Ideally, we would often want open boundary conditions here. lobc[1] = LinOpBCType::Dirichlet; hibc[1] = LinOpBCType::Dirichlet; + + // set flag so we know which dimensions to fix the potential for + dirichlet_flag[1] = true; + // parse the input file for the potential at the current time + getPhiBC(1, phi_bc_values_lo[1], phi_bc_values_hi[1]); } + else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, + "Field boundary conditions have to be either periodic or PEC " + "when using the electrostatic solver" + ); + } + + // set the boundary potential values if needed + setPhiBC(phi, dirichlet_flag, phi_bc_values_lo, phi_bc_values_hi); // Define the linear operator (Poisson operator) MLNodeLaplacian linop( geom_scaled, boxArray(), dmap ); @@ -282,23 +302,22 @@ WarpX::computePhiRZ (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho linop.setSigma( lev, *sigma[lev] ); } + for (int lev=0; lev < rho.size(); lev++){ + rho[lev]->mult(-1._rt/PhysConst::ep0); + } + // Solve the Poisson equation linop.setDomainBC( lobc, hibc ); MLMG mlmg(linop); mlmg.setVerbose(2); mlmg.setMaxIter(max_iters); mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), required_precision, 0.0); - - // Normalize by the correct physical constant - for (int lev=0; lev < rho.size(); lev++){ - phi[lev]->mult(-1._rt/PhysConst::ep0); - } } #else /* Compute the potential `phi` in Cartesian geometry by solving the Poisson equation - hwith `rho` as a source, assuming that the source moves at a constant + with `rho` as a source, assuming that the source moves at a constant speed \f$\vec{\beta}\f$. This uses the amrex solver. @@ -321,20 +340,39 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> // Define the boundary conditions Array<LinOpBCType,AMREX_SPACEDIM> lobc, hibc; + std::array<bool,AMREX_SPACEDIM> dirichlet_flag; + Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo, phi_bc_values_hi; for (int idim=0; idim<AMREX_SPACEDIM; idim++){ - if ( Geom(0).isPeriodic(idim) ) { + if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic + && WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic ) { lobc[idim] = LinOpBCType::Periodic; hibc[idim] = LinOpBCType::Periodic; - } else { - // Use Dirichlet boundary condition by default. + dirichlet_flag[idim] = false; + } else if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC + && WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC ) { // Ideally, we would often want open boundary conditions here. lobc[idim] = LinOpBCType::Dirichlet; hibc[idim] = LinOpBCType::Dirichlet; + + // set flag so we know which dimensions to fix the potential for + dirichlet_flag[idim] = true; + // parse the input file for the potential at the current time + getPhiBC(idim, phi_bc_values_lo[idim], phi_bc_values_hi[idim]); + } + else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, + "Field boundary conditions have to be either periodic or PEC " + "when using the electrostatic solver" + ); } } + // set the boundary potential values if needed + setPhiBC(phi, dirichlet_flag, phi_bc_values_lo, phi_bc_values_hi); + // Define the linear operator (Poisson operator) MLNodeTensorLaplacian linop( Geom(), boxArray(), DistributionMap() ); + // Set the value of beta amrex::Array<amrex::Real,AMREX_SPACEDIM> beta_solver = #if (AMREX_SPACEDIM==2) @@ -346,18 +384,145 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> // Solve the Poisson equation linop.setDomainBC( lobc, hibc ); + + for (int lev=0; lev < rho.size(); lev++){ + rho[lev]->mult(-1._rt/PhysConst::ep0); + } + MLMG mlmg(linop); mlmg.setVerbose(2); mlmg.setMaxIter(max_iters); mlmg.solve( GetVecOfPtrs(phi), GetVecOfConstPtrs(rho), required_precision, 0.0); +} +#endif - // Normalize by the correct physical constant - for (int lev=0; lev < rho.size(); lev++){ - phi[lev]->mult(-1._rt/PhysConst::ep0); +/* \bried Set Dirichlet boundary conditions for the electrostatic solver. + + The given potential's values are fixed on the boundaries of the given + dimension according to the desired values from the simulation input file, + boundary.potential_lo and boundary.potential_hi. + + \param[inout] phi The electrostatic potential + \param[in] idim The dimension for which the Dirichlet boundary condition is set +*/ +void +WarpX::setPhiBC( amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi, + std::array<bool,AMREX_SPACEDIM> dirichlet_flag, + Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_lo, + Array<amrex::Real,AMREX_SPACEDIM> phi_bc_values_hi ) const +{ + // check if any dimension has Dirichlet boundary conditions + bool has_Dirichlet = false; + for (int idim=0; idim<AMREX_SPACEDIM; idim++){ + if (dirichlet_flag[idim]) { + has_Dirichlet = true; + } } + if (!has_Dirichlet) return; + + // loop over all mesh refinement levels and set the boundary values + for (int lev=0; lev <= max_level; lev++) { + + amrex::Box domain = Geom(lev).Domain(); + domain.surroundingNodes(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + // Extract the potential + auto phi_arr = phi[lev]->array(mfi); + // Extract tileboxes for which to loop + const Box& tb = mfi.tilebox( phi[lev]->ixType().toIntVect() ); + + // loop over dimensions + for (int idim=0; idim<AMREX_SPACEDIM; idim++){ + // check if the boundary in this dimension should be set + if (!dirichlet_flag[idim]) continue; + + // a check can be added below to test if the boundary values + // are already correct, in which case the ParallelFor over the + // cells can be skipped + + if (!domain.strictly_contains(tb)) { + amrex::ParallelFor( tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + + IntVect iv(AMREX_D_DECL(i,j,k)); + + if (iv[idim] == domain.smallEnd(idim)){ + phi_arr(i,j,k) = phi_bc_values_lo[idim]; + } + if (iv[idim] == domain.bigEnd(idim)) { + phi_arr(i,j,k) = phi_bc_values_hi[idim]; + } + + } // loop ijk + ); + } + } // idim + }} // lev & MFIter } + +/* \bried Utility function to parse input file for boundary potentials. + + The input values are parsed to allow math expressions for the potentials + that specify time dependence. + + \param[in] idim The dimension for which the potential is queried + \param[inout] pot_lo The specified value of `phi` on the lower boundary. + \param[inout] pot_hi The specified value of `phi` on the upper boundary. +*/ +void +WarpX::getPhiBC( const int idim, amrex::Real &pot_lo, amrex::Real &pot_hi ) const +{ + // set default potentials to zero in order for current tests to pass + // but forcing the user to specify a potential might be better + std::string potential_lo_str = "0"; + std::string potential_hi_str = "0"; + + // Get the boundary potentials specified in the simulation input file + // first as strings and then parse those for possible math expressions + ParmParse pp_boundary("boundary"); + +#ifdef WARPX_DIM_RZ + if (idim == 1) { + pp_boundary.query("potential_lo_z", potential_lo_str); + pp_boundary.query("potential_hi_z", potential_hi_str); + } + else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, + "Field boundary condition values can currently only be specified " + "for z when using RZ geometry." + ); + } +#else + if (idim == 0) { + pp_boundary.query("potential_lo_x", potential_lo_str); + pp_boundary.query("potential_hi_x", potential_hi_str); + } + else if (idim == 1){ + if (AMREX_SPACEDIM == 2){ + pp_boundary.query("potential_lo_z", potential_lo_str); + pp_boundary.query("potential_hi_z", potential_hi_str); + } + else { + pp_boundary.query("potential_lo_y", potential_lo_str); + pp_boundary.query("potential_hi_y", potential_hi_str); + } + } + else { + pp_boundary.query("potential_lo_z", potential_lo_str); + pp_boundary.query("potential_hi_z", potential_hi_str); + } #endif + auto parser_lo = makeParser(potential_lo_str, {"t"}); + pot_lo = parser_lo.eval(gett_new(0)); + auto parser_hi = makeParser(potential_hi_str, {"t"}); + pot_hi = parser_hi.eval(gett_new(0)); +} + /* \bried Compute the electric field that corresponds to `phi`, and add it to the set of MultiFab `E`. |