From 74cffc29f41ff424fd987c81d4fb71ddfbfb711b Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 2 May 2019 17:01:31 -0700 Subject: Start implementation of spectral PML --- Source/FieldSolver/SpectralSolver/SpectralSolver.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralSolver.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index c21c3cfb1..a91fcbc47 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -30,6 +30,7 @@ SpectralSolver::SpectralSolver( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); // - Initialize arrays for fields in spectral space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); + field_data = SpectralFieldData( realspace_ba, k_space, dm, + algorithm->getRequiredNumberOfFields() ); }; -- cgit v1.2.3 From 257e71c2eaca05416122e8b748c13a324e82d831 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 3 May 2019 15:45:09 -0700 Subject: Prepare spectralsolver for PML --- Source/BoundaryConditions/PML.H | 9 +++++++++ Source/BoundaryConditions/PML.cpp | 1 - Source/FieldSolver/SpectralSolver/SpectralSolver.H | 3 ++- Source/FieldSolver/SpectralSolver/SpectralSolver.cpp | 20 ++++++++++++++++---- 4 files changed, 27 insertions(+), 6 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralSolver.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 25dfc7996..39dfa5a25 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -6,6 +6,10 @@ #include #include +#ifdef WARPX_USE_PSATD +#include +#endif + #if (AMREX_SPACEDIM == 3) #define WRPX_PML_TO_FORTRAN(x) \ @@ -161,6 +165,11 @@ private: std::unique_ptr sigba_fp; std::unique_ptr sigba_cp; +#ifdef WARPX_USE_PSATD + std::unique_ptr spectral_solver_fp; + std::unique_ptr spectral_solver_cp; +#endif + static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index d8052c4fa..cac43faa5 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -425,7 +425,6 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); } - } BoxArray diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index d4019a9a3..c570b017b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -23,7 +23,8 @@ class SpectralSolver const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ); + const amrex::RealVect dx, const amrex::Real dt, + const bool pml=false ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index a91fcbc47..80555a7b3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -7,13 +7,22 @@ * This function selects the spectral algorithm to be used, allocates the * corresponding coefficients for the discretized field update equation, * and prepares the structures that store the fields in spectral space. + * + * \param norder_x Order of accuracy of the spatial derivatives along x + * \param norder_y Order of accuracy of the spatial derivatives along y + * \param norder_z Order of accuracy of the spatial derivatives along z + * \param nodal Whether the solver is applied to a nodal or staggered grid + * \param dx Cell size along each dimension + * \param dt Time step + * \param pml Whether the boxes in which the solver is applied are PML boxes */ SpectralSolver::SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { + const amrex::RealVect dx, const amrex::Real dt, + const bool pml ) { // Initialize all structures using the same distribution mapping dm @@ -24,10 +33,13 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - // TODO: Add more algorithms + selection depending on input parameters - // For the moment, this only uses the standard PsatdAlgorithm - algorithm = std::unique_ptr( new PsatdAlgorithm( + if (pml) { + algorithm = std::unique_ptr( new PMLPsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } else { + algorithm = std::unique_ptr( new PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } // - Initialize arrays for fields in spectral space + FFT plans field_data = SpectralFieldData( realspace_ba, k_space, dm, -- cgit v1.2.3 From f71b9c0c409ebf1750b9c40e861bfeefac4247ca Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 4 May 2019 17:57:14 -0700 Subject: Initialize spectral solver --- Source/BoundaryConditions/PML.H | 4 +++- Source/BoundaryConditions/PML.cpp | 18 +++++++++++++++++- .../FieldSolver/SpectralSolver/SpectralSolver.cpp | 1 + Source/Initialization/WarpXInitData.cpp | 21 ++++++++++++--------- 4 files changed, 33 insertions(+), 11 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralSolver.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 39dfa5a25..15c70c10b 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -104,7 +104,9 @@ class PML public: PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::Geometry* geom, const amrex::Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window); + int ncell, int delta, int ref_ratio, amrex::Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window); void ComputePMLFactors (amrex::Real dt); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index cac43faa5..831339a2f 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -329,7 +329,9 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window) + int ncell, int delta, int ref_ratio, Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window) : m_geom(geom), m_cgeom(cgeom) { @@ -390,6 +392,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); +#ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))}; + spectral_solver_fp.reset( new SpectralSolver( ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, dx, dt, in_pml ) ); +#endif + if (cgeom) { #ifndef WARPX_USE_PSATD @@ -424,6 +433,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + + #ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; + spectral_solver_fp.reset( new SpectralSolver( cba, cdm, + nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); + #endif } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 80555a7b3..4b9def013 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -1,6 +1,7 @@ #include #include #include +#include /* \brief Initialize the spectral Maxwell solver * diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..f59f5d81f 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -87,7 +87,7 @@ WarpX::InitDiagnostics () { const Real* current_lo = geom[0].ProbLo(); const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; - + // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); @@ -96,7 +96,7 @@ WarpX::InitDiagnostics () { zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, gamma_boost, - t_new[0], dt_boost, + t_new[0], dt_boost, moving_window_dir, geom[0])); } } @@ -117,10 +117,10 @@ WarpX::InitFromScratch () InitPML(); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC if (do_electrostatic) { getLevelMasks(masks); - + // the plus one is to convert from num_cells to num_nodes getLevelMasks(gather_masks, n_buffer + 1); } @@ -133,13 +133,16 @@ WarpX::InitPML () if (do_pml) { pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, - pml_ncell, pml_delta, 0, do_dive_cleaning, do_moving_window)); + pml_ncell, pml_delta, 0, dt[0], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), - pml_ncell, pml_delta, refRatio(lev-1)[0], do_dive_cleaning, - do_moving_window)); + pml_ncell, pml_delta, refRatio(lev-1)[0], dt[lev], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); } } } @@ -222,7 +225,7 @@ WarpX::InitOpenbc () Vector alllohi(6*nprocs,100000); MPI_Allgather(lohi, 6, MPI_INT, alllohi.data(), 6, MPI_INT, ParallelDescriptor::Communicator()); - + BoxList bl{IndexType::TheNodeType()}; for (int i = 0; i < nprocs; ++i) { @@ -248,7 +251,7 @@ WarpX::InitOpenbc () rho_openbc.copy(*rho, 0, 0, 1, rho->nGrow(), 0, gm.periodicity(), FabArrayBase::ADD); const Real* dx = gm.CellSize(); - + warpx_openbc_potential(rho_openbc[myproc].dataPtr(), phi_openbc[myproc].dataPtr(), dx); BoxArray nba = boxArray(lev); -- cgit v1.2.3