From 9035ee165054e25aedf98d97f16786d8d5f2965a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 10 Jan 2020 10:52:31 -0800 Subject: Started implementing finite difference solver Added Yee algorithm --- .../FieldSolver/FiniteDifferenceSolver/EvolveB.cpp | 68 ++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp (limited to 'Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp') diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp new file mode 100644 index 000000000..e5514c697 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp @@ -0,0 +1,68 @@ +// TODO include statements + +FiniteDifferenceSolver::EvolveB ( VectorField Bfield, + ConstVectorField Efield, + amrex::Real dt ) { + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) + if (fdtd_algo == MaxwellSolverAlgo::Yee){ + EvolveB ( Bfield, Efield, dt ); + } else if (fdtd_algo == MaxwellSolverAlgo::CKC) { + EvolveB ( Bfield, Efield, dt ); + } else { + amrex::Abort("Unknown algorithm"); + } +) + +template +FiniteDifferenceSolver::EvolveB ( VectorField Bfield, + ConstVectorField Efield, + amrex::Real dt ) { + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + auto const& Bx = Bfield[0]->array(mfi); + auto const& By = Bfield[1]->array(mfi); + auto const& Bz = Bfield[2]->array(mfi); + auto const& Ex = Efield[0]->array(mfi); + auto const& Ey = Efield[1]->array(mfi); + auto const& Ez = Efield[2]->array(mfi); + + // Extract stencil coefficients + Real const* AMREX_RESTRICT coefs_x = stencil_coefs_x.dataPtr(); + Real const* AMREX_RESTRICT coefs_y = stencil_coefs_y.dataPtr(); + Real const* AMREX_RESTRICT coefs_z = stencil_coefs_z.dataPtr(); + + // Extract tileboxes for which to loop + const Box& tbx = mfi.tilebox(Bx_nodal_flag); + const Box& tby = mfi.tilebox(By_nodal_flag); + const Box& tbz = mfi.tilebox(Bz_nodal_flag); + + // Loop over the cells and update the fields + amrex::ParallelFor(tbx, tby, tbz, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bx(i, j, k) += dt * algo::UpwardDz( Ey, i, j, k, coefs_z) + - dt * algo::UpwardDy( Ez, i, j, k, coefs_y); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + By(i, j, k) += dt * algo::UpwardDx( Ez, i, j, k, coefs_x) + - dt * algo::UpwardDz( Ex, i, j, k, coefs_z); + }, + + [=] AMREX_GPU_DEVICE (int i, int j, int k){ + Bz(i, j, k) += dt * algo::UpwardDy( Ex, i, j, k, coefs_y) + - dt * algo::UpwardDx( Ey, i, j, k, coefs_x); + } + + ); + + } + +}; -- cgit v1.2.3