From 2e9f320008b0f03d6d91ba6349f93f96196484dd Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 11 Nov 2024 14:02:50 -0800 Subject: [PATCH] Remove SWFFT, FFTW, and heFFTe tutorials and documentation (#141) * remove SWFFT, FFTW, and heFFTe tutorials and documentation; they are superceded by amrex::FFT * fix broken figure --- Docs/source/FFTW_Tutorial.rst | 48 -- Docs/source/MUI_Tutorial.rst | 2 +- Docs/source/SWFFT_Tutorial.rst | 86 ---- Docs/source/{SWFFT => figs}/iface_rect.png | Bin Docs/source/heFFTe_Tutorial.rst | 48 -- Docs/source/index.rst | 12 - ExampleCodes/FFTW/Basic/GNUmakefile | 22 - ExampleCodes/FFTW/Basic/Make.package | 3 - ExampleCodes/FFTW/Basic/inputs | 10 - ExampleCodes/FFTW/Basic/main.cpp | 420 ---------------- ExampleCodes/FFTW/Basic/myfunc.H | 18 - ExampleCodes/FFTW/Basic/myfunc.cpp | 91 ---- ExampleCodes/FFTW/Poisson/GNUmakefile | 22 - ExampleCodes/FFTW/Poisson/Make.package | 3 - ExampleCodes/FFTW/Poisson/inputs | 14 - ExampleCodes/FFTW/Poisson/main.cpp | 447 ----------------- ExampleCodes/FFTW/Poisson/myfunc.H | 16 - ExampleCodes/FFTW/Poisson/myfunc.cpp | 20 - ExampleCodes/SWFFT/SWFFT_poisson/GNUmakefile | 22 - ExampleCodes/SWFFT/SWFFT_poisson/Make.package | 7 - ExampleCodes/SWFFT/SWFFT_poisson/README | 19 - ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.H | 33 -- .../SWFFT/SWFFT_poisson/SWFFT_Test.cpp | 173 ------- .../SWFFT/SWFFT_poisson/SWFFT_Test_F.F90 | 87 ---- .../SWFFT/SWFFT_poisson/SWFFT_Test_F.H | 21 - ExampleCodes/SWFFT/SWFFT_poisson/inputs.128 | 3 - ExampleCodes/SWFFT/SWFFT_poisson/inputs.32 | 3 - ExampleCodes/SWFFT/SWFFT_poisson/inputs.64 | 3 - ExampleCodes/SWFFT/SWFFT_poisson/main.cpp | 16 - ExampleCodes/SWFFT/SWFFT_poisson/run_me | 3 - .../SWFFT/SWFFT_poisson/swfft_solver.cpp | 169 ------- ExampleCodes/SWFFT/SWFFT_simple/GNUmakefile | 22 - ExampleCodes/SWFFT/SWFFT_simple/Make.package | 7 - ExampleCodes/SWFFT/SWFFT_simple/README | 36 -- ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.H | 33 -- .../SWFFT/SWFFT_simple/SWFFT_Test.cpp | 209 -------- .../SWFFT/SWFFT_simple/SWFFT_Test_F.F90 | 133 ----- .../SWFFT/SWFFT_simple/SWFFT_Test_F.H | 21 - .../SWFFT/SWFFT_simple/inputs.multipleGrids | 13 - .../SWFFT/SWFFT_simple/inputs.oneGrid | 13 - ExampleCodes/SWFFT/SWFFT_simple/main.cpp | 16 - ExampleCodes/SWFFT/SWFFT_simple/run_me_2d | 2 - ExampleCodes/SWFFT/SWFFT_simple/run_me_3d | 2 - .../SWFFT/SWFFT_simple/swfft_compute.cpp | 162 ------ ExampleCodes/heFFTe/Basic/GNUmakefile | 35 -- ExampleCodes/heFFTe/Basic/Make.package | 1 - ExampleCodes/heFFTe/Basic/README | 41 -- ExampleCodes/heFFTe/Basic/inputs | 15 - ExampleCodes/heFFTe/Basic/main.cpp | 462 ------------------ ExampleCodes/heFFTe/Poisson/GNUmakefile | 35 -- ExampleCodes/heFFTe/Poisson/Make.package | 1 - ExampleCodes/heFFTe/Poisson/README | 41 -- ExampleCodes/heFFTe/Poisson/inputs | 18 - ExampleCodes/heFFTe/Poisson/main.cpp | 307 ------------ ExampleCodes/heFFTe/Poisson2/GNUmakefile | 35 -- ExampleCodes/heFFTe/Poisson2/Make.package | 1 - ExampleCodes/heFFTe/Poisson2/inputs | 11 - ExampleCodes/heFFTe/Poisson2/main.cpp | 314 ------------ 58 files changed, 1 insertion(+), 3826 deletions(-) delete mode 100644 Docs/source/FFTW_Tutorial.rst delete mode 100644 Docs/source/SWFFT_Tutorial.rst rename Docs/source/{SWFFT => figs}/iface_rect.png (100%) delete mode 100644 Docs/source/heFFTe_Tutorial.rst delete mode 100644 ExampleCodes/FFTW/Basic/GNUmakefile delete mode 100644 ExampleCodes/FFTW/Basic/Make.package delete mode 100644 ExampleCodes/FFTW/Basic/inputs delete mode 100644 ExampleCodes/FFTW/Basic/main.cpp delete mode 100644 ExampleCodes/FFTW/Basic/myfunc.H delete mode 100644 ExampleCodes/FFTW/Basic/myfunc.cpp delete mode 100644 ExampleCodes/FFTW/Poisson/GNUmakefile delete mode 100644 ExampleCodes/FFTW/Poisson/Make.package delete mode 100644 ExampleCodes/FFTW/Poisson/inputs delete mode 100644 ExampleCodes/FFTW/Poisson/main.cpp delete mode 100644 ExampleCodes/FFTW/Poisson/myfunc.H delete mode 100644 ExampleCodes/FFTW/Poisson/myfunc.cpp delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/GNUmakefile delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/Make.package delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/README delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.H delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.cpp delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.F90 delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.H delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/inputs.128 delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/inputs.32 delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/inputs.64 delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/main.cpp delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/run_me delete mode 100644 ExampleCodes/SWFFT/SWFFT_poisson/swfft_solver.cpp delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/GNUmakefile delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/Make.package delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/README delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.H delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.cpp delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.F90 delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.H delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/inputs.multipleGrids delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/inputs.oneGrid delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/main.cpp delete mode 100755 ExampleCodes/SWFFT/SWFFT_simple/run_me_2d delete mode 100755 ExampleCodes/SWFFT/SWFFT_simple/run_me_3d delete mode 100644 ExampleCodes/SWFFT/SWFFT_simple/swfft_compute.cpp delete mode 100644 ExampleCodes/heFFTe/Basic/GNUmakefile delete mode 100644 ExampleCodes/heFFTe/Basic/Make.package delete mode 100644 ExampleCodes/heFFTe/Basic/README delete mode 100644 ExampleCodes/heFFTe/Basic/inputs delete mode 100644 ExampleCodes/heFFTe/Basic/main.cpp delete mode 100644 ExampleCodes/heFFTe/Poisson/GNUmakefile delete mode 100644 ExampleCodes/heFFTe/Poisson/Make.package delete mode 100644 ExampleCodes/heFFTe/Poisson/README delete mode 100644 ExampleCodes/heFFTe/Poisson/inputs delete mode 100644 ExampleCodes/heFFTe/Poisson/main.cpp delete mode 100644 ExampleCodes/heFFTe/Poisson2/GNUmakefile delete mode 100644 ExampleCodes/heFFTe/Poisson2/Make.package delete mode 100644 ExampleCodes/heFFTe/Poisson2/inputs delete mode 100644 ExampleCodes/heFFTe/Poisson2/main.cpp diff --git a/Docs/source/FFTW_Tutorial.rst b/Docs/source/FFTW_Tutorial.rst deleted file mode 100644 index 19370a1b..00000000 --- a/Docs/source/FFTW_Tutorial.rst +++ /dev/null @@ -1,48 +0,0 @@ -.. role:: cpp(code) - :language: c++ - -.. role:: fortran(code) - :language: fortran - -.. _tutorials_fftw: - -FFTW -========================== - -These tutorials demonstrate how to call fftw3 (CPU) or cuFFT (GPU) to solve for and manipulate Fourier transform data using a single MPI rank. - -There are two FFTW tutorials, ``Basic`` and ``Poisson``: - -- ``Basic`` tutorial: The tutorial found in ``amrex-tutorials/ExampleCodes/FFTW/Basic`` is - useful if the objective is to simply take a forward FFT of data, and the DFT's ordering in k-space - matters to the user. This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, - and then redistributes the data in k-space where the center cell in the domain corresponds to the k=0 mode. - The results are written to a plot file. - -- ``Poisson`` tutorial: This tutorial: ``amrex-tutorials/ExampleCodes/FFTW/Poisson`` - solves a Poisson equation with periodic boundary conditions. In it, both a forward FFT and reverse FFT - are called to solve the equation, however, no reordering of the DFT data in k-space is performed. - -We note that both fftw and cufft assume a row-major ordering of data; since a :cpp:`MultiFab` is column major, -the output to the spectral array is spatially-transposed. - -.. toctree:: - :maxdepth: 1 - -.. _section:fftw_tutorial:fftw_basic: - -Basic --------------------------- - -This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, -and then redistributes the data in k-space where the center cell in the domain corresponds to the k=0 mode. -The results are written to a plot file. - -.. _section:fftw_tutorial:fftw_pois: - -Poisson --------------------------- - -In this test case we set up a right hand side (rhs), call the forward transform, -modify the coefficients, then call the backward solver and output the solution -to the discrete Poisson equation. diff --git a/Docs/source/MUI_Tutorial.rst b/Docs/source/MUI_Tutorial.rst index ded1738f..f1de5aab 100644 --- a/Docs/source/MUI_Tutorial.rst +++ b/Docs/source/MUI_Tutorial.rst @@ -34,7 +34,7 @@ as I am aware, the code works for any AMReX grid structure. Details of how to bu The figure below shows one possible grid structure of the 2D (red grid) and 3D (multicolored blocks) setup. -.. |a| image:: ./SWFFT/iface_rect.png +.. |a| image:: ./figs/iface_rect.png :width: 60% .. table:: diff --git a/Docs/source/SWFFT_Tutorial.rst b/Docs/source/SWFFT_Tutorial.rst deleted file mode 100644 index 15782193..00000000 --- a/Docs/source/SWFFT_Tutorial.rst +++ /dev/null @@ -1,86 +0,0 @@ -.. role:: cpp(code) - :language: c++ - -.. role:: fortran(code) - :language: fortran - -.. _tutorials_swfft: - -SWFFT -========================== - -This Tutorial demonstrates how to call the SWFFT wrapper to the FFTW3 solver. - -Note that the SWFFT source code was developed by Adrian Pope and colleagues and -is available at: - -https://git.cels.anl.gov/hacc/SWFFT - -Please refer to the AMReX documentation at :ref:`amrex_docs:swfftdoc` for a brief explanation -of how the SWFFT redistributes data into pencil grids. - -AMReX contains two SWFFT tutorials, ``SWFFT_poisson`` and ``SWFFT_simple``: - -- ``SWFFT_poisson`` tutorial: The tutorial found in ``amrex-tutorials/ExampleCodes/SWFFT/SWFFT_poisson`` solves - a Poisson equation with periodic boundary conditions. In it, both a forward FFT and reverse FFT - are called to solve the equation, however, no reordering of the DFT data in k-space is performed. - -- ``SWFFT_simple`` tutorial: This tutorial: ``amrex-tutorials/ExampleCodes/SWFFT/SWFFT_simple``, is - useful if the objective is to simply take a forward FFT of data, and the DFT's ordering in k-space - matters to the user. This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, - and then redistributes the data in k-space back to the "correct," 0 to :math:`2\pi`, ordering. - The results are written to a plot file. - -.. toctree:: - :maxdepth: 1 - -.. _section:swfft_tutorial:swfft_pois: - -SWFFT_poisson --------------------------- - -In this test case we set up a right hand side (rhs), call the forward transform, -modify the coefficients, then call the backward solver and output the solution -to the discrete Poisson equation. - -To build the code, type 'make' in ``amrex-tutorials/ExampleCodes/SWFFT/SWFFT_poisson``. This -will include code from ``amrex/Src/Extern/SWFFT`` and you will need to -link to the FFT solvers themselves (on NERSC's Cori machine, for example, -you would need to "module load fft") - -To run the code, type 'main3d.gnu.MPI.ex inputs' in this directory - -To visualize the output, set the bool write_data to true, then -use amrvis3d (source available at https://github.com/AMReX-Codes/Amrvis): - -amrvis3d -mf RHS SOL_EXACT SOL_COMP - -to visualize the rhs, the exact solution and the computed solution. - -The max norm of the difference between the exact and computed solution is also printed. - -For instructions on how to take a forward FFT only using SWFFT, please refer to :ref:`section:swfft_tutorial:swfft_simple`. - -.. _section:swfft_tutorial:swfft_simple: - -SWFFT_simple --------------------------- - -This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, and then redistributes the data in k-space back to the "correct," 0 to :math:`2\pi`, ordering. The results are written to a plot file. - -In a similar fashion to the ``SWFFT_poisson`` tutorial: - -To build the code, type 'make' in ``amrex-tutorials/ExampleCodes/SWFFT/SWFFT_simple``. This -will include code from ``amrex/Src/Extern/SWFFT`` and you will need to -link to the FFT solvers themselves (on NERSC's Cori machine, for example, -you would need to "module load fft") - -To run the code, type 'main*.ex inputs.oneGrid' in this directory to run the code in serial. To run the code in parallel, type 'mpiexec -n $N main*.ex inputs.multipleGrids' instead, where ``N`` holds the number of MPI processes (equal to the number of grids). ``run_me_2d`` and ``run_me_3d`` also provide examples of how to run the code. - -Use amrvis2d or amrvis3d to visualize the output (source available at https://github.com/AMReX-Codes/Amrvis): - -amrvis${dims}d plt_fft* - -where ``dims`` specifies ``AMREX_SPACEDIM``. The DFT of the data and the original data are labeled as ``FFT_of_phi`` and ``phi`` within the plot file. - -The :ref:`section:swfft_tutorial:swfft_pois` tutorial provides an example of solving a Poisson equation using a discrete spectral method, in which a forward and reverse FFT of a :cpp:`MultiFab` are computed. diff --git a/Docs/source/SWFFT/iface_rect.png b/Docs/source/figs/iface_rect.png similarity index 100% rename from Docs/source/SWFFT/iface_rect.png rename to Docs/source/figs/iface_rect.png diff --git a/Docs/source/heFFTe_Tutorial.rst b/Docs/source/heFFTe_Tutorial.rst deleted file mode 100644 index 7bd081c0..00000000 --- a/Docs/source/heFFTe_Tutorial.rst +++ /dev/null @@ -1,48 +0,0 @@ -.. role:: cpp(code) - :language: c++ - -.. role:: fortran(code) - :language: fortran - -.. _tutorials_heffte: - -heFFTe -========================== - -These tutorials demonstrate how to heFFTe to solve for and manipulate Fourier transform data using multiple MPI ranks/GPUs. - -There are two heFFTe tutorials, ``Basic`` and ``Poisson``: - -- ``Basic`` tutorial: The tutorial found in ``amrex-tutorials/ExampleCodes/heFFTe/Basic`` is - useful if the objective is to simply take a forward FFT of data, and the DFT's ordering in k-space - matters to the user. This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, - and then redistributes the data in k-space where the center cell in the domain corresponds to the k=0 mode. - The results are written to a plot file. - -- ``Poisson`` tutorial: This tutorial: ``amrex-tutorials/ExampleCodes/heFFTe/Poisson`` - solves a Poisson equation with periodic boundary conditions. In it, both a forward FFT and reverse FFT - are called to solve the equation, however, no reordering of the DFT data in k-space is performed. - -We note that both backends, fftw and cufft, assume a row-major ordering of data; since a :cpp:`MultiFab` is column major, -the output to the spectral array is spatially-transposed. - -.. toctree:: - :maxdepth: 1 - -.. _section:heffte_tutorial:heffte_basic: - -Basic --------------------------- - -This tutorial initializes a 3D or 2D :cpp:`MultiFab`, takes a forward FFT, -and then redistributes the data in k-space where the center cell in the domain corresponds to the k=0 mode. -The results are written to a plot file. - -.. _section:heffte_tutorial:heffte_pois: - -Poisson --------------------------- - -In this test case we set up a right hand side (rhs), call the forward transform, -modify the coefficients, then call the backward solver and output the solution -to the discrete Poisson equation. diff --git a/Docs/source/index.rst b/Docs/source/index.rst index b4180b10..f185fc6a 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -46,10 +46,8 @@ sorted by the following categories: and analysis tool. - :ref:`EB` -- Examples of embedded boundaries. - :ref:`FFT` -- Demonstration on how to use the amrex::FFT class for forward/inverse FFT and data manipulation -- :ref:`FFTW` -- FFTW and cuFFT single-rank tutorials. - :ref:`ForkJoin` -- Parallel execution and subgrouping of MPI ranks. - :ref:`GPU` -- Offload work to the GPUs using AMReX tools. -- :ref:`heFFTe` -- heFFTe distributed tutorials. - :ref:`Linear Solvers` -- Examples of several linear solvers. - :ref:`ML/PYTORCH` -- Use of pytorch models to replace point-wise computational kernels. - :ref:`MPMD` -- Usage of AMReX-MPMD (Multiple Program Multiple Data) framework. @@ -60,8 +58,6 @@ sorted by the following categories: to solve a scalar advection-diffusion-reaction equation. - :ref:`SENSEI` -- In situ data analysis and visualization through a unified interface. - :ref:`SUNDIALS` -- Time integration with SUNDIALS backend and native AMReX types. -- :ref:`SWFFT` -- Demonstrates how to call the SWFFT wrapper to the FFTW3 (A distributed memory - implementation of the discrete Fourier transform) solver. .. toctree:: @@ -72,10 +68,8 @@ sorted by the following categories: Blueprint_Tutorial EB_Tutorial FFT_Tutorial - FFTW_Tutorial ForkJoin_Tutorial GPU_Tutorial - heFFTe_Tutorial LinearSolvers_Tutorial ML_Tutorial MPMD_Tutorials @@ -85,7 +79,6 @@ sorted by the following categories: SDC_Tutorial SENSEI_Tutorial SUNDIALS_Tutorial - SWFFT_Tutorial .. _`AMR`: AMR_Tutorial.html @@ -98,14 +91,10 @@ sorted by the following categories: .. _`FFT`: FFT_Tutorial.html -.. _`FFTW`: FFTW_Tutorial.html - .. _`ForkJoin`: ForkJoin_Tutorial.html .. _`GPU`: GPU_Tutorial.html -.. _`heFFTe`: heFFTe_Tutorial.html - .. _`Linear Solvers`: LinearSolvers_Tutorial.html .. _`MPMD`: MPMD_Tutorials.html @@ -120,7 +109,6 @@ sorted by the following categories: .. _`SENSEI`: SENSEI_Tutorial.html -.. _`SWFFT`: SWFFT_Tutorial.html diff --git a/ExampleCodes/FFTW/Basic/GNUmakefile b/ExampleCodes/FFTW/Basic/GNUmakefile deleted file mode 100644 index e6680cc3..00000000 --- a/ExampleCodes/FFTW/Basic/GNUmakefile +++ /dev/null @@ -1,22 +0,0 @@ -# AMREX_HOME defines the directory in which we will find all the AMReX code. -AMREX_HOME ?= ../../../../amrex - -USE_MPI = TRUE -DEBUG = FALSE -USE_OMP = FALSE -COMP = gnu -DIM = 2 - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package - -include $(AMREX_HOME)/Src/Base/Make.package - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules - -ifeq ($(USE_CUDA),TRUE) - LIBRARIES += -lcufft -else - LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3 -endif diff --git a/ExampleCodes/FFTW/Basic/Make.package b/ExampleCodes/FFTW/Basic/Make.package deleted file mode 100644 index 46fa6f45..00000000 --- a/ExampleCodes/FFTW/Basic/Make.package +++ /dev/null @@ -1,3 +0,0 @@ -CEXE_sources += main.cpp myfunc.cpp -CEXE_headers = myfunc.H - diff --git a/ExampleCodes/FFTW/Basic/inputs b/ExampleCodes/FFTW/Basic/inputs deleted file mode 100644 index e48ea95b..00000000 --- a/ExampleCodes/FFTW/Basic/inputs +++ /dev/null @@ -1,10 +0,0 @@ -amrex.the_arena_is_managed = 1 - -max_grid_size = 32 -n_cell_x = 64 -n_cell_y = 64 -n_cell_z = 64 - -prob_hi_x = 1. -prob_hi_y = 1. -prob_hi_z = 1. diff --git a/ExampleCodes/FFTW/Basic/main.cpp b/ExampleCodes/FFTW/Basic/main.cpp deleted file mode 100644 index 454b366c..00000000 --- a/ExampleCodes/FFTW/Basic/main.cpp +++ /dev/null @@ -1,420 +0,0 @@ -#include - -#include -#include - -#ifdef AMREX_USE_CUDA -#include -#else -#include -#include -#endif - -#include "myfunc.H" - -using namespace amrex; - -int main (int argc, char* argv[]) -{ - - Initialize(argc,argv); - { - - // store the current time so we can later compute total run time. - Real start_time = ParallelDescriptor::second(); - - // ********************************** - // DECLARE SIMULATION PARAMETERS - // ********************************** - - // number of cells on each side of the domain - int n_cell_x; - int n_cell_y; - int n_cell_z; - - // dimensions of the domain - Real prob_hi_x; - Real prob_hi_y; - Real prob_hi_z; - - // This is the largest size a grid can be - int max_grid_size; - - // ********************************** - // READ PARAMETER VALUES FROM INPUTS FILE - // ********************************** - { - // ParmParse is way of reading inputs from the inputs file - // pp.get means we require the inputs file to have it - // pp.query means we optionally need the inputs file to have it - but you should supply a default value above - ParmParse pp; - - pp.get("n_cell_x",n_cell_x); - pp.get("n_cell_y",n_cell_y); - pp.get("n_cell_z",n_cell_z); - - pp.get("prob_hi_x",prob_hi_x); - pp.get("prob_hi_y",prob_hi_y); - pp.get("prob_hi_z",prob_hi_z); - - pp.get("max_grid_size",max_grid_size); - } - - // ********************************** - // DEFINE SIMULATION SETUP AND GEOMETRY - // ********************************** - // make BoxArray and Geometry - // ba will contain a list of boxes that cover the domain - // geom contains information such as the physical domain size, - // number of points in the domain, and periodicity - BoxArray ba; - Geometry geom; - - // define lower and upper indices - IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); - IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1)); - - // Make a single box that is the entire domain - Box domain(dom_lo, dom_hi); - - // Initialize the boxarray "ba" from the single box "domain" - ba.define(domain); - - // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction - ba.maxSize(max_grid_size); - - // How Boxes are distrubuted among MPI processes - DistributionMapping dm(ba); - - // This defines the physical box size in each direction - RealBox real_box({ AMREX_D_DECL( 0., 0., 0.)}, - { AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} ); - - // periodic in all direction - Array is_periodic{AMREX_D_DECL(1,1,1)}; - - // This defines a Geometry object - geom.define(domain, real_box, CoordSys::cartesian, is_periodic); - - // extract dx from the geometry object - GpuArray dx = geom.CellSizeArray(); - - // MultiFab storage for phi, and the real and imaginary parts of the dft - MultiFab phi (ba, dm, 1, 0); - MultiFab phi_dft_real(ba, dm, 1, 0); - MultiFab phi_dft_imag(ba, dm, 1, 0); - - // ********************************** - // INITIALIZE PHI - // ********************************** - - Real omega = M_PI/2.0; - - // loop over boxes - for (MFIter mfi(phi); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - - const Array4& phi_ptr = phi.array(mfi); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - // ********************************** - // SET VALUES FOR EACH CELL - // ********************************** - - Real x = (i+0.5) * dx[0]; - Real y = (AMREX_SPACEDIM>=2) ? (j+0.5) * dx[1] : 0.; - Real z = (AMREX_SPACEDIM==3) ? (k+0.5) * dx[2] : 0.; - phi_ptr(i,j,k) = std::sin(4*M_PI*x/prob_hi_x + omega); - if (AMREX_SPACEDIM >= 2) { - phi_ptr(i,j,k) *= std::sin(6*M_PI*y/prob_hi_y + omega); - } - if (AMREX_SPACEDIM == 3) { - phi_ptr(i,j,k) *= std::sin(8*M_PI*z/prob_hi_z + omega); - } - }); - } - - // ********************************** - // COPY PHI INTO A MULTIFAB WITH ONE BOX - // ********************************** - - // create a new BoxArray and DistributionMapping for a MultiFab with 1 box - BoxArray ba_onegrid(geom.Domain()); - DistributionMapping dm_onegrid(ba_onegrid); - - // storage for phi and the dft - MultiFab phi_onegrid (ba_onegrid, dm_onegrid, 1, 0); - MultiFab phi_dft_real_onegrid(ba_onegrid, dm_onegrid, 1, 0); - MultiFab phi_dft_imag_onegrid(ba_onegrid, dm_onegrid, 1, 0); - - // copy phi into phi_onegrid - phi_onegrid.ParallelCopy(phi, 0, 0, 1); - - // ********************************** - // COMPUTE FFT - // ********************************** - -#ifdef AMREX_USE_CUDA - using FFTplan = cufftHandle; - using FFTcomplex = cuDoubleComplex; -#else - using FFTplan = fftw_plan; - using FFTcomplex = fftw_complex; -#endif - - // number of points in the domain - long npts = domain.numPts(); - Real sqrtnpts = std::sqrt(npts); - - // contain to store FFT - note it is shrunk by approximately a half in x - Vector > > > spectral_field; - - Vector forward_plan; - - for (MFIter mfi(phi_onegrid); mfi.isValid(); ++mfi) { - - // grab a single box - Box realspace_bx = mfi.fabbox(); - - // size of box - IntVect fft_size = realspace_bx.length(); // This will be different for FFTs of complex data - - // this is the size of the box, except the 0th component is 'halved plus 1' - IntVect spectral_bx_size = fft_size; - spectral_bx_size[0] = fft_size[0]/2 + 1; - - // spectral box - Box spectral_bx = Box(IntVect(0), spectral_bx_size - IntVect(1)); - - spectral_field.emplace_back(new BaseFab >(spectral_bx,1, - The_Device_Arena())); - spectral_field.back()->setVal(0.0); // touch the memory - - FFTplan fplan; - -#ifdef AMREX_USE_CUDA - -#if (AMREX_SPACEDIM == 1) - cufftResult result = cufftPlan1d(&fplan, fft_size[0], CUFFT_D2Z, 1); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan1d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 2) - cufftResult result = cufftPlan2d(&fplan, fft_size[1], fft_size[0], CUFFT_D2Z); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan2d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 3) - cufftResult result = cufftPlan3d(&fplan, fft_size[2], fft_size[1], fft_size[0], CUFFT_D2Z); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan3d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#endif - -#else // host - -#if (AMREX_SPACEDIM == 1) - fplan = fftw_plan_dft_r2c_1d(fft_size[0], - phi_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 2) - fplan = fftw_plan_dft_r2c_2d(fft_size[1], fft_size[0], - phi_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 3) - fplan = fftw_plan_dft_r2c_3d(fft_size[2], fft_size[1], fft_size[0], - phi_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#endif - -#endif - - forward_plan.push_back(fplan); - } - - ParallelDescriptor::Barrier(); - - // ForwardTransform - for (MFIter mfi(phi_onegrid); mfi.isValid(); ++mfi) { - int i = mfi.LocalIndex(); -#ifdef AMREX_USE_CUDA - cufftSetStream(forward_plan[i], Gpu::gpuStream()); - cufftResult result = cufftExecD2Z(forward_plan[i], - phi_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field[i]->dataPtr())); - if (result != CUFFT_SUCCESS) { - AllPrint() << " forward transform using cufftExec failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#else - fftw_execute(forward_plan[i]); -#endif - } - - // copy data to a full-sized MultiFab - // this involves copying the complex conjugate from the half-sized field - // into the appropriate place in the full MultiFab - for (MFIter mfi(phi_dft_real_onegrid); mfi.isValid(); ++mfi) { - - Array4< GpuComplex > spectral = (*spectral_field[0]).array(); - - Array4 const& realpart = phi_dft_real_onegrid.array(mfi); - Array4 const& imagpart = phi_dft_imag_onegrid.array(mfi); - - Box bx = mfi.fabbox(); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - /* - Unpacking rules: - - For domains from (0,0,0) to (Nx-1,Ny-1,Nz-1) - - For any cells with i index > Nx/2, these values are complex conjugates of the corresponding - entry where (Nx-i,Ny-j,Nz-k) UNLESS that index is zero, in which case you use 0. - - e.g. for an 8^3 domain, any cell with i index - - Cell (6,2,3) is complex conjugate of (2,6,5) - - Cell (4,1,0) is complex conjugate of (4,7,0) (note that the FFT is computed for 0 <= i <= Nx/2) - */ - if (i <= bx.length(0)/2) { - // copy value - realpart(i,j,k) = spectral(i,j,k).real(); - imagpart(i,j,k) = spectral(i,j,k).imag(); - } else { - // copy complex conjugate - int iloc = bx.length(0)-i; - int jloc = 0; - int kloc = 0; -#if (AMREX_SPACEDIM >= 2) - jloc = (j == 0) ? 0 : bx.length(1)-j; -#endif -#if (AMREX_SPACEDIM == 3) - kloc = (k == 0) ? 0 : bx.length(2)-k; -#endif - - realpart(i,j,k) = spectral(iloc,jloc,kloc).real(); - imagpart(i,j,k) = -spectral(iloc,jloc,kloc).imag(); - } - - realpart(i,j,k) /= sqrtnpts; - imagpart(i,j,k) /= sqrtnpts; - }); - } - - // destroy fft plan - for (int i = 0; i < forward_plan.size(); ++i) { -#ifdef AMREX_USE_CUDA - cufftDestroy(forward_plan[i]); -#else - fftw_destroy_plan(forward_plan[i]); -#endif - } - - // ********************************** - // SHIFT DATA - // ********************************** - - // zero_avg=0 means set the k=0 value to zero, - // otherwise it sets the k=0 value to the average value of the signal in real space - int zero_avg = 1; - - // shift data - ShiftFFT(phi_dft_real_onegrid,geom,zero_avg); - ShiftFFT(phi_dft_imag_onegrid,geom,zero_avg); - - // ********************************** - // COPY DFT INTO THE DISTRIBUTED MULTIFAB - // ********************************** - - phi_dft_real.ParallelCopy(phi_dft_real_onegrid, 0, 0, 1); - phi_dft_imag.ParallelCopy(phi_dft_imag_onegrid, 0, 0, 1); - - // ********************************** - // WRITE DATA AND FFT TO PLOT FILE - // ********************************** - - // storage for magnitude and phase angle - MultiFab phi_dft_magn(ba, dm, 1, 0); - MultiFab phi_dft_phase(ba, dm, 1, 0); - - for (MFIter mfi(phi_dft_real); mfi.isValid(); ++mfi) - { - // Pointers to the magnitude, phase, real, and imaginary data - const Array4& phi_dft_magn_ptr = phi_dft_magn.array(mfi); - const Array4& phi_dft_phase_ptr = phi_dft_phase.array(mfi); - const Array4& phi_dft_real_ptr = phi_dft_real.array(mfi); - const Array4& phi_dft_imag_ptr = phi_dft_imag.array(mfi); - - const Box& bx = mfi.validbox(); - - // Set the value of the magnitude and phase angle using the real and imaginary parts of the dft - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - double re = phi_dft_real_ptr(i,j,k); - double im = phi_dft_imag_ptr(i,j,k); - phi_dft_magn_ptr(i,j,k) = std::sqrt(re*re + im*im); // Here we want to store the values of the magnitude - - // Avoid division by zero - if (re == 0.0) { - if (im == 0.0){ - phi_dft_phase_ptr(i,j,k) = 0.0; - } else if (im > 0.0) { - phi_dft_phase_ptr(i,j,k) = M_PI/2.0; - } else { - phi_dft_phase_ptr(i,j,k) = -M_PI/2.0; - } - } else { - phi_dft_phase_ptr(i,j,k) = std::atan(im/re); // Here we want to store the values of the phase angle - } - }); - } - - // storage for variables to write to plotfile - MultiFab plotfile(ba, dm, 5, 0); - - // copy phi, phi_dft_real, and phi_dft_imag, phi_dft_magn, and phi_dft_phase into plotfile - MultiFab::Copy(plotfile, phi , 0, 0, 1, 0); - MultiFab::Copy(plotfile, phi_dft_real, 0, 1, 1, 0); - MultiFab::Copy(plotfile, phi_dft_imag, 0, 2, 1, 0); - MultiFab::Copy(plotfile, phi_dft_magn, 0, 3, 1, 0); - MultiFab::Copy(plotfile, phi_dft_phase, 0, 4, 1, 0); - - // time and step are dummy variables required to WriteSingleLevelPlotfile - Real time = 0.; - int step = 0; - - // arguments - // 1: name of plotfile - // 2: MultiFab containing data to plot - // 3: variables names - // 4: geometry object - // 5: "time" of plotfile; not relevant in this example - // 6: "time step" of plotfile; not relevant in this example - WriteSingleLevelPlotfile("plt", plotfile, {"phi", "phi_dft_real", "phi_dft_imag","phi_dft_magn","phi_dft_phase"}, geom, time, step); - - // Call the timer again and compute the maximum difference between the start time - // and stop time over all processors - Real stop_time = ParallelDescriptor::second() - start_time; - ParallelDescriptor::ReduceRealMax(stop_time); - Print() << "Run time = " << stop_time << std::endl; - - } - Finalize(); -} - - diff --git a/ExampleCodes/FFTW/Basic/myfunc.H b/ExampleCodes/FFTW/Basic/myfunc.H deleted file mode 100644 index a70fe3c1..00000000 --- a/ExampleCodes/FFTW/Basic/myfunc.H +++ /dev/null @@ -1,18 +0,0 @@ -#ifndef MYFUNC_H_ -#define MYFUNC_H_ - -#include - -#ifdef AMREX_USE_CUDA -#include -#endif - -using namespace amrex; - -void ShiftFFT(MultiFab& dft_out, const Geometry& geom, const int& zero_avg); - -#ifdef AMREX_USE_CUDA -std::string cufftErrorToString (const cufftResult& err); -#endif - -#endif diff --git a/ExampleCodes/FFTW/Basic/myfunc.cpp b/ExampleCodes/FFTW/Basic/myfunc.cpp deleted file mode 100644 index 03e9ffbd..00000000 --- a/ExampleCodes/FFTW/Basic/myfunc.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include "myfunc.H" - -void ShiftFFT(MultiFab& dft_onegrid, const Geometry& geom, const int& zero_avg) { - - /* - Shifting rules: - - For each direction d, shift the data in the +d direction by n_cell[d]/2 and then modulo by n_cell[d]. - - e.g. for an 8^3 domain - - Cell (7,2,3) is shifted to ( (7+4)%8, (2+4)%8, (3+4)%8 ) = (3,6,7) - - */ - MultiFab dft_onegrid_temp; - dft_onegrid_temp.define(dft_onegrid.boxArray(), dft_onegrid.DistributionMap(), 1, 0); - - MultiFab::Copy(dft_onegrid_temp,dft_onegrid,0,0,1,0); - - // zero k=0 mode - if (zero_avg == 1) { - for (MFIter mfi(dft_onegrid); mfi.isValid(); ++mfi) { - - const Box& bx = mfi.tilebox(); - - const Array4& dft_temp = dft_onegrid_temp.array(mfi); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (i == 0 && j == 0 && k == 0) { - dft_temp(i,j,k) = 0.; - } - }); - } - } - - // Shift DFT by N/2+1 (pi) - for (MFIter mfi(dft_onegrid); mfi.isValid(); ++mfi) { - - const Box& bx = mfi.tilebox(); - - const Array4& dft = dft_onegrid.array(mfi); - const Array4& dft_temp = dft_onegrid_temp.array(mfi); - - int nx = bx.length(0); - int nxh = nx/2; -#if (AMREX_SPACEDIM >= 2) - int ny = bx.length(1); - int nyh = ny/2; -#endif -#if (AMREX_SPACEDIM == 3) - int nz = bx.length(2); - int nzh = nz/2; -#endif - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ip = (i+nxh)%nx; - int jp = 0; - int kp = 0; -#if (AMREX_SPACEDIM >= 2) - jp = (j+nyh)%ny; -#endif -#if (AMREX_SPACEDIM == 3) - kp = (k+nzh)%nz; -#endif - dft(ip,jp,kp) = dft_temp(i,j,k); - }); - - } - -} - -#ifdef AMREX_USE_CUDA -std::string cufftErrorToString (const cufftResult& err) -{ - switch (err) { - case CUFFT_SUCCESS: return "CUFFT_SUCCESS"; - case CUFFT_INVALID_PLAN: return "CUFFT_INVALID_PLAN"; - case CUFFT_ALLOC_FAILED: return "CUFFT_ALLOC_FAILED"; - case CUFFT_INVALID_TYPE: return "CUFFT_INVALID_TYPE"; - case CUFFT_INVALID_VALUE: return "CUFFT_INVALID_VALUE"; - case CUFFT_INTERNAL_ERROR: return "CUFFT_INTERNAL_ERROR"; - case CUFFT_EXEC_FAILED: return "CUFFT_EXEC_FAILED"; - case CUFFT_SETUP_FAILED: return "CUFFT_SETUP_FAILED"; - case CUFFT_INVALID_SIZE: return "CUFFT_INVALID_SIZE"; - case CUFFT_UNALIGNED_DATA: return "CUFFT_UNALIGNED_DATA"; - default: return std::to_string(err) + " (unknown error code)"; - } -} -#endif diff --git a/ExampleCodes/FFTW/Poisson/GNUmakefile b/ExampleCodes/FFTW/Poisson/GNUmakefile deleted file mode 100644 index e6680cc3..00000000 --- a/ExampleCodes/FFTW/Poisson/GNUmakefile +++ /dev/null @@ -1,22 +0,0 @@ -# AMREX_HOME defines the directory in which we will find all the AMReX code. -AMREX_HOME ?= ../../../../amrex - -USE_MPI = TRUE -DEBUG = FALSE -USE_OMP = FALSE -COMP = gnu -DIM = 2 - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package - -include $(AMREX_HOME)/Src/Base/Make.package - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules - -ifeq ($(USE_CUDA),TRUE) - LIBRARIES += -lcufft -else - LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3 -endif diff --git a/ExampleCodes/FFTW/Poisson/Make.package b/ExampleCodes/FFTW/Poisson/Make.package deleted file mode 100644 index 46fa6f45..00000000 --- a/ExampleCodes/FFTW/Poisson/Make.package +++ /dev/null @@ -1,3 +0,0 @@ -CEXE_sources += main.cpp myfunc.cpp -CEXE_headers = myfunc.H - diff --git a/ExampleCodes/FFTW/Poisson/inputs b/ExampleCodes/FFTW/Poisson/inputs deleted file mode 100644 index 0f403d40..00000000 --- a/ExampleCodes/FFTW/Poisson/inputs +++ /dev/null @@ -1,14 +0,0 @@ -amrex.the_arena_is_managed = 1 - -max_grid_size = 32 -n_cell_x = 64 -n_cell_y = 64 -n_cell_z = 64 - -prob_lo_x = 0. -prob_lo_y = 0. -prob_lo_z = 0. - -prob_hi_x = 1. -prob_hi_y = 1. -prob_hi_z = 1. diff --git a/ExampleCodes/FFTW/Poisson/main.cpp b/ExampleCodes/FFTW/Poisson/main.cpp deleted file mode 100644 index 6f6dbf01..00000000 --- a/ExampleCodes/FFTW/Poisson/main.cpp +++ /dev/null @@ -1,447 +0,0 @@ -#include - -#include -#include - -#ifdef AMREX_USE_CUDA -#include -#else -#include -#include -#endif - -#include "myfunc.H" - -using namespace amrex; - -int main (int argc, char* argv[]) -{ - - Initialize(argc,argv); - { - - // store the current time so we can later compute total run time. - Real start_time = ParallelDescriptor::second(); - - // ********************************** - // DECLARE SIMULATION PARAMETERS - // ********************************** - - // number of cells on each side of the domain - int n_cell_x; - int n_cell_y; - int n_cell_z; - - // dimensions of the domain - Real prob_lo_x; - Real prob_lo_y; - Real prob_lo_z; - - Real prob_hi_x; - Real prob_hi_y; - Real prob_hi_z; - - // This is the largest size a grid can be - int max_grid_size; - - // ********************************** - // READ PARAMETER VALUES FROM INPUTS FILE - // ********************************** - { - // ParmParse is way of reading inputs from the inputs file - // pp.get means we require the inputs file to have it - // pp.query means we optionally need the inputs file to have it - but you should supply a default value above - ParmParse pp; - - pp.get("n_cell_x",n_cell_x); - pp.get("n_cell_y",n_cell_y); - pp.get("n_cell_z",n_cell_z); - - pp.get("prob_lo_x",prob_lo_x); - pp.get("prob_lo_y",prob_lo_y); - pp.get("prob_lo_z",prob_lo_z); - - pp.get("prob_hi_x",prob_hi_x); - pp.get("prob_hi_y",prob_hi_y); - pp.get("prob_hi_z",prob_hi_z); - - pp.get("max_grid_size",max_grid_size); - } - - // Determine the domain length in each direction - Real L_x = std::abs(prob_hi_x - prob_lo_x); - Real L_y = std::abs(prob_hi_y - prob_lo_y); - Real L_z = std::abs(prob_hi_z - prob_lo_z); - - // ********************************** - // DEFINE SIMULATION SETUP AND GEOMETRY - // ********************************** - // make BoxArray and Geometry - // ba will contain a list of boxes that cover the domain - // geom contains information such as the physical domain size, - // number of points in the domain, and periodicity - BoxArray ba; - Geometry geom; - - // define lower and upper indices - IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); - IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1)); - - // Make a single box that is the entire domain - Box domain(dom_lo, dom_hi); - - // Initialize the boxarray "ba" from the single box "domain" - ba.define(domain); - - // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction - ba.maxSize(max_grid_size); - - // How Boxes are distrubuted among MPI processes - DistributionMapping dm(ba); - - // This defines the physical box size in each direction - RealBox real_box({ AMREX_D_DECL(prob_lo_x, prob_lo_y, prob_lo_z)}, - { AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} ); - - // periodic in all direction - Array is_periodic{AMREX_D_DECL(1,1,1)}; - - // This defines a Geometry object - geom.define(domain, real_box, CoordSys::cartesian, is_periodic); - - // extract dx from the geometry object - GpuArray dx = geom.CellSizeArray(); - - // MultiFab storage for rhs and the solution to the Poisson equation, lap(soln) = rhs - MultiFab rhs(ba, dm, 1, 0); - MultiFab soln(ba, dm, 1, 0); - - // ********************************** - // INITIALIZE RHS - // ********************************** - - // loop over boxes - for (MFIter mfi(rhs); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - - const Array4& rhs_ptr = rhs.array(mfi); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - // ********************************** - // SET VALUES FOR EACH CELL - // ********************************** - - Real x = (i+0.5) * dx[0]; - Real y = (AMREX_SPACEDIM>=2) ? (j+0.5) * dx[1] : 0.; - Real z = (AMREX_SPACEDIM==3) ? (k+0.5) * dx[2] : 0.; - - rhs_ptr(i,j,k) = std::exp(-10.*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))); - - }); - } - - // ********************************** - // COPY RHS INTO A MULTIFAB WITH ONE BOX - // ********************************** - - // create a new BoxArray and DistributionMapping for a MultiFab with 1 box - BoxArray ba_onegrid(geom.Domain()); - DistributionMapping dm_onegrid(ba_onegrid); - - // storage for rhs and soln - MultiFab rhs_onegrid(ba_onegrid, dm_onegrid, 1, 0); - MultiFab soln_onegrid(ba_onegrid, dm_onegrid, 1, 0); - - // copy rhs into rhs_onegrid - rhs_onegrid.ParallelCopy(rhs, 0, 0, 1); - - // ********************************** - // COMPUTE FFT - // ********************************** - -#ifdef AMREX_USE_CUDA - using FFTplan = cufftHandle; - using FFTcomplex = cuDoubleComplex; -#else - using FFTplan = fftw_plan; - using FFTcomplex = fftw_complex; -#endif - - // number of points in the domain - long npts = domain.numPts(); - Real sqrtnpts = std::sqrt(npts); - - // contain to store FFT - note it is shrunk by approximately a half in x - Vector > > > spectral_field; - - Vector forward_plan; - - for (MFIter mfi(rhs_onegrid); mfi.isValid(); ++mfi) { - - // grab a single box - Box realspace_bx = mfi.fabbox(); - - // size of box - IntVect fft_size = realspace_bx.length(); // This will be different for FFTs of complex data - - // this is the size of the box, except the 0th component is 'halved plus 1' - IntVect spectral_bx_size = fft_size; - spectral_bx_size[0] = fft_size[0]/2 + 1; - - // spectral box - Box spectral_bx = Box(IntVect(0), spectral_bx_size - IntVect(1)); - - spectral_field.emplace_back(new BaseFab >(spectral_bx,1, - The_Device_Arena())); - spectral_field.back()->setVal(0.0); // touch the memory - - FFTplan fplan; - -#ifdef AMREX_USE_CUDA - -#if (AMREX_SPACEDIM == 1) - cufftResult result = cufftPlan1d(&fplan, fft_size[0], CUFFT_D2Z, 1); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan1d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 2) - cufftResult result = cufftPlan2d(&fplan, fft_size[1], fft_size[0], CUFFT_D2Z); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan2d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 3) - cufftResult result = cufftPlan3d(&fplan, fft_size[2], fft_size[1], fft_size[0], CUFFT_D2Z); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan3d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#endif - -#else // host - -#if (AMREX_SPACEDIM == 1) - fplan = fftw_plan_dft_r2c_1d(fft_size[0], - rhs_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 2) - fplan = fftw_plan_dft_r2c_2d(fft_size[1], fft_size[0], - rhs_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 3) - fplan = fftw_plan_dft_r2c_3d(fft_size[2], fft_size[1], fft_size[0], - rhs_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field.back()->dataPtr()), - FFTW_ESTIMATE); -#endif - -#endif - - forward_plan.push_back(fplan); - } - - ParallelDescriptor::Barrier(); - - // ForwardTransform - for (MFIter mfi(rhs_onegrid); mfi.isValid(); ++mfi) { - int i = mfi.LocalIndex(); -#ifdef AMREX_USE_CUDA - cufftSetStream(forward_plan[i], Gpu::gpuStream()); - cufftResult result = cufftExecD2Z(forward_plan[i], - rhs_onegrid[mfi].dataPtr(), - reinterpret_cast - (spectral_field[i]->dataPtr())); - if (result != CUFFT_SUCCESS) { - AllPrint() << " forward transform using cufftExec failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#else - fftw_execute(forward_plan[i]); -#endif - } - - // Now we take the standard FFT and scale it by 1/k^2 - for (MFIter mfi(rhs_onegrid); mfi.isValid(); ++mfi) - { - Array4< GpuComplex > spectral = (*spectral_field[0]).array(); - - const Box& bx = mfi.fabbox(); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - // the "spectral" data only exists for i <= nx/2 - if (i <= bx.length(0)/2) { - - Real a = 2.*M_PI*i / L_x; - Real b = 2.*M_PI*j / L_y; - Real c = 2.*M_PI*k / L_z; - - // the values in the upper-half of the spectral array in y and z are here interpreted as negative wavenumbers - if (j >= n_cell_y/2) b = 2.*M_PI*(n_cell_y-j) / L_y; - if (k >= n_cell_z/2) c = 2.*M_PI*(n_cell_z-k) / L_z; - -#if (AMREX_SPACEDIM == 1) - Real k2 = -a*a; -#elif (AMREX_SPACEDIM == 2) - Real k2 = -(a*a + b*b); -#elif (AMREX_SPACEDIM == 3) - Real k2 = -(a*a + b*b + c*c); -#endif - - if (k2 != 0.) { - spectral(i,j,k) /= k2; - } else { - spectral(i,j,k) *= 0.; // interpretation here is that the average value of the solution is zero - } - } - - }); - } - - // Now we have completed the fft and scaled each value by 1/k^2 - // The scaled fft is inside spectral_field - // Take inverse fft of spectral_field and put it in soln_onegrid - Vector backward_plan; - - for (MFIter mfi(soln_onegrid); mfi.isValid(); ++mfi) { - - // grab a single box including ghost cell range - Box realspace_bx = mfi.fabbox(); - - // size of box including ghost cell range - IntVect fft_size = realspace_bx.length(); // This will be different for hybrid FFT - - FFTplan bplan; - -#ifdef AMREX_USE_CUDA - -#if (AMREX_SPACEDIM == 1) - cufftResult result = cufftPlan1d(&bplan, fft_size[0], CUFFT_Z2D, 1); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan1d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 2) - cufftResult result = cufftPlan2d(&bplan, fft_size[1], fft_size[0], CUFFT_Z2D); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan2d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#elif (AMREX_SPACEDIM == 3) - cufftResult result = cufftPlan3d(&bplan, fft_size[2], fft_size[1], fft_size[0], CUFFT_Z2D); - if (result != CUFFT_SUCCESS) { - AllPrint() << " cufftplan3d forward failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#endif - -#else // host - -#if (AMREX_SPACEDIM == 1) - bplan = fftw_plan_dft_c2r_1d(fft_size[0], - reinterpret_cast - (spectral_field.back()->dataPtr()), - soln_onegrid[mfi].dataPtr(), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 2) - bplan = fftw_plan_dft_c2r_2d(fft_size[1], fft_size[0], - reinterpret_cast - (spectral_field.back()->dataPtr()), - soln_onegrid[mfi].dataPtr(), - FFTW_ESTIMATE); -#elif (AMREX_SPACEDIM == 3) - bplan = fftw_plan_dft_c2r_3d(fft_size[2], fft_size[1], fft_size[0], - reinterpret_cast - (spectral_field.back()->dataPtr()), - soln_onegrid[mfi].dataPtr(), - FFTW_ESTIMATE); -#endif - -#endif - - backward_plan.push_back(bplan);// This adds an instance of bplan to the end of backward_plan - } - - for (MFIter mfi(soln_onegrid); mfi.isValid(); ++mfi) { - int i = mfi.LocalIndex(); -#ifdef AMREX_USE_CUDA - cufftSetStream(backward_plan[i], Gpu::gpuStream()); - cufftResult result = cufftExecZ2D(backward_plan[i], - reinterpret_cast - (spectral_field[i]->dataPtr()), - soln_onegrid[mfi].dataPtr()); - if (result != CUFFT_SUCCESS) { - AllPrint() << " inverse transform using cufftExec failed! Error: " - << cufftErrorToString(result) << "\n"; - } -#else - fftw_execute(backward_plan[i]); -#endif - - } - - // copy contents of soln_onegrid into soln - soln.ParallelCopy(soln_onegrid, 0, 0, 1); - - // Must divide each point by the total number of points in the domain for properly scaled inverse FFT - Real volinv = 1./npts; - soln.mult(volinv); - - // destroy fft plan - for (int i = 0; i < forward_plan.size(); ++i) { -#ifdef AMREX_USE_CUDA - cufftDestroy(forward_plan[i]); -#else - fftw_destroy_plan(forward_plan[i]); -#endif - } - - // destroy ifft plan - for (int i = 0; i < backward_plan.size(); ++i) { -#ifdef AMREX_USE_CUDA - cufftDestroy(backward_plan[i]); -#else - fftw_destroy_plan(backward_plan[i]); -#endif - - } - - // storage for variables to write to plotfile - MultiFab plotfile(ba, dm, 2, 0); - - // copy rhs and soln into plotfile - MultiFab::Copy(plotfile, rhs , 0, 0, 1, 0); - MultiFab::Copy(plotfile, soln, 0, 1, 1, 0); - - // time and step are dummy variables required to WriteSingleLevelPlotfile - Real time = 0.; - int step = 0; - - // arguments - // 1: name of plotfile - // 2: MultiFab containing data to plot - // 3: variables names - // 4: geometry object - // 5: "time" of plotfile; not relevant in this example - // 6: "time step" of plotfile; not relevant in this example - WriteSingleLevelPlotfile("plt", plotfile, {"rhs", "soln"}, geom, time, step); - - // Call the timer again and compute the maximum difference between the start time - // and stop time over all processors - Real stop_time = ParallelDescriptor::second() - start_time; - ParallelDescriptor::ReduceRealMax(stop_time); - Print() << "Run time = " << stop_time << std::endl; - - } - Finalize(); -} - - diff --git a/ExampleCodes/FFTW/Poisson/myfunc.H b/ExampleCodes/FFTW/Poisson/myfunc.H deleted file mode 100644 index b637f57d..00000000 --- a/ExampleCodes/FFTW/Poisson/myfunc.H +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef MYFUNC_H_ -#define MYFUNC_H_ - -#include - -#ifdef AMREX_USE_CUDA -#include -#endif - -using namespace amrex; - -#ifdef AMREX_USE_CUDA -std::string cufftErrorToString (const cufftResult& err); -#endif - -#endif diff --git a/ExampleCodes/FFTW/Poisson/myfunc.cpp b/ExampleCodes/FFTW/Poisson/myfunc.cpp deleted file mode 100644 index 7f98e091..00000000 --- a/ExampleCodes/FFTW/Poisson/myfunc.cpp +++ /dev/null @@ -1,20 +0,0 @@ -#include "myfunc.H" - -#ifdef AMREX_USE_CUDA -std::string cufftErrorToString (const cufftResult& err) -{ - switch (err) { - case CUFFT_SUCCESS: return "CUFFT_SUCCESS"; - case CUFFT_INVALID_PLAN: return "CUFFT_INVALID_PLAN"; - case CUFFT_ALLOC_FAILED: return "CUFFT_ALLOC_FAILED"; - case CUFFT_INVALID_TYPE: return "CUFFT_INVALID_TYPE"; - case CUFFT_INVALID_VALUE: return "CUFFT_INVALID_VALUE"; - case CUFFT_INTERNAL_ERROR: return "CUFFT_INTERNAL_ERROR"; - case CUFFT_EXEC_FAILED: return "CUFFT_EXEC_FAILED"; - case CUFFT_SETUP_FAILED: return "CUFFT_SETUP_FAILED"; - case CUFFT_INVALID_SIZE: return "CUFFT_INVALID_SIZE"; - case CUFFT_UNALIGNED_DATA: return "CUFFT_UNALIGNED_DATA"; - default: return std::to_string(err) + " (unknown error code)"; - } -} -#endif diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/GNUmakefile b/ExampleCodes/SWFFT/SWFFT_poisson/GNUmakefile deleted file mode 100644 index 068a3f11..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/GNUmakefile +++ /dev/null @@ -1,22 +0,0 @@ -AMREX_HOME ?= ../../../../amrex - -DIM = 3 - -USE_MPI = TRUE - -DEBUG = FALSE - -COMP = gcc - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(AMREX_HOME)/Src/Boundary/Make.package - -include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package -INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT -VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT - -LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3_omp -lfftw3 -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/Make.package b/ExampleCodes/SWFFT/SWFFT_poisson/Make.package deleted file mode 100644 index c2602745..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/Make.package +++ /dev/null @@ -1,7 +0,0 @@ - -CEXE_sources += main.cpp SWFFT_Test.cpp -CEXE_sources += swfft_solver.cpp -CEXE_headers += SWFFT_Test.H SWFFT_TestF.H - -F90EXE_sources += SWFFT_Test_F.F90 - diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/README b/ExampleCodes/SWFFT/SWFFT_poisson/README deleted file mode 100644 index 0f1200cb..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/README +++ /dev/null @@ -1,19 +0,0 @@ -This Tutorial demonstrates how to call the SWFFT wrapper to the FFTW3 solver. - -In this test case we set up a right hand side (rhs), call the forward transform, -modify the coefficients, then call the backward solver and output the solution -to the discrete Poisson equation. - -To build the code, type 'make' in this directory. - -To run the code, type 'main3d.gnu.MPI.ex inputs' in this directory - -To visualize the output, set the bool write_data to true, then -use amrvis3d (source available at https://github.com/AMReX-Codes/Amrvis): - -amrvis3d -mf RHS SOL_EXACT SOL_COMP - -to visualize the rhs, the exact solution and the computed solution. - -The max norm of the difference between the exact and computed solution is also printed. - diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.H b/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.H deleted file mode 100644 index 535cd7f2..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.H +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef MY_SWS_H_ -#define MY_SWS_H_ - -#include -#include - -// We solve Lap(soln) = rhs -class SWFFT_Test -{ -public: - SWFFT_Test (); - - void solve (); - -private: - - // Runtime parameters - amrex::IntVect max_grid_size; // Maximum allowable grid size in each coord direction - amrex::IntVect n_cell; // Number of cells in each coord direction - - int verbose = 2; - - // data - amrex::MultiFab rhs; - amrex::MultiFab soln; - amrex::MultiFab the_soln; - amrex::Geometry geom; - - void init_rhs (); - void comp_the_solution (); -}; - -#endif diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.cpp b/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.cpp deleted file mode 100644 index 45c3b7df..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test.cpp +++ /dev/null @@ -1,173 +0,0 @@ - -#include -#include - -#include -#include -#include - -// These are for SWFFT -#include -#include -#include - -#include - -#define ALIGN 16 - -using namespace amrex; - -void swfft_solver(MultiFab& rhs, MultiFab& soln, Geometry& geom, int verbose); - -SWFFT_Test::SWFFT_Test () -{ - static_assert(AMREX_SPACEDIM == 3, "3D only"); - - // runtime parameters - { - ParmParse pp; - - // Read in n_cell. Use defaults if not explicitly defined. - int cnt = pp.countval("n_cell"); - if (cnt > 1) { - Vector ncs; - pp.getarr("n_cell",ncs); - n_cell = IntVect{ncs[0],ncs[1],ncs[2]}; - } else if (cnt > 0) { - int ncs; - pp.get("n_cell",ncs); - n_cell = IntVect{ncs,ncs,ncs}; - } else { - n_cell = IntVect{32,32,32}; - } - - // Read in max_grid_size. Use defaults if not explicitly defined. - cnt = pp.countval("max_grid_size"); - if (cnt > 1) { - Vector mgs; - pp.getarr("max_grid_size",mgs); - max_grid_size = IntVect{mgs[0],mgs[1],mgs[2]}; - } else if (cnt > 0) { - int mgs; - pp.get("max_grid_size",mgs); - max_grid_size = IntVect{mgs,mgs,mgs}; - } else { - max_grid_size = IntVect{32,32,32}; - } - - pp.query("verbose", verbose); - } - - BoxArray ba; - { - // Make up a dx that is not 1 - Real dx = 1./double(n_cell[0]); - - IntVect dom_lo(0,0,0); - IntVect dom_hi(n_cell[0]-1,n_cell[1]-1,n_cell[2]-1); - Box domain(dom_lo,dom_hi); - ba.define(domain); - ba.maxSize(max_grid_size); - - // We assume a unit box (just for convenience) - Real x_hi = n_cell[0]*dx; - Real y_hi = n_cell[1]*dx; - Real z_hi = n_cell[2]*dx; - RealBox real_box({0.0,0.0,0.0}, {x_hi,y_hi,z_hi}); - - // The FFT assumes fully periodic boundaries - std::array is_periodic {1,1,1}; - - geom.define(domain, &real_box, CoordSys::cartesian, is_periodic.data()); - } - - // Make sure we define both the soln and the rhs with the same DistributionMapping - DistributionMapping dmap{ba}; - - // Note that we are defining rhs with NO ghost cells - rhs.define(ba, dmap, 1, 0); - init_rhs(); - - // Note that we are defining soln with NO ghost cells - soln.define(ba, dmap, 1, 0); - the_soln.define(ba, dmap, 1, 0); - - comp_the_solution(); -} - -void -SWFFT_Test::init_rhs () -{ - const Real* dx = geom.CellSize(); - Box domain(geom.Domain()); - - for (MFIter mfi(rhs,true); mfi.isValid(); ++mfi) - { - const Box& tbx = mfi.tilebox(); - fort_init_rhs(BL_TO_FORTRAN_BOX(tbx), - BL_TO_FORTRAN_ANYD(rhs[mfi]), - BL_TO_FORTRAN_BOX(domain), - geom.ProbLo(), geom.ProbHi(), dx); - } - - Real sum_rhs = rhs.sum(); - amrex::Print() << "Sum of rhs over the domain was " << sum_rhs << std::endl; - - sum_rhs = sum_rhs / domain.numPts(); - rhs.plus(-sum_rhs,0,1,0); - - sum_rhs = rhs.sum(); - amrex::Print() << "Sum of rhs over the domain is now " << sum_rhs << std::endl; -} - -void -SWFFT_Test::comp_the_solution () -{ - const Real* dx = geom.CellSize(); - Box domain(geom.Domain()); - - for (MFIter mfi(the_soln); mfi.isValid(); ++mfi) - { - const Box& tbx = mfi.tilebox(); - fort_comp_asol(BL_TO_FORTRAN_BOX(tbx), - BL_TO_FORTRAN_ANYD(the_soln[mfi]), - BL_TO_FORTRAN_BOX(domain), - geom.ProbLo(), geom.ProbHi(), dx); - } -} - -void -SWFFT_Test::solve () -{ - Real start_time = amrex::second(); - swfft_solver(rhs, soln, geom, verbose); - Real total_time = amrex::second() - start_time; - - bool write_data = false; - - if (write_data) - { - VisMF::Write(rhs,"RHS"); - VisMF::Write(soln,"SOL_COMP"); - VisMF::Write(the_soln,"SOL_EXACT"); - } - - if (verbose) - { - amrex::Print() << "MAX / MIN VALUE OF COMP SOLN " << soln.max(0) << " " - << soln.min(0) << std::endl; - amrex::Print() << "MAX / MIN VALUE OF EXACT SOLN " << the_soln.max(0) << " " - << the_soln.min(0) << std::endl; - } - - BoxArray ba = rhs.boxArray(); - const DistributionMapping& dm = rhs.DistributionMap(); - MultiFab diff(ba, dm, 1, 0); - MultiFab::Copy(diff, soln, 0, 0, 1, 0); - MultiFab::Subtract(diff, the_soln, 0, 0, 1, 0); - amrex::Print() << "\nMax-norm of the error is " << diff.norm0() << "\n"; - amrex::Print() << "Time spent in solve: " << total_time << std::endl; - - if (write_data) - VisMF::Write(diff,"DIFF"); -} diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.F90 b/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.F90 deleted file mode 100644 index 8db9c361..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.F90 +++ /dev/null @@ -1,87 +0,0 @@ -#include "AMReX_LO_BCTYPES.H" - -module abl_module - - use amrex_fort_module - use amrex_error_module - implicit none - -contains - - subroutine fort_init_rhs (lo, hi, rhs, rlo, rhi, domlo, domhi, problo, probhi, dx) & - bind(c,name="fort_init_rhs") - - integer, intent(in) :: lo(3), hi(3), rlo(3), rhi(3), domlo(3), domhi(3) - real(amrex_real), intent(inout) :: rhs(rlo(1):rhi(1),rlo(2):rhi(2),rlo(3):rhi(3)) - real(amrex_real), intent(in) :: problo(3), probhi(3) - real(amrex_real), intent(in) :: dx(3) - - integer :: i,j,k - real(amrex_real) :: x,y,z - real(amrex_real) :: Lx,Ly,Lz - real(amrex_real) :: pi, fpi, tpi, fac - - Lx = probhi(1)-problo(1) - Ly = probhi(2)-problo(2) - Lz = probhi(3)-problo(3) - - pi = 4.d0 * atan(1.d0) - tpi = 2.0d0 * pi - fpi = 4.0d0 * pi - fac = 3.0d0 * tpi**2 / (Lx**2 * Ly**2 * Lz**2) - - do k = lo(3), hi(3) - z = (dble(k)+0.5d0)*dx(3)/Lz - - do j = lo(2), hi(2) - y = (dble(j)+0.5d0)*dx(2)/Ly - - do i = lo(1), hi(1) - x = (dble(i)+0.5d0)*dx(1)/Lx - - rhs(i,j,k) = -fac * (sin(tpi*x) * sin(tpi*y) * sin(tpi*z)) & - & -fac * (sin(fpi*x) * sin(fpi*y) * sin(fpi*z)) - end do - end do - end do - - end subroutine fort_init_rhs - - subroutine fort_comp_asol (lo, hi, soln, slo, shi, domlo, domhi, problo, probhi, dx) & - bind(c,name='fort_comp_asol') - - integer, intent(in) :: lo(3), hi(3), slo(3), shi(3), domlo(3), domhi(3) - real(amrex_real), intent(inout) :: soln(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3)) - real(amrex_real), intent(in) :: problo(3), probhi(3) - real(amrex_real), intent(in) :: dx(3) - - integer :: i,j,k - real(amrex_real) :: x,y,z - real(amrex_real) :: Lx,Ly,Lz - real(amrex_real) :: pi, fpi, tpi - - pi = 4.d0 * atan(1.d0) - tpi = 2.0d0 * pi - fpi = 4.0d0 * pi - - Lx = probhi(1)-problo(1) - Ly = probhi(1)-problo(2) - Lz = probhi(3)-problo(3) - - do k = lo(3), hi(3) - z = (dble(k)+0.5d0)*dx(3) - do j = lo(2), hi(2) - y = (dble(j)+0.5d0)*dx(2) - do i = lo(1), hi(1) - x = (dble(i)+0.5d0)*dx(1) - - soln(i,j,k) = 1.d0 * (sin(tpi*x/Lx) * sin(tpi*y/Ly) * sin(tpi*z/Lz)) & - & + .25d0 * (sin(fpi*x/Lx) * sin(fpi*y/Ly) * sin(fpi*z/Lz)) - - end do - end do - end do - - end subroutine fort_comp_asol - -end module abl_module diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.H b/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.H deleted file mode 100644 index df5387f5..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/SWFFT_Test_F.H +++ /dev/null @@ -1,21 +0,0 @@ -#include - -#ifdef __cplusplus -extern "C" { -#endif - - void fort_init_rhs (const int* lo, const int* hi, - amrex_real* rhs, const int* rlo, const int* rhi, - const int* dom_lo, const int* dom_hi, - const amrex_real* problo, const amrex_real* probhi, - const amrex_real* dx); - - void fort_comp_asol (const int* lo, const int* hi, - const amrex_real* soln, const int* slo, const int* shi, - const int* dom_lo, const int* dom_hi, - const amrex_real* problo, const amrex_real* probhi, - const amrex_real* dx); - -#ifdef __cplusplus -} -#endif diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.128 b/ExampleCodes/SWFFT/SWFFT_poisson/inputs.128 deleted file mode 100644 index c1df9916..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.128 +++ /dev/null @@ -1,3 +0,0 @@ -n_cell = 128 128 128 -max_grid_size = 128 128 128 -verbose = 0 diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.32 b/ExampleCodes/SWFFT/SWFFT_poisson/inputs.32 deleted file mode 100644 index 309ba023..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.32 +++ /dev/null @@ -1,3 +0,0 @@ -n_cell = 32 32 32 -max_grid_size = 64 64 64 -verbose = 0 diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.64 b/ExampleCodes/SWFFT/SWFFT_poisson/inputs.64 deleted file mode 100644 index a9778979..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/inputs.64 +++ /dev/null @@ -1,3 +0,0 @@ -n_cell = 64 64 64 -max_grid_size = 64 64 64 -verbose = 0 diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/main.cpp b/ExampleCodes/SWFFT/SWFFT_poisson/main.cpp deleted file mode 100644 index 64828f10..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/main.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include -#include - -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - BL_PROFILE("main()"); - SWFFT_Test sw_test; - sw_test.solve(); - } - - amrex::Finalize(); - return 0; -} diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/run_me b/ExampleCodes/SWFFT/SWFFT_poisson/run_me deleted file mode 100644 index df674bcb..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/run_me +++ /dev/null @@ -1,3 +0,0 @@ -mpirun -n 4 main*ex inputs.32 -mpirun -n 4 main*ex inputs.64 -mpirun -n 4 main*ex inputs.128 diff --git a/ExampleCodes/SWFFT/SWFFT_poisson/swfft_solver.cpp b/ExampleCodes/SWFFT/SWFFT_poisson/swfft_solver.cpp deleted file mode 100644 index f4c81deb..00000000 --- a/ExampleCodes/SWFFT/SWFFT_poisson/swfft_solver.cpp +++ /dev/null @@ -1,169 +0,0 @@ - -#include -#include -#include - -// These are for SWFFT -#include -#include -#include - -#include - -#define ALIGN 16 - -using namespace amrex; - -void -swfft_solver(MultiFab& rhs, MultiFab& soln, Geometry& geom, int verbose) -{ - const BoxArray& ba = soln.boxArray(); - amrex::Print() << "BA " << ba << std::endl; - const DistributionMapping& dm = soln.DistributionMap(); - - if (rhs.nGrow() != 0 || soln.nGrow() != 0) - amrex::Error("Current implementation requires that both rhs and soln have no ghost cells"); - - // Define pi and (two pi) here - Real pi = 4 * std::atan(1.0); - Real tpi = 2 * pi; - - // We assume that all grids have the same size hence - // we have the same nx,ny,nz on all ranks - int nx = ba[0].size()[0]; - int ny = ba[0].size()[1]; - int nz = ba[0].size()[2]; - - Box domain(geom.Domain()); - - int nbx = domain.length(0) / nx; - int nby = domain.length(1) / ny; - int nbz = domain.length(2) / nz; - int nboxes = nbx * nby * nbz; - if (nboxes != ba.size()) - amrex::Error("NBOXES NOT COMPUTED CORRECTLY"); - amrex::Print() << "Number of boxes:\t" << nboxes << std::endl; - - Vector rank_mapping; - rank_mapping.resize(nboxes); - - DistributionMapping dmap = rhs.DistributionMap(); - - for (int ib = 0; ib < nboxes; ++ib) - { - int i = ba[ib].smallEnd(0) / nx; - int j = ba[ib].smallEnd(1) / ny; - int k = ba[ib].smallEnd(2) / nz; - - int local_index = k*nbx*nby + j*nbx + i; - - rank_mapping[local_index] = dmap[ib]; - if (verbose) - amrex::Print() << "LOADING RANK NUMBER " << dmap[ib] << " FOR GRID NUMBER " << ib - << " WHICH IS LOCAL NUMBER " << local_index << std::endl; - } - - - Real h = geom.CellSize(0); - Real hsq = h*h; - - Real start_time = amrex::second(); - - // Assume for now that nx = ny = nz - int Ndims[3] = { nbz, nby, nbx }; - int n[3] = {domain.length(2), domain.length(1), domain.length(0)}; - hacc::Distribution d(MPI_COMM_WORLD,n,Ndims,&rank_mapping[0]); - hacc::Dfft dfft(d); - - for (MFIter mfi(rhs,false); mfi.isValid(); ++mfi) - { - int gid = mfi.index(); - - size_t local_size = dfft.local_size(); - - std::vector > a; - std::vector > b; - - a.resize(nx*ny*nz); - b.resize(nx*ny*nz); - - dfft.makePlans(&a[0],&b[0],&a[0],&b[0]); - - // ******************************************* - // Copy real data from Rhs into real part of a -- no ghost cells and - // put into C++ ordering (not Fortran) - // ******************************************* - complex_t zero(0.0, 0.0); - size_t local_indx = 0; - for(size_t k=0; k<(size_t)nz; k++) { - for(size_t j=0; j<(size_t)ny; j++) { - for(size_t i=0; i<(size_t)nx; i++) { - - complex_t temp(rhs[mfi].dataPtr()[local_indx],0.); - a[local_indx] = temp; - local_indx++; - - } - } - } - -// ******************************************* -// Compute the forward transform -// ******************************************* - dfft.forward(&a[0]); - -// ******************************************* -// Now divide the coefficients of the transform -// ******************************************* - local_indx = 0; - const int *self = dfft.self_kspace(); - const int *local_ng = dfft.local_ng_kspace(); - const int *global_ng = dfft.global_ng(); - - for(size_t i=0; i<(size_t)local_ng[0]; i++) { - size_t global_i = local_ng[0]*self[0] + i; - - for(size_t j=0; j<(size_t)local_ng[1]; j++) { - size_t global_j = local_ng[1]*self[1] + j; - - for(size_t k=0; k<(size_t)local_ng[2]; k++) { - size_t global_k = local_ng[2]*self[2] + k; - - if (global_i == 0 && global_j == 0 & global_k == 0) { - a[local_indx] = 0; - } else { - double fac = 2. * ( - (cos(tpi*double(global_i)/double(global_ng[0])) - 1.) + - (cos(tpi*double(global_j)/double(global_ng[1])) - 1.) + - (cos(tpi*double(global_k)/double(global_ng[2])) - 1.) ); - - a[local_indx] = a[local_indx] / fac; - } - local_indx++; - - } - } - } - -// ******************************************* -// Compute the backward transform -// ******************************************* - dfft.backward(&a[0]); - - size_t global_size = dfft.global_size(); - double fac = hsq / global_size; - - local_indx = 0; - for(size_t k=0; k<(size_t)nz; k++) { - for(size_t j=0; j<(size_t)ny; j++) { - for(size_t i=0; i<(size_t)nx; i++) { - - // Divide by 2 pi N - soln[mfi].dataPtr()[local_indx] = fac * std::real(a[local_indx]); - local_indx++; - - } - } - } - } -} diff --git a/ExampleCodes/SWFFT/SWFFT_simple/GNUmakefile b/ExampleCodes/SWFFT/SWFFT_simple/GNUmakefile deleted file mode 100644 index 068a3f11..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/GNUmakefile +++ /dev/null @@ -1,22 +0,0 @@ -AMREX_HOME ?= ../../../../amrex - -DIM = 3 - -USE_MPI = TRUE - -DEBUG = FALSE - -COMP = gcc - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ./Make.package -include $(AMREX_HOME)/Src/Base/Make.package -include $(AMREX_HOME)/Src/Boundary/Make.package - -include $(AMREX_HOME)/Src/Extern/SWFFT/Make.package -INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT -VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/SWFFT - -LIBRARIES += -L$(FFTW_DIR) -lfftw3_mpi -lfftw3_omp -lfftw3 -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/SWFFT/SWFFT_simple/Make.package b/ExampleCodes/SWFFT/SWFFT_simple/Make.package deleted file mode 100644 index c924d73d..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/Make.package +++ /dev/null @@ -1,7 +0,0 @@ - -CEXE_sources += main.cpp SWFFT_Test.cpp -CEXE_sources += swfft_compute.cpp -CEXE_headers += SWFFT_Test.H SWFFT_TestF.H - -F90EXE_sources += SWFFT_Test_F.F90 - diff --git a/ExampleCodes/SWFFT/SWFFT_simple/README b/ExampleCodes/SWFFT/SWFFT_simple/README deleted file mode 100644 index a7be85b3..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/README +++ /dev/null @@ -1,36 +0,0 @@ -This Tutorial demonstrates how to call the SWFFT wrapper to the FFTW3 solver. - -In this test case we call the forward transform on an array of data. - -To build the code, type 'make' in this directory. - -To run the code, type 'main3d.gnu.MPI.ex inputs' in this directory - -To visualize the output, use amrvis3d (source available at -https://github.com/AMReX-Codes/Amrvis): - -amrvis3d plt_fft - -The DFT of the data as well as the original array can be viewed as separate -variables in the plotfile "plt_fft" - -Notes: this tutorial is very similar to amrex/Tutorials/SWFFT, -with these primary exceptions: - - This tutorial is simply performing an FFT on a multifab, unlike the SWFFT tutorial, - which solves a Poisson equation with periodic boundary conditions using - spectral methods - - Only a forward FFT is called. No reverse FFT is called, which is why the data - is redistributed from z-pencils back to blocks in k-space. In addition, - this means that grid ordering (rank_mapping) is not transposed when constructing the - hacc::Dfft class, to account for the difference between C++ and Fortran array indexing - - This tutorial works in 2D, by "hacking" hacc::Dfft, which is meant for data - arranged on a 3D domain - - The FFT of the data, as well as the original multifab, are printed to a single - plotfile - -More notes: - - The DFT is ordered in k-space from {0 to 2pi} (actualy {0 to 2pi*(N-1)/N}). In other words, - the kx, ky, and kz wave-numbers range from {0 to N-1}, as opposed to - {-(N/2-1) to N/2} or {-pi to pi} - - The FFT is not scaled - e.g. DFT of a discrete "delta" function (phi(0,0,0)=1, phi=0 - everywhere else) results in a constant value of phi_dft=1 throughout the domain \ No newline at end of file diff --git a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.H b/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.H deleted file mode 100644 index 87908050..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.H +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef MY_SWS_H_ -#define MY_SWS_H_ - -#include -#include - -// We solve Lap(phi_dft) = phi_spatial -class SWFFT_Test -{ -public: - SWFFT_Test (); - - void computeFFT (); - -private: - - // Runtime parameters - amrex::IntVect max_grid_size; // Maximum allowable grid size in each coord direction - amrex::IntVect n_cell; // Number of cells in each coord direction - - int verbose = 2; - int prob_type = 0; - - // data - amrex::MultiFab phi_spatial; - amrex::MultiFab phi_dft; - amrex::Geometry geom; - - void init_phi_spatial (); - void WritePlotFile(const int step=0, const amrex::Real time=0.0); -}; - -#endif diff --git a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.cpp b/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.cpp deleted file mode 100644 index 60d0872c..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test.cpp +++ /dev/null @@ -1,209 +0,0 @@ - -#include -#include - -#include -#include -#include - -#include - -// These are for SWFFT -#include -#include -#include - -#include - -#define ALIGN 16 - -using namespace amrex; - -void swfft_compute(MultiFab& phi_spatial, MultiFab& phi_dft, Geometry& geom, int verbose); - -SWFFT_Test::SWFFT_Test () -{ - - // runtime parameters - { - ParmParse pp; - - // Read in n_cell. Use defaults if not explicitly defined. - int cnt = pp.countval("n_cell"); - if (cnt > 1) { - Vector ncs; - pp.getarr("n_cell",ncs); -#if (AMREX_SPACEDIM == 2) - n_cell = IntVect{ncs[0],ncs[1]}; -#elif (AMREX_SPACEDIM == 3) - n_cell = IntVect{ncs[0],ncs[1],ncs[2]}; -#endif - } else if (cnt > 0) { - int ncs; - pp.get("n_cell",ncs); -#if (AMREX_SPACEDIM == 2) - n_cell = IntVect{ncs,ncs}; -#elif (AMREX_SPACEDIM == 3) - n_cell = IntVect{ncs,ncs,ncs}; -#endif - } else { -#if (AMREX_SPACEDIM == 2) - n_cell = IntVect{32,32}; -#elif (AMREX_SPACEDIM == 3) - n_cell = IntVect{32,32,32}; -#endif - } - - // Read in max_grid_size. Use defaults if not explicitly defined. - cnt = pp.countval("max_grid_size"); - if (cnt > 1) { - Vector mgs; - pp.getarr("max_grid_size",mgs); -#if (AMREX_SPACEDIM == 2) - max_grid_size = IntVect{mgs[0],mgs[1]}; -#elif (AMREX_SPACEDIM == 3) - max_grid_size = IntVect{mgs[0],mgs[1],mgs[2]}; -#endif - } else if (cnt > 0) { - int mgs; - pp.get("max_grid_size",mgs); -#if (AMREX_SPACEDIM == 2) - max_grid_size = IntVect{mgs,mgs}; -#elif (AMREX_SPACEDIM == 3) - max_grid_size = IntVect{mgs,mgs,mgs}; -#endif - } else { -#if (AMREX_SPACEDIM == 2) - max_grid_size = IntVect{32,32}; -#elif (AMREX_SPACEDIM == 3) - max_grid_size = IntVect{32,32,32}; -#endif - } - - pp.query("verbose", verbose); - pp.query("prob_type", prob_type); - } - - BoxArray ba; - { - // Make up a dx that is not 1 - Real dx = 1./double(n_cell[0]); - -#if (AMREX_SPACEDIM == 2) - IntVect dom_lo(0,0); - IntVect dom_hi(n_cell[0]-1,n_cell[1]-1); -#elif (AMREX_SPACEDIM == 3) - IntVect dom_lo(0,0,0); - IntVect dom_hi(n_cell[0]-1,n_cell[1]-1,n_cell[2]-1); -#endif - Box domain(dom_lo,dom_hi); - ba.define(domain); - ba.maxSize(max_grid_size); - - // We assume a unit box (just for convenience) - Real x_hi = n_cell[0]*dx; - Real y_hi = n_cell[1]*dx; -#if (AMREX_SPACEDIM == 2) - RealBox real_box({0.0,0.0}, {x_hi,y_hi}); -#elif (AMREX_SPACEDIM == 3) - Real z_hi = n_cell[2]*dx; - RealBox real_box({0.0,0.0,0.0}, {x_hi,y_hi,z_hi}); -#endif - - // The FFT assumes fully periodic boundaries -#if (AMREX_SPACEDIM == 2) - std::array is_periodic {1,1}; -#elif (AMREX_SPACEDIM == 3) - std::array is_periodic {1,1,1}; -#endif - geom.define(domain, &real_box, CoordSys::cartesian, is_periodic.data()); - } - - // Make sure we define both the phi_dft and the phi_spatial with the same DistributionMapping - DistributionMapping dmap{ba}; - - // Note that we are defining phi_spatial with NO ghost cells - phi_spatial.define(ba, dmap, 1, 0); - init_phi_spatial(); - - // Note that we are defining phi_dft with NO ghost cells - phi_dft.define(ba, dmap, 1, 0); - -} - -void -SWFFT_Test::init_phi_spatial () -{ - const Real* dx = geom.CellSize(); - Box domain(geom.Domain()); - - for (MFIter mfi(phi_spatial,true); mfi.isValid(); ++mfi) - { - const Box& tbx = mfi.tilebox(); - fort_init_phi_spatial(BL_TO_FORTRAN_BOX(tbx), - BL_TO_FORTRAN_ANYD(phi_spatial[mfi]), - BL_TO_FORTRAN_BOX(domain), - geom.ProbLo(), geom.ProbHi(), dx, &prob_type); - } -} - -void -SWFFT_Test::computeFFT () -{ - Real start_time = amrex::second(); - swfft_compute(phi_spatial, phi_dft, geom, verbose); - Real total_time = amrex::second() - start_time; - - WritePlotFile(); - bool write_data = false; - - if (write_data) - { - VisMF::Write(phi_spatial,"PHI_SPATIAL"); - VisMF::Write(phi_dft,"PHI_DFT"); - } - - if (verbose) - { - amrex::Print() << "MAX / MIN VALUE OF DFT " << phi_dft.max(0) << " " - << phi_dft.min(0) << std::endl; - } - - amrex::Print() << "Time spent taking FFT: " << total_time << std::endl; -} - -void -SWFFT_Test::WritePlotFile (const int step, const amrex::Real time) -{ - - MultiFab plotfile; - Vector varNames; - int nPlot; - - ////////////////////////////////////////////////////////////////////////////////// - // Write out FFT to plot file - ////////////////////////////////////////////////////////////////////////////////// - const std::string plotfilename = "plt_fft"; - nPlot = 2; - plotfile.define(phi_dft.boxArray(), phi_dft.DistributionMap(), nPlot, 0); - varNames.resize(nPlot); - - // keep a counter for plotfile variables - int cnt = 0; - - varNames[cnt++] = "FFT_of_phi"; - varNames[cnt++] = "phi"; - - // reset plotfile variable counter - cnt = 0; - - // copy into plotfile - MultiFab::Copy(plotfile, phi_dft, 0, cnt, 1, 0); - cnt++; - - MultiFab::Copy(plotfile, phi_spatial, 0, cnt, 1, 0); - cnt++; - - // write a plotfile - WriteSingleLevelPlotfile(plotfilename,plotfile,varNames,geom,time,step); -} diff --git a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.F90 b/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.F90 deleted file mode 100644 index 0d8022f7..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.F90 +++ /dev/null @@ -1,133 +0,0 @@ -#include "AMReX_LO_BCTYPES.H" - -module abl_module - - use amrex_fort_module - use amrex_error_module - implicit none - -contains - -#if (AMREX_SPACEDIM == 2) - - subroutine fort_init_phi_spatial (lo, hi, phi_spatial, rlo, rhi, domlo, domhi, problo, probhi, dx, prob_type) & - bind(c,name="fort_init_phi_spatial") - - integer, intent(in) :: prob_type - integer, intent(in) :: lo(2), hi(2), rlo(2), rhi(2), domlo(2), domhi(2) - real(amrex_real), intent(inout) :: phi_spatial(rlo(1):rhi(1),rlo(2):rhi(2)) - real(amrex_real), intent(in) :: problo(2), probhi(2) - real(amrex_real), intent(in) :: dx(2) - - integer :: i,j - real(amrex_real) :: x,y - real(amrex_real) :: Lx,Ly - real(amrex_real) :: pi, fpi, tpi, fac - - Lx = probhi(1)-problo(1) - Ly = probhi(2)-problo(2) - - pi = 4.d0 * atan(1.d0) - tpi = 2.0d0 * pi - fpi = 4.0d0 * pi - fac = 1.d0 - - do j = lo(2), hi(2) - y = (dble(j)+0.5d0)*dx(2)/Ly - - do i = lo(1), hi(1) - x = (dble(i)+0.5d0)*dx(1)/Lx - - SELECT CASE (prob_type) - CASE (0) - phi_spatial(i,j) = -fac * sin(tpi*x) - CASE (1) - phi_spatial(i,j) = -fac * (sin(tpi*x) * sin(tpi*y)) & - & -fac * (sin(fpi*x) * sin(fpi*y)) - CASE (2) - phi_spatial(i,j) = -fac * (sin(1.d0*tpi*x) * sin(1.d0*tpi*y)) & - & -fac * (sin(2.d0*tpi*x) * sin(2.d0*tpi*y)) & - & -fac * (sin(3.d0*tpi*x) * sin(3.d0*tpi*y)) & - & -fac * (sin(4.d0*tpi*x) * sin(4.d0*tpi*y)) - CASE (3) - phi_spatial = 0.d0 - if (SUM(lo-domlo) == 0) then - phi_spatial(lo(1),lo(2)) = 1.d0 - endif - CASE (4) - phi_spatial(i,j) = exp(-5.d2*((x/Lx-0.5d0)**2 + (y/Ly-0.5d0)**2)) - CASE DEFAULT - phi_spatial = 0.d0 - END SELECT - - end do - end do - - end subroutine fort_init_phi_spatial - -#elif (AMREX_SPACEDIM == 3) - - subroutine fort_init_phi_spatial (lo, hi, phi_spatial, rlo, rhi, domlo, domhi, problo, probhi, dx, prob_type) & - bind(c,name="fort_init_phi_spatial") - - integer, intent(in) :: prob_type - integer, intent(in) :: lo(3), hi(3), rlo(3), rhi(3), domlo(3), domhi(3) - real(amrex_real), intent(inout) :: phi_spatial(rlo(1):rhi(1),rlo(2):rhi(2),rlo(3):rhi(3)) - real(amrex_real), intent(in) :: problo(3), probhi(3) - real(amrex_real), intent(in) :: dx(3) - - integer :: i,j,k - real(amrex_real) :: x,y,z - real(amrex_real) :: Lx,Ly,Lz - real(amrex_real) :: pi, fpi, tpi, fac - - Lx = probhi(1)-problo(1) - Ly = probhi(2)-problo(2) - Lz = probhi(3)-problo(3) - - pi = 4.d0 * atan(1.d0) - tpi = 2.0d0 * pi - fpi = 4.0d0 * pi - ! fac = 3.0d0 * tpi**2 / (Lx**2 * Ly**2 * Lz**2) - fac = 1.0d0 - - do k = lo(3), hi(3) - z = (dble(k)+0.5d0)*dx(3)/Lz - - do j = lo(2), hi(2) - y = (dble(j)+0.5d0)*dx(2)/Ly - - do i = lo(1), hi(1) - x = (dble(i)+0.5d0)*dx(1)/Lx - - SELECT CASE (prob_type) - CASE (0) - phi_spatial(i,j,k) = -fac * sin(tpi*x) - CASE (1) - phi_spatial(i,j,k) = -fac * (sin(tpi*x) * sin(tpi*y) * sin(tpi*z)) & - & -fac * (sin(fpi*x) * sin(fpi*y) * sin(fpi*z)) - CASE (2) - phi_spatial(i,j,k) = -fac * (sin(1.d0*tpi*x) * sin(1.d0*tpi*y) * sin(1.d0*tpi*z)) & - & -fac * (sin(2.d0*tpi*x) * sin(2.d0*tpi*y) * sin(2.d0*tpi*z)) & - & -fac * (sin(3.d0*tpi*x) * sin(3.d0*tpi*y) * sin(3.d0*tpi*z)) & - & -fac * (sin(4.d0*tpi*x) * sin(4.d0*tpi*y) * sin(4.d0*tpi*z)) - CASE (3) - phi_spatial = 0.d0 - if (SUM(lo-domlo) == 0) then - phi_spatial(lo(1),lo(2),lo(3)) = 1.d0 - endif - CASE (4) - phi_spatial(i,j,k) = exp(-5.d2*((x/Lx-0.5d0)**2 + (y/Ly-0.5d0)**2 + (z/Lz-0.5d0)**2)) - CASE DEFAULT - phi_spatial = 0.d0 - END SELECT - - end do - end do - end do - - end subroutine fort_init_phi_spatial - -#endif - -end module abl_module diff --git a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.H b/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.H deleted file mode 100644 index 932ba485..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/SWFFT_Test_F.H +++ /dev/null @@ -1,21 +0,0 @@ -#include - -#ifdef __cplusplus -extern "C" { -#endif - - void fort_init_phi_spatial (const int* lo, const int* hi, - amrex_real* phi_spatial, const int* rlo, const int* rhi, - const int* dom_lo, const int* dom_hi, - const amrex_real* problo, const amrex_real* probhi, - const amrex_real* dx, int* prob_type); - - void fort_comp_asol (const int* lo, const int* hi, - const amrex_real* phi_dft, const int* slo, const int* shi, - const int* dom_lo, const int* dom_hi, - const amrex_real* problo, const amrex_real* probhi, - const amrex_real* dx); - -#ifdef __cplusplus -} -#endif diff --git a/ExampleCodes/SWFFT/SWFFT_simple/inputs.multipleGrids b/ExampleCodes/SWFFT/SWFFT_simple/inputs.multipleGrids deleted file mode 100644 index 2d339fff..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/inputs.multipleGrids +++ /dev/null @@ -1,13 +0,0 @@ -n_cell = 32 32 32 -max_grid_size = 16 16 16 -verbose = 0 - -/////////////////////////////////////// -// prob_type: -// prob_type = 0 - Fourier mode A -// prob_type = 1 - Fourier mode B -// prob_type = 2 - Fourier mode C -// prob_type = 3 - Delta -// prob_type = 4 - Gaussian -/////////////////////////////////////// -prob_type = 4 \ No newline at end of file diff --git a/ExampleCodes/SWFFT/SWFFT_simple/inputs.oneGrid b/ExampleCodes/SWFFT/SWFFT_simple/inputs.oneGrid deleted file mode 100644 index bbde9e6b..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/inputs.oneGrid +++ /dev/null @@ -1,13 +0,0 @@ -n_cell = 32 32 32 -max_grid_size = 32 32 32 -verbose = 0 - -/////////////////////////////////////// -// prob_type: -// prob_type = 0 - Fourier mode A -// prob_type = 1 - Fourier mode B -// prob_type = 2 - Fourier mode C -// prob_type = 3 - Delta -// prob_type = 4 - Gaussian -/////////////////////////////////////// -prob_type = 4 \ No newline at end of file diff --git a/ExampleCodes/SWFFT/SWFFT_simple/main.cpp b/ExampleCodes/SWFFT/SWFFT_simple/main.cpp deleted file mode 100644 index 3285fb5e..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/main.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include -#include - -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc,argv); - - { - BL_PROFILE("main()"); - SWFFT_Test sw_test; - sw_test.computeFFT(); - } - - amrex::Finalize(); - return 0; -} diff --git a/ExampleCodes/SWFFT/SWFFT_simple/run_me_2d b/ExampleCodes/SWFFT/SWFFT_simple/run_me_2d deleted file mode 100755 index 791ca4ad..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/run_me_2d +++ /dev/null @@ -1,2 +0,0 @@ -mpirun -n 1 ./main2*ex inputs.oneGrid -mpirun -n 4 ./main2*ex inputs.multipleGrids diff --git a/ExampleCodes/SWFFT/SWFFT_simple/run_me_3d b/ExampleCodes/SWFFT/SWFFT_simple/run_me_3d deleted file mode 100755 index e5383749..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/run_me_3d +++ /dev/null @@ -1,2 +0,0 @@ -mpirun -n 1 ./main3*ex inputs.oneGrid -mpirun -n 8 ./main3*ex inputs.multipleGrids diff --git a/ExampleCodes/SWFFT/SWFFT_simple/swfft_compute.cpp b/ExampleCodes/SWFFT/SWFFT_simple/swfft_compute.cpp deleted file mode 100644 index 948a4317..00000000 --- a/ExampleCodes/SWFFT/SWFFT_simple/swfft_compute.cpp +++ /dev/null @@ -1,162 +0,0 @@ - -#include -#include -#include -#include - -// These are for SWFFT -#include -#include -#include - -#include - -#define ALIGN 16 - -using namespace amrex; - -void -swfft_compute(MultiFab& phi_spatial, MultiFab& phi_dft, Geometry& geom, int verbose) -{ - const BoxArray& ba = phi_dft.boxArray(); - amrex::Print() << "BA " << ba << std::endl; - const DistributionMapping& dm = phi_dft.DistributionMap(); - - if (ba.size() != ParallelDescriptor::NProcs()) { - amrex::Error("Need same number of MPI processes as grids"); - exit(0); - } - - if (phi_spatial.nGrow() != 0 || phi_dft.nGrow() != 0) - amrex::Error("Current implementation requires that both phi_spatial and phi_dft have no ghost cells"); - - // We assume that all grids have the same size hence - // we have the same nx,ny,nz on all ranks - int nx = ba[0].size()[0]; - int ny = ba[0].size()[1]; -#if (AMREX_SPACEDIM == 2) - int nz = 1; -#elif (AMREX_SPACEDIM == 3) - int nz = ba[0].size()[2]; -#endif - - Box domain(geom.Domain()); - - int nbx = domain.length(0) / nx; - int nby = domain.length(1) / ny; -#if (AMREX_SPACEDIM == 2) - int nbz = 1; -#elif (AMREX_SPACEDIM == 3) - int nbz = domain.length(2) / nz; -#endif - int nboxes = nbx * nby * nbz; - if (nboxes != ba.size()) - amrex::Error("NBOXES NOT COMPUTED CORRECTLY"); - amrex::Print() << "Number of boxes:\t" << nboxes << std::endl; - - Vector rank_mapping; - rank_mapping.resize(nboxes); - - DistributionMapping dmap = phi_spatial.DistributionMap(); - - for (int ib = 0; ib < nboxes; ++ib) - { - int i = ba[ib].smallEnd(0) / nx; - int j = ba[ib].smallEnd(1) / ny; -#if (AMREX_SPACEDIM == 2) - int k = 0; -#elif (AMREX_SPACEDIM == 3) - int k = ba[ib].smallEnd(2) / nz; -#endif - - // This would be the "correct" local index if the data wasn't being transformed - int local_index = k*nbx*nby + j*nbx + i; - - // This is what we pass to dfft to compensate for the Fortran ordering - // of amrex data in MultiFabs. - // int local_index = i*nby*nbz + j*nbz + k; - - rank_mapping[local_index] = dmap[ib]; - - if (verbose) - amrex::Print() << "LOADING RANK NUMBER " << dmap[ib] << " FOR GRID NUMBER " << ib - << " WHICH IS LOCAL NUMBER " << local_index << std::endl; - } - - Real start_time = amrex::second(); - - // Assume for now that nx = ny = nz -#if (AMREX_SPACEDIM == 2) - int Ndims[3] = { 1, nby, nbx}; - int n[3] = { 1, domain.length(1), domain.length(0)}; -#elif (AMREX_SPACEDIM == 3) - int Ndims[3] = { nbz, nby, nbx }; - int n[3] = { domain.length(2), domain.length(1), domain.length(0)}; -#endif - hacc::Distribution d(MPI_COMM_WORLD,n,Ndims,&rank_mapping[0]); - hacc::Dfft dfft(d); - - for (MFIter mfi(phi_spatial,false); mfi.isValid(); ++mfi) - { - int gid = mfi.index(); - - size_t local_size = dfft.local_size(); - - std::vector > a; - std::vector > b; - - a.resize(nx*ny*nz); - b.resize(nx*ny*nz); - - dfft.makePlans(&a[0],&b[0],&a[0],&b[0]); - - // ******************************************* - // Copy real data from Rhs into real part of a -- no ghost cells and - // put into C++ ordering (not Fortran) - // ******************************************* - complex_t zero(0.0, 0.0); - size_t local_indx = 0; - for(size_t k=0; k<(size_t)nz; k++) { - for(size_t j=0; j<(size_t)ny; j++) { - for(size_t i=0; i<(size_t)nx; i++) { - - complex_t temp(phi_spatial[mfi].dataPtr()[local_indx],0.); - a[local_indx] = temp; - local_indx++; - - } - } - } - -// ******************************************* -// Compute the forward transform -// ******************************************* - dfft.forward(&a[0]); - -// ******************************************* -// Redistribute from z-pencils to blocks -// ******************************************* - d.redistribute_2_to_3(&a[0],&b[0],2); - - size_t global_size = dfft.global_size(); - double fac; - - // fac = sqrt(1.0 / (double)global_size); - fac = 1.0; // Overwrite fac - - local_indx = 0; - // int local_indx_p = 0; - for(size_t k=0; k<(size_t)nz; k++) { - for(size_t j=0; j<(size_t)ny; j++) { - for(size_t i=0; i<(size_t)nx; i++) { - - // Divide by 2 pi N - phi_dft[mfi].dataPtr()[local_indx] = fac * std::abs(b[local_indx]); - local_indx++; - - } - } - } - } - -} diff --git a/ExampleCodes/heFFTe/Basic/GNUmakefile b/ExampleCodes/heFFTe/Basic/GNUmakefile deleted file mode 100644 index 26827791..00000000 --- a/ExampleCodes/heFFTe/Basic/GNUmakefile +++ /dev/null @@ -1,35 +0,0 @@ -AMREX_HOME ?= ../../../../amrex -HEFFTE_HOME ?= ../../../../heffte/build - -DEBUG = FALSE -DIM = 3 -COMP = gcc -TINY_PROFILE = TRUE -USE_MPI = TRUE -USE_CUDA = FALSE -USE_HIP = FALSE - -BL_NO_FORT = TRUE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs -include $(AMREX_HOME)/Src/Base/Make.package -include Make.package - -VPATH_LOCATIONS += $(HEFFTE_HOME)/include -INCLUDE_LOCATIONS += $(HEFFTE_HOME)/include -LIBRARY_LOCATIONS += $(HEFFTE_HOME)/lib - -libraries += -lheffte - -ifeq ($(USE_CUDA),TRUE) - libraries += -lcufft -else ifeq ($(USE_HIP),TRUE) - # Use rocFFT. ROC_PATH is defined in amrex - INCLUDE_LOCATIONS += $(ROC_PATH)/rocfft/include - LIBRARY_LOCATIONS += $(ROC_PATH)/rocfft/lib - LIBRARIES += -L$(ROC_PATH)/rocfft/lib -lrocfft -else - libraries += -lfftw3_mpi -lfftw3f -lfftw3 -endif - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/heFFTe/Basic/Make.package b/ExampleCodes/heFFTe/Basic/Make.package deleted file mode 100644 index 6b4b865e..00000000 --- a/ExampleCodes/heFFTe/Basic/Make.package +++ /dev/null @@ -1 +0,0 @@ -CEXE_sources += main.cpp diff --git a/ExampleCodes/heFFTe/Basic/README b/ExampleCodes/heFFTe/Basic/README deleted file mode 100644 index 4ec6cdae..00000000 --- a/ExampleCodes/heFFTe/Basic/README +++ /dev/null @@ -1,41 +0,0 @@ -heFFTe instaallation guide: -https://icl-utk-edu.github.io/heffte/md_doxygen_installation.html - -Installation - ->> git clone git@github.com:icl-utk-edu/heffte.git -OR ->> git clone https://github.com/icl-utk-edu/heffte.git - ->> cd heffte ->> mkdir build ->> cd build - -###################### -HOST BUILD -###################### - -# NOTE: -DCMAKE_INSTALL_PREFIX can be a different location /path/to/DCMAKE_INSTALL_PREFIX - ->> cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=OFF .. - ->> make -j4 ->> make install - ->> cd /path/to/amrex-tutorials/ExampleCodes/heFFTe/Basic ->> make -j4 HEFFTE_HOME=/path/to/DCMAKE_INSTALL_PREFIX - -###################### -NVIDIA/CUDA BUILD -###################### - -# NOTE: -DCMAKE_INSTALL_PREFIX can be a different location /path/to/DCMAKE_INSTALL_PREFIX -# NOTE: If you configure with -DHeffte_DISABLE_GPU_AWARE_MPI=OFF, you must use --gpu-bind=none in your slurm script - ->> cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=ON -DHeffte_DISABLE_GPU_AWARE_MPI=OFF .. - ->> make -j4 ->> make install - ->> cd /path/to/amrex-tutorials/ExampleCodes/heFFTe/Basic/ ->> make -j4 USE_CUDA=TRUE HEFFTE_HOME=/path/to/DCMAKE_INSTALL_PREFIX diff --git a/ExampleCodes/heFFTe/Basic/inputs b/ExampleCodes/heFFTe/Basic/inputs deleted file mode 100644 index 0e81402b..00000000 --- a/ExampleCodes/heFFTe/Basic/inputs +++ /dev/null @@ -1,15 +0,0 @@ -n_cell_x = 64 -n_cell_y = 64 -n_cell_z = 64 - -max_grid_size_x = 32 -max_grid_size_y = 32 -max_grid_size_z = 64 - -prob_lo_x = 0. -prob_lo_y = 0. -prob_lo_z = 0. - -prob_hi_x = 1. -prob_hi_y = 1. -prob_hi_z = 1. diff --git a/ExampleCodes/heFFTe/Basic/main.cpp b/ExampleCodes/heFFTe/Basic/main.cpp deleted file mode 100644 index acb44cde..00000000 --- a/ExampleCodes/heFFTe/Basic/main.cpp +++ /dev/null @@ -1,462 +0,0 @@ -#include - -#include -#include -#include -#include -#include - -using namespace amrex; - -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc, argv); { - - BL_PROFILE("main"); - - // ********************************** - // DECLARE SIMULATION PARAMETERS - // ********************************** - - // number of cells on each side of the domain - int n_cell_x; - int n_cell_y; - int n_cell_z; - - // maximum grid size in each direction - int max_grid_size_x; - int max_grid_size_y; - int max_grid_size_z; - - // physical dimensions of the domain - Real prob_lo_x = 0.; - Real prob_lo_y = 0.; - Real prob_lo_z = 0.; - Real prob_hi_x = 1.; - Real prob_hi_y = 1.; - Real prob_hi_z = 1.; - - // ********************************** - // READ PARAMETER VALUES FROM INPUTS FILE - // ********************************** - { - // ParmParse is way of reading inputs from the inputs file - // pp.get means we require the inputs file to have it - // pp.query means we optionally need the inputs file to have it - but you should supply a default value above - - ParmParse pp; - - pp.get("n_cell_x",n_cell_x); - pp.get("n_cell_y",n_cell_y); - pp.get("n_cell_z",n_cell_z); - - pp.get("max_grid_size_x",max_grid_size_x); - pp.get("max_grid_size_y",max_grid_size_y); - pp.get("max_grid_size_z",max_grid_size_z); - - pp.query("prob_lo_x",prob_lo_x); - pp.query("prob_lo_y",prob_lo_y); - pp.query("prob_lo_z",prob_lo_z); - - pp.query("prob_hi_x",prob_hi_x); - pp.query("prob_hi_y",prob_hi_y); - pp.query("prob_hi_z",prob_hi_z); - } - - - // define lower and upper indices of domain - IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); - IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1)); - - // Make a single box that is the entire domain - Box domain(dom_lo, dom_hi); - - // number of points in the domain - long npts = domain.numPts(); - Real sqrtnpts = std::sqrt(npts); - - // Initialize the boxarray "ba" from the single box "domain" - BoxArray ba(domain); - - // create IntVect of max_grid_size - IntVect max_grid_size(AMREX_D_DECL(max_grid_size_x,max_grid_size_y,max_grid_size_z)); - - // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction - ba.maxSize(max_grid_size); - - // How Boxes are distrubuted among MPI processes - DistributionMapping dm(ba); - - // This defines the physical box size in each direction - RealBox real_box({ AMREX_D_DECL(prob_lo_x, prob_lo_y, prob_lo_z)}, - { AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} ); - - // periodic in all direction - Array is_periodic{AMREX_D_DECL(1,1,1)}; - - // geometry object for real data - Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic); - - // extract dx from the geometry object - GpuArray dx = geom.CellSizeArray(); - - MultiFab phi(ba,dm,1,0); - - // check to make sure each MPI rank has exactly 1 box - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(phi.local_size() == 1, "Must have one Box per MPI process"); - - Real omega = M_PI/2.0; - - for (MFIter mfi(phi); mfi.isValid(); ++mfi) { - - Array4 const& fab = phi.array(mfi); - - const Box& bx = mfi.fabbox(); - - amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, RandomEngine const& engine) noexcept - { - - // ********************************** - // SET VALUES FOR EACH CELL - // ********************************** - - Real x = prob_lo_x + (i+0.5) * dx[0]; - Real y = (AMREX_SPACEDIM>=2) ? prob_lo_y + (j+0.5) * dx[1] : 0.; - Real z = (AMREX_SPACEDIM==3) ? prob_lo_z + (k+0.5) * dx[2] : 0.; - fab(i,j,k) = std::sin(4*M_PI*x/prob_hi_x + omega); - if (AMREX_SPACEDIM >= 2) { - fab(i,j,k) *= std::sin(6*M_PI*y/prob_hi_y + omega); - } - if (AMREX_SPACEDIM == 3) { - fab(i,j,k) *= std::sin(8*M_PI*z/prob_hi_z + omega); - } - - // fab(i,j,k) = std::sqrt(x)*std::sqrt(y); - - // fab(i,j,k) = amrex::Random(engine); - - }); - } - - Real time = 0.; - int step = 0; - - // write out phi to plotfile - WriteSingleLevelPlotfile("phi", phi, {"phi"}, geom, time, step); - - // since there is 1 MPI rank per box, here each MPI rank obtains its local box and the associated boxid - Box local_box; - int local_boxid; - { - for (int i = 0; i < ba.size(); ++i) { - Box b = ba[i]; - // each MPI rank has its own local_box Box and local_boxid ID - if (ParallelDescriptor::MyProc() == dm[i]) { - local_box = b; - local_boxid = i; - } - } - } - - // now each MPI rank works on its own box - // for real->complex fft's, the fft is stored in an (nx/2+1) x ny x nz dataset - - // start by coarsening each box by 2 in the x-direction - Box c_local_box = amrex::coarsen(local_box, IntVect(AMREX_D_DECL(2,1,1))); - - // if the coarsened box's high-x index is even, we shrink the size in 1 in x - // this avoids overlap between coarsened boxes - if (c_local_box.bigEnd(0) * 2 == local_box.bigEnd(0)) { - c_local_box.setBig(0,c_local_box.bigEnd(0)-1); - } - // for any boxes that touch the hi-x domain we - // increase the size of boxes by 1 in x - // this makes the overall fft dataset have size (Nx/2+1 x Ny x Nz) - if (local_box.bigEnd(0) == geom.Domain().bigEnd(0)) { - c_local_box.growHi(0,1); - } - - // each MPI rank gets storage for its piece of the fft - BaseFab > spectral_field(c_local_box, 1, The_Device_Arena()); - - // create real->complex fft objects with the appropriate backend and data about - // the domain size and its local box size -#if (AMREX_SPACEDIM==2) - -#ifdef AMREX_USE_CUDA - heffte::fft2d_r2c fft -#elif AMREX_USE_HIP - heffte::fft2d_r2c fft -#else - heffte::fft2d_r2c fft -#endif - ({{local_box.smallEnd(0),local_box.smallEnd(1),0}, - {local_box.bigEnd(0) ,local_box.bigEnd(1) ,0}}, - {{c_local_box.smallEnd(0),c_local_box.smallEnd(1),0}, - {c_local_box.bigEnd(0) ,c_local_box.bigEnd(1) ,0}}, - 0, ParallelDescriptor::Communicator()); - -#elif (AMREX_SPACEDIM==3) - -#ifdef AMREX_USE_CUDA - heffte::fft3d_r2c fft -#elif AMREX_USE_HIP - heffte::fft3d_r2c fft -#else - heffte::fft3d_r2c fft -#endif - ({{local_box.smallEnd(0),local_box.smallEnd(1),local_box.smallEnd(2)}, - {local_box.bigEnd(0) ,local_box.bigEnd(1) ,local_box.bigEnd(2)}}, - {{c_local_box.smallEnd(0),c_local_box.smallEnd(1),c_local_box.smallEnd(2)}, - {c_local_box.bigEnd(0) ,c_local_box.bigEnd(1) ,c_local_box.bigEnd(2)}}, - 0, ParallelDescriptor::Communicator()); - -#endif - - using heffte_complex = typename heffte::fft_output::type; - heffte_complex* spectral_data = (heffte_complex*) spectral_field.dataPtr(); - - { BL_PROFILE("HEFFTE-total"); - { - BL_PROFILE("ForwardTransform"); - fft.forward(phi[local_boxid].dataPtr(), spectral_data); - } - { - BL_PROFILE("BackwardTransform"); - fft.backward(spectral_data, phi[local_boxid].dataPtr()); - } - } - - // scale by 1/npts (both forward and inverse need 1/sqrtnpts scaling so I am doing it all here) - phi.mult(1./npts); - - // ********************************** - // diagnostics - // ********************************** - - // create a BoxArray containing the fft boxes - // by construction, these boxes correlate to the associated spectral_data - // this we can copy the spectral data into this multifab since we know they are owned by the same MPI rank - BoxArray fft_ba; - { - BoxList bl; - bl.reserve(ba.size()); - - for (int i = 0; i < ba.size(); ++i) { - Box b = ba[i]; - - Box r_box = b; - Box c_box = amrex::coarsen(r_box, IntVect(AMREX_D_DECL(2,1,1))); - - // this avoids overlap for the cases when one or more r_box's - // have an even cell index in the hi-x cell - if (c_box.bigEnd(0) * 2 == r_box.bigEnd(0)) { - c_box.setBig(0,c_box.bigEnd(0)-1); - } - - // increase the size of boxes touching the hi-x domain by 1 in x - // this is an (Nx x Ny x Nz) -> (Nx/2+1 x Ny x Nz) real-to-complex sizing - if (b.bigEnd(0) == geom.Domain().bigEnd(0)) { - c_box.growHi(0,1); - } - bl.push_back(c_box); - - } - fft_ba.define(std::move(bl)); - } - - // storage for real, imaginary, magnitude, and phase - MultiFab fft_data(fft_ba,dm,4,0); - - // this copies the spectral data into a distributed MultiFab - for (MFIter mfi(fft_data); mfi.isValid(); ++mfi) { - - Array4 const& data = fft_data.array(mfi); - Array4< GpuComplex > spectral = spectral_field.array(); - - const Box& bx = mfi.fabbox(); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real re = spectral(i,j,k).real() / sqrtnpts; - Real im = spectral(i,j,k).imag() / sqrtnpts; - - data(i,j,k,0) = re; - data(i,j,k,1) = im; - data(i,j,k,2) = std::sqrt(re*re + im*im); - - // Here we want to store the values of the phase angle - // Avoid division by zero - if (re == 0.0) { - if (im == 0.0){ - data(i,j,k,3) = 0.0; - } else if (im > 0.0) { - data(i,j,k,3) = M_PI/2.0; - } else { - data(i,j,k,3) = -M_PI/2.0; - } - } else { - data(i,j,k,3) = std::atan(im/re); - } - }); - } - - // domain for fft data used to contruct a geometry object - Box domain_fft = amrex::coarsen(domain, IntVect(AMREX_D_DECL(2,1,1))); - // shrink by 1 in x in case there are an odd number of cells in the x-direction in domain - if (domain_fft.bigEnd(0) * 2 == domain.bigEnd(0)) { - domain_fft.setBig(0,domain_fft.bigEnd(0)-1); - } - // grow by 1 in the x-direction to match the size of the FFT - domain_fft.growHi(0,1); - - Geometry geom_fft(domain_fft, real_box, CoordSys::cartesian, is_periodic); - - WriteSingleLevelPlotfile("fft_data", fft_data, {"real", "imag", "magitude", "phase"}, geom_fft, time, step); - - // ********************************** - // unpack data onto a 'full'-sized MultiFab - // shift data so k=0 mode is at the center - // ********************************** - - BoxArray ba_onegrid(domain); - DistributionMapping dm_onegrid(ba_onegrid); - - // real, imaginary, magnitude, phase - MultiFab fft_data_onegrid(ba_onegrid, dm_onegrid, 4, 0); - MultiFab fft_data_onegrid_shifted(ba_onegrid, dm_onegrid, 4, 0); - - // copy in real and imaginary parts - fft_data_onegrid.ParallelCopy(fft_data, 0, 0, 2); - - // unpack data into a 'full'-sized MultiFab - // this involves copying the complex conjugate from the half-sized field - // into the appropriate place in the full MultiFab - for (MFIter mfi(fft_data_onegrid); mfi.isValid(); ++mfi) { - - Array4 const& data = fft_data_onegrid.array(mfi); - - const Box& bx = mfi.fabbox(); - - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - /* - Unpacking rules: - - For domains from (0,0,0) to (Nx-1,Ny-1,Nz-1) - - For any cells with i index > Nx/2, these values are complex conjugates of the corresponding - entry where (Nx-i,Ny-j,Nz-k) UNLESS that index is zero, in which case you use 0. - - e.g. for an 8^3 domain, any cell with i index - - Cell (6,2,3) is complex conjugate of (2,6,5) - - Cell (4,1,0) is complex conjugate of (4,7,0) (note that the FFT is computed for 0 <= i <= Nx/2) - */ - if (i <= bx.length(0)/2) { - // do nothing, data is already here - } else { - // copy complex conjugate - int iloc = bx.length(0)-i; - int jloc = 0; - int kloc = 0; -#if (AMREX_SPACEDIM >= 2) - jloc = (j == 0) ? 0 : bx.length(1)-j; -#endif -#if (AMREX_SPACEDIM == 3) - kloc = (k == 0) ? 0 : bx.length(2)-k; -#endif - - data(i,j,k,0) = data(iloc,jloc,kloc,0); - data(i,j,k,1) = -data(iloc,jloc,kloc,1); - } - - Real re = data(i,j,k,0); - Real im = data(i,j,k,1); - - data(i,j,k,2) = std::sqrt(re*re + im*im); - - // Here we want to store the values of the phase angle - // Avoid division by zero - if (re == 0.0) { - if (im == 0.0){ - data(i,j,k,3) = 0.0; - } else if (im > 0.0) { - data(i,j,k,3) = M_PI/2.0; - } else { - data(i,j,k,3) = -M_PI/2.0; - } - } else { - data(i,j,k,3) = std::atan(im/re); - } - - }); - } - - // zero_avg=0 means set the k=0 value to zero, - // otherwise it sets the k=0 value to the average value of the signal in real space - int zero_avg = 1; - - // zero k=0 mode - if (zero_avg == 1) { - for (MFIter mfi(fft_data_onegrid); mfi.isValid(); ++mfi) { - - const Box& bx = mfi.tilebox(); - - Array4 const& data = fft_data_onegrid.array(mfi); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (i == 0 && j == 0 && k == 0) { - for (int comp=0; comp<4; ++comp) { - data(i,j,k,comp) = 0.; - } - } - }); - } - } - - // SHIFT DATA - for (MFIter mfi(fft_data_onegrid); mfi.isValid(); ++mfi) { - - const Box& bx = mfi.tilebox(); - - const Array4& dft = fft_data_onegrid_shifted.array(mfi); - const Array4& dft_temp = fft_data_onegrid.array(mfi); - - int nx = bx.length(0); - int nxh = nx/2; -#if (AMREX_SPACEDIM >= 2) - int ny = bx.length(1); - int nyh = ny/2; -#endif -#if (AMREX_SPACEDIM == 3) - int nz = bx.length(2); - int nzh = nz/2; -#endif - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - int ip = (i+nxh)%nx; - int jp = 0; - int kp = 0; -#if (AMREX_SPACEDIM >= 2) - jp = (j+nyh)%ny; -#endif -#if (AMREX_SPACEDIM == 3) - kp = (k+nzh)%nz; -#endif - for (int comp=0; comp<4; ++comp) { - dft(ip,jp,kp,comp) = dft_temp(i,j,k,comp); - } - }); - - } - - WriteSingleLevelPlotfile("plt", fft_data_onegrid_shifted, - {"phi_dft_real", "phi_dft_imag", "phi_dft_magitude", "phi_dft_phase"}, - geom, time, step); - - } amrex::Finalize(); -} diff --git a/ExampleCodes/heFFTe/Poisson/GNUmakefile b/ExampleCodes/heFFTe/Poisson/GNUmakefile deleted file mode 100644 index 26827791..00000000 --- a/ExampleCodes/heFFTe/Poisson/GNUmakefile +++ /dev/null @@ -1,35 +0,0 @@ -AMREX_HOME ?= ../../../../amrex -HEFFTE_HOME ?= ../../../../heffte/build - -DEBUG = FALSE -DIM = 3 -COMP = gcc -TINY_PROFILE = TRUE -USE_MPI = TRUE -USE_CUDA = FALSE -USE_HIP = FALSE - -BL_NO_FORT = TRUE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs -include $(AMREX_HOME)/Src/Base/Make.package -include Make.package - -VPATH_LOCATIONS += $(HEFFTE_HOME)/include -INCLUDE_LOCATIONS += $(HEFFTE_HOME)/include -LIBRARY_LOCATIONS += $(HEFFTE_HOME)/lib - -libraries += -lheffte - -ifeq ($(USE_CUDA),TRUE) - libraries += -lcufft -else ifeq ($(USE_HIP),TRUE) - # Use rocFFT. ROC_PATH is defined in amrex - INCLUDE_LOCATIONS += $(ROC_PATH)/rocfft/include - LIBRARY_LOCATIONS += $(ROC_PATH)/rocfft/lib - LIBRARIES += -L$(ROC_PATH)/rocfft/lib -lrocfft -else - libraries += -lfftw3_mpi -lfftw3f -lfftw3 -endif - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/heFFTe/Poisson/Make.package b/ExampleCodes/heFFTe/Poisson/Make.package deleted file mode 100644 index 6b4b865e..00000000 --- a/ExampleCodes/heFFTe/Poisson/Make.package +++ /dev/null @@ -1 +0,0 @@ -CEXE_sources += main.cpp diff --git a/ExampleCodes/heFFTe/Poisson/README b/ExampleCodes/heFFTe/Poisson/README deleted file mode 100644 index 126d04bb..00000000 --- a/ExampleCodes/heFFTe/Poisson/README +++ /dev/null @@ -1,41 +0,0 @@ -heFFTe instaallation guide: -https://icl-utk-edu.github.io/heffte/md_doxygen_installation.html - -Installation - ->> git clone git@github.com:icl-utk-edu/heffte.git -OR ->> git clone https://github.com/icl-utk-edu/heffte.git - ->> cd heffte ->> mkdir build ->> cd build - -###################### -HOST BUILD -###################### - -# NOTE: -DCMAKE_INSTALL_PREFIX can be a different location /path/to/DCMAKE_INSTALL_PREFIX - ->> cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=OFF .. - ->> make -j4 ->> make install - ->> cd /path/to/amrex-tutorials/ExampleCodes/heFFTe/Poisson ->> make -j4 HEFFTE_HOME=/path/to/DCMAKE_INSTALL_PREFIX - -###################### -NVIDIA/CUDA BUILD -###################### - -# NOTE: -DCMAKE_INSTALL_PREFIX can be a different location /path/to/DCMAKE_INSTALL_PREFIX -# NOTE: If you configure with -DHeffte_DISABLE_GPU_AWARE_MPI=OFF, you must use --gpu-bind=none in your slurm script - ->> cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_STANDARD=17 -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=. -DHeffte_ENABLE_FFTW=ON -DHeffte_ENABLE_CUDA=ON -DHeffte_DISABLE_GPU_AWARE_MPI=ON .. - ->> make -j4 ->> make install - ->> cd /path/to/amrex-tutorials/ExampleCodes/heFFTe/Poisson/ ->> make -j4 USE_CUDA=TRUE HEFFTE_HOME=/path/to/DCMAKE_INSTALL_PREFIX diff --git a/ExampleCodes/heFFTe/Poisson/inputs b/ExampleCodes/heFFTe/Poisson/inputs deleted file mode 100644 index 89815e3b..00000000 --- a/ExampleCodes/heFFTe/Poisson/inputs +++ /dev/null @@ -1,18 +0,0 @@ -n_cell_x = 64 -n_cell_y = 64 -n_cell_z = 64 - -max_grid_size_x = 32 -max_grid_size_y = 32 -max_grid_size_z = 64 - -prob_lo_x = 0. -prob_lo_y = 0. -prob_lo_z = 0. - -prob_hi_x = 1. -prob_hi_y = 1. -prob_hi_z = 1. - -laplacian_type = exact -laplacian_type = discrete diff --git a/ExampleCodes/heFFTe/Poisson/main.cpp b/ExampleCodes/heFFTe/Poisson/main.cpp deleted file mode 100644 index 7bbbb249..00000000 --- a/ExampleCodes/heFFTe/Poisson/main.cpp +++ /dev/null @@ -1,307 +0,0 @@ -#include - -#include -#include -#include -#include -#include - -using namespace amrex; - -AMREX_ENUM(LaplacianType, - Unknown, Exact, Discrete -); - -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc, argv); { - - BL_PROFILE("main"); - - // ********************************** - // DECLARE SIMULATION PARAMETERS - // ********************************** - - // number of cells on each side of the domain - int n_cell_x; - int n_cell_y; - int n_cell_z; - - // maximum grid size in each direction - int max_grid_size_x; - int max_grid_size_y; - int max_grid_size_z; - - // physical dimensions of the domain - Real prob_lo_x = 0.; - Real prob_lo_y = 0.; - Real prob_lo_z = 0.; - Real prob_hi_x = 1.; - Real prob_hi_y = 1.; - Real prob_hi_z = 1.; - - LaplacianType laplacian_type = LaplacianType::Unknown; - - // ********************************** - // READ PARAMETER VALUES FROM INPUTS FILE - // ********************************** - { - // ParmParse is way of reading inputs from the inputs file - // pp.get means we require the inputs file to have it - // pp.query means we optionally need the inputs file to have it - but you should supply a default value above - - ParmParse pp; - - pp.get("n_cell_x",n_cell_x); - pp.get("n_cell_y",n_cell_y); - pp.get("n_cell_z",n_cell_z); - - pp.get("max_grid_size_x",max_grid_size_x); - pp.get("max_grid_size_y",max_grid_size_y); - pp.get("max_grid_size_z",max_grid_size_z); - - pp.query("prob_lo_x",prob_lo_x); - pp.query("prob_lo_y",prob_lo_y); - pp.query("prob_lo_z",prob_lo_z); - - pp.query("prob_hi_x",prob_hi_x); - pp.query("prob_hi_y",prob_hi_y); - pp.query("prob_hi_z",prob_hi_z); - - pp.query_enum_case_insensitive("laplacian_type",laplacian_type); - } - - if (laplacian_type == LaplacianType::Unknown) { - amrex::Abort("Must specify exact or discrete Laplacian"); - } - - - // Determine the domain length in each direction - Real L_x = std::abs(prob_hi_x - prob_lo_x); - Real L_y = std::abs(prob_hi_y - prob_lo_y); - Real L_z = std::abs(prob_hi_z - prob_lo_z); - - // define lower and upper indices of domain - IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); - IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1)); - - // Make a single box that is the entire domain - Box domain(dom_lo, dom_hi); - - // number of points in the domain - long npts = domain.numPts(); - - // Initialize the boxarray "ba" from the single box "domain" - BoxArray ba(domain); - - // create IntVect of max_grid_size - IntVect max_grid_size(AMREX_D_DECL(max_grid_size_x,max_grid_size_y,max_grid_size_z)); - - // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction - ba.maxSize(max_grid_size); - - // How Boxes are distrubuted among MPI processes - DistributionMapping dm(ba); - - // This defines the physical box size in each direction - RealBox real_box({ AMREX_D_DECL(prob_lo_x, prob_lo_y, prob_lo_z)}, - { AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} ); - - // periodic in all direction - Array is_periodic{AMREX_D_DECL(1,1,1)}; - - // geometry object for real data - Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic); - - // extract dx from the geometry object - GpuArray dx = geom.CellSizeArray(); - - MultiFab rhs(ba,dm,1,0); - MultiFab soln(ba,dm,1,0); - - // check to make sure each MPI rank has exactly 1 box - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rhs.local_size() == 1, "Must have one Box per MPI process"); - - Real omega = M_PI/2.0; - - for (MFIter mfi(rhs); mfi.isValid(); ++mfi) { - - Array4 const& rhs_ptr = rhs.array(mfi); - - const Box& bx = mfi.fabbox(); - - amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, RandomEngine const& engine) noexcept - { - - // ********************************** - // SET VALUES FOR EACH CELL - // ********************************** - - Real x = (i+0.5) * dx[0]; - Real y = (AMREX_SPACEDIM>=2) ? (j+0.5) * dx[1] : 0.; - Real z = (AMREX_SPACEDIM==3) ? (k+0.5) * dx[2] : 0.; - - rhs_ptr(i,j,k) = std::exp(-10.*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))); - - }); - } - - // since there is 1 MPI rank per box, here each MPI rank obtains its local box and the associated boxid - Box local_box; - int local_boxid; - { - for (int i = 0; i < ba.size(); ++i) { - Box b = ba[i]; - // each MPI rank has its own local_box Box and local_boxid ID - if (ParallelDescriptor::MyProc() == dm[i]) { - local_box = b; - local_boxid = i; - } - } - } - - // now each MPI rank works on its own box - // for real->complex fft's, the fft is stored in an (nx/2+1) x ny x nz dataset - - // start by coarsening each box by 2 in the x-direction - Box c_local_box = amrex::coarsen(local_box, IntVect(AMREX_D_DECL(2,1,1))); - - // if the coarsened box's high-x index is even, we shrink the size in 1 in x - // this avoids overlap between coarsened boxes - if (c_local_box.bigEnd(0) * 2 == local_box.bigEnd(0)) { - c_local_box.setBig(0,c_local_box.bigEnd(0)-1); - } - // for any boxes that touch the hi-x domain we - // increase the size of boxes by 1 in x - // this makes the overall fft dataset have size (Nx/2+1 x Ny x Nz) - if (local_box.bigEnd(0) == geom.Domain().bigEnd(0)) { - c_local_box.growHi(0,1); - } - - // each MPI rank gets storage for its piece of the fft - BaseFab > spectral_field(c_local_box, 1, The_Device_Arena()); - - // create real->complex fft objects with the appropriate backend and data about - // the domain size and its local box size -#if (AMREX_SPACEDIM==2) -#ifdef AMREX_USE_CUDA - heffte::fft2d_r2c fft -#elif AMREX_USE_HIP - heffte::fft2d_r2c fft -#else - heffte::fft2d_r2c fft -#endif - ({{local_box.smallEnd(0),local_box.smallEnd(1),0}, - {local_box.bigEnd(0) ,local_box.bigEnd(1) ,0}}, - {{c_local_box.smallEnd(0),c_local_box.smallEnd(1),0}, - {c_local_box.bigEnd(0) ,c_local_box.bigEnd(1) ,0}}, - 0, ParallelDescriptor::Communicator()); - -#elif (AMREX_SPACEDIM==3) - -#ifdef AMREX_USE_CUDA - heffte::fft3d_r2c fft -#elif AMREX_USE_HIP - heffte::fft3d_r2c fft -#else - heffte::fft3d_r2c fft -#endif - ({{local_box.smallEnd(0),local_box.smallEnd(1),local_box.smallEnd(2)}, - {local_box.bigEnd(0) ,local_box.bigEnd(1) ,local_box.bigEnd(2)}}, - {{c_local_box.smallEnd(0),c_local_box.smallEnd(1),c_local_box.smallEnd(2)}, - {c_local_box.bigEnd(0) ,c_local_box.bigEnd(1) ,c_local_box.bigEnd(2)}}, - 0, ParallelDescriptor::Communicator()); - -#endif - - Real start_step = static_cast(ParallelDescriptor::second()); - using heffte_complex = typename heffte::fft_output::type; - heffte_complex* spectral_data = (heffte_complex*) spectral_field.dataPtr(); - - { BL_PROFILE("HEFFTE-total"); - { - BL_PROFILE("ForwardTransform"); - fft.forward(rhs[local_boxid].dataPtr(), spectral_data); - } - - // Now we take the standard FFT and scale it by 1/k^2 - Array4< GpuComplex > spectral = spectral_field.array(); - - int use_discrete = 0; - if (laplacian_type == LaplacianType::Discrete) { - use_discrete = 1; - } - - ParallelFor(c_local_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - Real a = 2.*M_PI*i / L_x; - Real b = 2.*M_PI*j / L_y; - Real c = 2.*M_PI*k / L_z; - - // the values in the upper-half of the spectral array in y and z are here interpreted as negative wavenumbers - if (j >= n_cell_y/2) b = 2.*M_PI*(n_cell_y-j) / L_y; - if (k >= n_cell_z/2) c = 2.*M_PI*(n_cell_z-k) / L_z; - - Real k2; - - if (use_discrete == 0) { -#if (AMREX_SPACEDIM == 1) - k2 = -a*a; -#elif (AMREX_SPACEDIM == 2) - k2 = -(a*a + b*b); -#elif (AMREX_SPACEDIM == 3) - k2 = -(a*a + b*b + c*c); -#endif - } else { -#if (AMREX_SPACEDIM == 1) - k2 = 2*(std::cos(a*dx[0])-1.)/(dx[0]*dx[0]); -#elif (AMREX_SPACEDIM == 2) - k2 = 2*(std::cos(a*dx[0])-1.)/(dx[0]*dx[0]) + 2*(std::cos(b*dx[1])-1.)/(dx[1]*dx[1]); -#elif (AMREX_SPACEDIM == 3) - k2 = 2*(std::cos(a*dx[0])-1.)/(dx[0]*dx[0]) + 2*(std::cos(b*dx[1])-1.)/(dx[1]*dx[1]) + 2*(std::cos(c*dx[2])-1.)/(dx[2]*dx[2]); -#endif - } - - if (k2 != 0.) { - spectral(i,j,k) /= k2; - } else { - spectral(i,j,k) *= 0.; // interpretation here is that the average value of the solution is zero - } - - }); - - - { - BL_PROFILE("BackwardTransform"); - fft.backward(spectral_data, soln[local_boxid].dataPtr()); - } - } - - // scale by 1/npts (both forward and inverse need sqrt(npts) scaling so I am doing it all here) - soln.mult(1./npts); - - Real end_step = static_cast(ParallelDescriptor::second()); - // amrex::Print() << "TIME IN SOLVE " << end_step - start_step << std::endl; - - // storage for variables to write to plotfile - MultiFab plotfile(ba, dm, 2, 0); - - // copy rhs and soln into plotfile - MultiFab::Copy(plotfile, rhs , 0, 0, 1, 0); - MultiFab::Copy(plotfile, soln, 0, 1, 1, 0); - - // time and step are dummy variables required to WriteSingleLevelPlotfile - Real time = 0.; - int step = 0; - - // arguments - // 1: name of plotfile - // 2: MultiFab containing data to plot - // 3: variables names - // 4: geometry object - // 5: "time" of plotfile; not relevant in this example - // 6: "time step" of plotfile; not relevant in this example - WriteSingleLevelPlotfile("plt", plotfile, {"rhs", "soln"}, geom, time, step); - - } amrex::Finalize(); -} diff --git a/ExampleCodes/heFFTe/Poisson2/GNUmakefile b/ExampleCodes/heFFTe/Poisson2/GNUmakefile deleted file mode 100644 index 3b4e52fa..00000000 --- a/ExampleCodes/heFFTe/Poisson2/GNUmakefile +++ /dev/null @@ -1,35 +0,0 @@ -AMREX_HOME ?= ../../../../amrex -HEFFTE_HOME ?= ../../../../heffte/build - -DEBUG = FALSE -DIM = 3 -COMP = gcc -TINY_PROFILE = FALSE -USE_MPI = TRUE -USE_CUDA = FALSE -USE_HIP = FALSE - -BL_NO_FORT = TRUE - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs -include $(AMREX_HOME)/Src/Base/Make.package -include Make.package - -VPATH_LOCATIONS += $(HEFFTE_HOME)/include -INCLUDE_LOCATIONS += $(HEFFTE_HOME)/include -LIBRARY_LOCATIONS += $(HEFFTE_HOME)/lib - -libraries += -lheffte - -ifeq ($(USE_CUDA),TRUE) - libraries += -lcufft -else ifeq ($(USE_HIP),TRUE) - # Use rocFFT. ROC_PATH is defined in amrex - INCLUDE_LOCATIONS += $(ROC_PATH)/rocfft/include - LIBRARY_LOCATIONS += $(ROC_PATH)/rocfft/lib - LIBRARIES += -L$(ROC_PATH)/rocfft/lib -lrocfft -else - libraries += -lfftw3_mpi -lfftw3f -lfftw3 -endif - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/heFFTe/Poisson2/Make.package b/ExampleCodes/heFFTe/Poisson2/Make.package deleted file mode 100644 index 6b4b865e..00000000 --- a/ExampleCodes/heFFTe/Poisson2/Make.package +++ /dev/null @@ -1 +0,0 @@ -CEXE_sources += main.cpp diff --git a/ExampleCodes/heFFTe/Poisson2/inputs b/ExampleCodes/heFFTe/Poisson2/inputs deleted file mode 100644 index 9b888704..00000000 --- a/ExampleCodes/heFFTe/Poisson2/inputs +++ /dev/null @@ -1,11 +0,0 @@ -n_cell_x = 64 -n_cell_y = 64 -n_cell_z = 64 - -prob_lo_x = 0. -prob_lo_y = 0. -prob_lo_z = 0. - -prob_hi_x = 1. -prob_hi_y = 1. -prob_hi_z = 1. diff --git a/ExampleCodes/heFFTe/Poisson2/main.cpp b/ExampleCodes/heFFTe/Poisson2/main.cpp deleted file mode 100644 index 34727023..00000000 --- a/ExampleCodes/heFFTe/Poisson2/main.cpp +++ /dev/null @@ -1,314 +0,0 @@ -#include - -#include -#include -#include -#include -#include - -using namespace amrex; - -static_assert(AMREX_SPACEDIM == 3); - -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc, argv); { - - BL_PROFILE("main"); - - // ********************************** - // DECLARE SIMULATION PARAMETERS - // ********************************** - - // number of cells on each side of the domain - int n_cell_x; - int n_cell_y; - int n_cell_z; - - // physical dimensions of the domain - Real prob_lo_x = 0.; - Real prob_lo_y = 0.; - Real prob_lo_z = 0.; - Real prob_hi_x = 1.; - Real prob_hi_y = 1.; - Real prob_hi_z = 1.; - - // ********************************** - // READ PARAMETER VALUES FROM INPUTS FILE - // ********************************** - { - // ParmParse is way of reading inputs from the inputs file - // pp.get means we require the inputs file to have it - // pp.query means we optionally need the inputs file to have it - but you should supply a default value above - - ParmParse pp; - - pp.get("n_cell_x",n_cell_x); - pp.get("n_cell_y",n_cell_y); - pp.get("n_cell_z",n_cell_z); - - pp.query("prob_lo_x",prob_lo_x); - pp.query("prob_lo_y",prob_lo_y); - pp.query("prob_lo_z",prob_lo_z); - - pp.query("prob_hi_x",prob_hi_x); - pp.query("prob_hi_y",prob_hi_y); - pp.query("prob_hi_z",prob_hi_z); - } - - // Determine the domain length in each direction - Real L_x = std::abs(prob_hi_x - prob_lo_x); - Real L_y = std::abs(prob_hi_y - prob_lo_y); - - // define lower and upper indices of domain - IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); - IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1)); - - // Make a single box that is the entire domain - Box domain(dom_lo, dom_hi); - - // Initialize the boxarray "ba" from the single box "domain" There are - // exactly nprocs boxes. The domain decomposition is done in the x- and - // y-directions, but not the z-direction. - BoxArray ba = amrex::decompose(domain, ParallelDescriptor::NProcs(), - {AMREX_D_DECL(true,true,false)}); - - // How Boxes are distrubuted among MPI processes - DistributionMapping dm(ba); - - // This defines the physical box size in each direction - RealBox real_box({ AMREX_D_DECL(prob_lo_x, prob_lo_y, prob_lo_z)}, - { AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} ); - - // periodic in all direction - Array is_periodic{AMREX_D_DECL(1,1,1)}; - - // geometry object for real data - Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic); - - // extract dx from the geometry object - GpuArray dx = geom.CellSizeArray(); - - MultiFab rhs(ba,dm,1,0); - MultiFab soln(ba,dm,1,0); - - // check to make sure each MPI rank has exactly 1 box - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rhs.local_size() == 1, "Must have one Box per MPI process"); - - for (MFIter mfi(rhs); mfi.isValid(); ++mfi) { - - Array4 const& rhs_ptr = rhs.array(mfi); - - const Box& bx = mfi.fabbox(); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - // ********************************** - // SET VALUES FOR EACH CELL - // ********************************** - - Real x = (i+0.5) * dx[0]; - Real y = (AMREX_SPACEDIM>=2) ? (j+0.5) * dx[1] : 0.; - Real z = (AMREX_SPACEDIM==3) ? (k+0.5) * dx[2] : 0.; - - rhs_ptr(i,j,k) = std::exp(-10.*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))); - - }); - } - - // Shift rhs so that its sum is zero. - auto rhosum = rhs.sum(0); - rhs.plus(-rhosum/geom.Domain().d_numPts(), 0, 1); - - // since there is 1 MPI rank per box, here each MPI rank obtains its local box and the associated boxid - Box local_box; - int local_boxid; - { - for (int i = 0; i < ba.size(); ++i) { - Box b = ba[i]; - // each MPI rank has its own local_box Box and local_boxid ID - if (ParallelDescriptor::MyProc() == dm[i]) { - local_box = b; - local_boxid = i; - } - } - } - - // now each MPI rank works on its own box - // for real->complex fft's, the fft is stored in an (nx/2+1) x ny x nz dataset - - // start by coarsening each box by 2 in the x-direction - Box c_local_box = amrex::coarsen(local_box, IntVect(AMREX_D_DECL(2,1,1))); - - // if the coarsened box's high-x index is even, we shrink the size in 1 in x - // this avoids overlap between coarsened boxes - if (c_local_box.bigEnd(0) * 2 == local_box.bigEnd(0)) { - c_local_box.setBig(0,c_local_box.bigEnd(0)-1); - } - // for any boxes that touch the hi-x domain we - // increase the size of boxes by 1 in x - // this makes the overall fft dataset have size (Nx/2+1 x Ny x Nz) - if (local_box.bigEnd(0) == geom.Domain().bigEnd(0)) { - c_local_box.growHi(0,1); - } - - // each MPI rank gets storage for its piece of the fft - BaseFab > spectral_field(c_local_box, 1, The_Device_Arena()); - - // create real->complex fft objects with the appropriate backend and data about - // the domain size and its local box size - using fft_r2c_t = -#ifdef AMREX_USE_CUDA - heffte::fft2d_r2c; -#elif AMREX_USE_HIP - heffte::fft2d_r2c; -#else - heffte::fft2d_r2c; -#endif - - auto lo = amrex::lbound(local_box); - auto hi = amrex::ubound(local_box); - auto len = amrex::length(local_box); - auto clo = amrex::lbound(c_local_box); - auto chi = amrex::ubound(c_local_box); - auto clen = amrex::length(c_local_box); - - auto fft = std::make_unique - (heffte::box3d({ lo.x, lo.y, 0}, - { hi.x, hi.y, 0}), - heffte::box3d({clo.x, clo.y, 0}, - {chi.x, chi.y, 0}), - 0, ParallelDescriptor::Communicator()); - - Real start_step = static_cast(ParallelDescriptor::second()); - using heffte_complex = typename heffte::fft_output::type; - heffte_complex* spectral_data = (heffte_complex*) spectral_field.dataPtr(); - - int batch_size = n_cell_z; - Gpu::DeviceVector workspace(fft->size_workspace()*batch_size); - - { BL_PROFILE("HEFFTE-total"); - { - BL_PROFILE("ForwardTransform"); - fft->forward(batch_size, rhs[local_boxid].dataPtr(), spectral_data, - workspace.data()); - } - - // Now we take the standard FFT and scale it by 1/k^2 - Array4< GpuComplex > spectral = spectral_field.array(); - - FArrayBox tridiag_workspace(c_local_box,4); - auto const& ald = tridiag_workspace.array(0); - auto const& bd = tridiag_workspace.array(1); - auto const& cud = tridiag_workspace.array(2); - auto const& scratch = tridiag_workspace.array(3); - - Gpu::DeviceVector delzv(n_cell_z, dx[2]); - auto const* delz = delzv.data(); - - auto xybox = amrex::makeSlab(c_local_box, 2, 0); - ParallelFor(xybox, [=] AMREX_GPU_DEVICE(int i, int j, int) - { - Real a = 2.*M_PI*i / L_x; - Real b = 2.*M_PI*j / L_y; - - // the values in the upper-half of the spectral array in y and z are here interpreted as negative wavenumbers - if (j >= n_cell_y/2) b = 2.*M_PI*(n_cell_y-j) / L_y; - - Real k2 = 2*(std::cos(a*dx[0])-1.)/(dx[0]*dx[0]) + 2*(std::cos(b*dx[1])-1.)/(dx[1]*dx[1]); - - // Tridiagonal solve with homogeneous Neumann - for( int k=0; k= 0; k--) { - spectral(i,j,k) -= scratch(i,j,k) * spectral(i,j,k + 1); - } - }); - Gpu::streamSynchronize(); - - { - BL_PROFILE("BackwardTransform"); - fft->backward(batch_size, spectral_data, soln[local_boxid].dataPtr(), - heffte::scale::full); - } - } - - Real end_step = static_cast(ParallelDescriptor::second()); - // amrex::Print() << "TIME IN SOLVE " << end_step - start_step << std::endl; - - // storage for variables to write to plotfile - MultiFab plotfile(ba, dm, 2, 0); - - // copy rhs and soln into plotfile - MultiFab::Copy(plotfile, rhs , 0, 0, 1, 0); - MultiFab::Copy(plotfile, soln, 0, 1, 1, 0); - - // time and step are dummy variables required to WriteSingleLevelPlotfile - Real time = 0.; - int step = 0; - - // arguments - // 1: name of plotfile - // 2: MultiFab containing data to plot - // 3: variables names - // 4: geometry object - // 5: "time" of plotfile; not relevant in this example - // 6: "time step" of plotfile; not relevant in this example - WriteSingleLevelPlotfile("plt", plotfile, {"rhs", "soln"}, geom, time, step); - - { - MultiFab phi(soln.boxArray(), soln.DistributionMap(), 1, 1); - MultiFab res(soln.boxArray(), soln.DistributionMap(), 1, 0); - MultiFab::Copy(phi, soln, 0, 0, 1, 0); - phi.FillBoundary(geom.periodicity()); - auto const& res_ma = res.arrays(); - auto const& phi_ma = phi.const_arrays(); - auto const& rhs_ma = rhs.const_arrays(); - ParallelFor(res, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) - { - auto const& phia = phi_ma[b]; - auto lap = (phia(i-1,j,k)-2.*phia(i,j,k)+phia(i+1,j,k)) / (dx[0]*dx[0]) - + (phia(i,j-1,k)-2.*phia(i,j,k)+phia(i,j+1,k)) / (dx[1]*dx[1]); - if (k == 0) { - lap += (-phia(i,j,k)+phia(i,j,k+1)) / (dx[2]*dx[2]); - } else if (k == n_cell_z-1) { - lap += (phia(i,j,k-1)-phia(i,j,k)) / (dx[2]*dx[2]); - } else { - lap += (phia(i,j,k-1)-2.*phia(i,j,k)+phia(i,j,k+1)) / (dx[2]*dx[2]); - } - res_ma[b](i,j,k) = rhs_ma[b](i,j,k) - lap; - }); - amrex::Print() << " rhs.min & max: " << rhs.min(0) << " " << rhs.max(0) << "\n" - << " res.min & max: " << res.min(0) << " " << res.max(0) << "\n"; - } - - } amrex::Finalize(); -}