Skip to content

Commit c2175c2

Browse files
ajnonakaajnonaka
authored
Tutorials for Basic and Poisson using amrex::FFT (#140)
* simple Poisson example using amrex::FFT::Poisson * simple example of forward and inverse FFT and visualization of FFT data * remove discrete lap option * FFT tutorial documentation
1 parent 2ef516b commit c2175c2

File tree

10 files changed

+459
-0
lines changed

10 files changed

+459
-0
lines changed

‎Docs/source/FFT_Tutorial.rst‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
.. role:: cpp(code)
2+
:language: c++
3+
4+
.. role:: fortran(code)
5+
:language: fortran
6+
7+
.. _tutorials_fft:
8+
9+
FFT
10+
==========================
11+
12+
These tutorials demonstrate how to use the amrex::FFT classes to solve for and manipulate Fourier transform data.
13+
For more information on the amrex::FFT class refer to the AMReX documentation
14+
`here <https://amrex-codes.github.io/amrex/docs_html/FFT_Chapter.html>`_.
15+
16+
There are two FFT tutorials, ``Basic`` and ``Poisson``.
17+
18+
.. toctree::
19+
:maxdepth: 1
20+
21+
.. _section:fft_tutorial:fft_basic:
22+
23+
Basic
24+
--------------------------
25+
26+
The tutorial found in ``amrex-tutorials/ExampleCodes/FFT/Basic`` demonstrates how
27+
to take a forward FFT, manipulate or access the spectral data (in this example by copying the spectral
28+
data into a :cpp:`MultiFab` for plotfile visualziation) and then taking the inverse FFT.
29+
30+
.. _section:fft_tutorial:fft_pois:
31+
32+
Poisson
33+
--------------------------
34+
35+
This tutorial: ``amrex-tutorials/ExampleCodes/FFT/Poisson``
36+
solves a Poisson equation with periodic boundary conditions. It relies on the :cpp:`AMReX_FFT_Poisson.H`
37+
routines to solve the equation by using the ``forwardThenBackward`` function which takes the forward and inverse
38+
transform of the data with spectral data scaling in between as required by the Poisson equation.

‎Docs/source/index.rst‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,7 @@ sorted by the following categories:
4545
Conduit Mesh Blueprint for use with the ALPINE Ascent in situ visualization
4646
and analysis tool.
4747
- :ref:`EB<tutorials_eb>` -- Examples of embedded boundaries.
48+
- :ref:`FFT<tutorials_fft>` -- Demonstration on how to use the amrex::FFT class for forward/inverse FFT and data manipulation
4849
- :ref:`FFTW<tutorials_fftw>` -- FFTW and cuFFT single-rank tutorials.
4950
- :ref:`ForkJoin<tutorials_forkjoin>` -- Parallel execution and subgrouping of MPI ranks.
5051
- :ref:`GPU<tutorials_gpu>` -- Offload work to the GPUs using AMReX tools.
@@ -70,6 +71,7 @@ sorted by the following categories:
7071
Basic_Tutorial
7172
Blueprint_Tutorial
7273
EB_Tutorial
74+
FFT_Tutorial
7375
FFTW_Tutorial
7476
ForkJoin_Tutorial
7577
GPU_Tutorial
@@ -94,6 +96,8 @@ sorted by the following categories:
9496

9597
.. _`EB`: EB_Tutorial.html
9698

99+
.. _`FFT`: FFT_Tutorial.html
100+
97101
.. _`FFTW`: FFTW_Tutorial.html
98102

99103
.. _`ForkJoin`: ForkJoin_Tutorial.html

‎ExampleCodes/FFT/Basic/GNUmakefile‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
AMREX_HOME ?= ../../../../amrex
2+
3+
DEBUG = FALSE
4+
DIM = 3
5+
COMP = gcc
6+
TINY_PROFILE = TRUE
7+
USE_MPI = TRUE
8+
USE_CUDA = FALSE
9+
USE_HIP = FALSE
10+
USE_FFT = TRUE
11+
12+
BL_NO_FORT = TRUE
13+
14+
include $(AMREX_HOME)/Tools/GNUMake/Make.defs
15+
include $(AMREX_HOME)/Src/Base/Make.package
16+
include $(AMREX_HOME)/Src/FFT/Make.package
17+
include Make.package
18+
19+
20+
include $(AMREX_HOME)/Tools/GNUMake/Make.rules

‎ExampleCodes/FFT/Basic/Make.package‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
CEXE_sources += main.cpp

‎ExampleCodes/FFT/Basic/inputs‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
n_cell_x = 64
2+
n_cell_y = 64
3+
n_cell_z = 64
4+
5+
max_grid_size_x = 32
6+
max_grid_size_y = 32
7+
max_grid_size_z = 64
8+
9+
prob_lo_x = 0.
10+
prob_lo_y = 0.
11+
prob_lo_z = 0.
12+
13+
prob_hi_x = 1.
14+
prob_hi_y = 1.
15+
prob_hi_z = 1.

‎ExampleCodes/FFT/Basic/main.cpp‎

Lines changed: 190 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,190 @@
1+
#include <AMReX.H>
2+
#include <AMReX_MultiFab.H>
3+
#include <AMReX_ParmParse.H>
4+
#include <AMReX_PlotFileUtil.H>
5+
6+
#include <AMReX_FFT.H>
7+
8+
using namespace amrex;
9+
10+
int main (int argc, char* argv[])
11+
{
12+
amrex::Initialize(argc, argv); {
13+
14+
BL_PROFILE("main");
15+
16+
// **********************************
17+
// DECLARE SIMULATION PARAMETERS
18+
// **********************************
19+
20+
// number of cells on each side of the domain
21+
int n_cell_x;
22+
int n_cell_y;
23+
int n_cell_z;
24+
25+
// maximum grid size in each direction
26+
int max_grid_size_x;
27+
int max_grid_size_y;
28+
int max_grid_size_z;
29+
30+
// physical dimensions of the domain
31+
Real prob_lo_x = 0.;
32+
Real prob_lo_y = 0.;
33+
Real prob_lo_z = 0.;
34+
Real prob_hi_x = 1.;
35+
Real prob_hi_y = 1.;
36+
Real prob_hi_z = 1.;
37+
38+
// **********************************
39+
// READ PARAMETER VALUES FROM INPUTS FILE
40+
// **********************************
41+
{
42+
// ParmParse is way of reading inputs from the inputs file
43+
// pp.get means we require the inputs file to have it
44+
// pp.query means we optionally need the inputs file to have it - but you should supply a default value above
45+
46+
ParmParse pp;
47+
48+
pp.get("n_cell_x",n_cell_x);
49+
pp.get("n_cell_y",n_cell_y);
50+
pp.get("n_cell_z",n_cell_z);
51+
52+
pp.get("max_grid_size_x",max_grid_size_x);
53+
pp.get("max_grid_size_y",max_grid_size_y);
54+
pp.get("max_grid_size_z",max_grid_size_z);
55+
56+
pp.query("prob_lo_x",prob_lo_x);
57+
pp.query("prob_lo_y",prob_lo_y);
58+
pp.query("prob_lo_z",prob_lo_z);
59+
60+
pp.query("prob_hi_x",prob_hi_x);
61+
pp.query("prob_hi_y",prob_hi_y);
62+
pp.query("prob_hi_z",prob_hi_z);
63+
}
64+
65+
// Determine the domain length in each direction
66+
Real L_x = std::abs(prob_hi_x - prob_lo_x);
67+
Real L_y = std::abs(prob_hi_y - prob_lo_y);
68+
Real L_z = std::abs(prob_hi_z - prob_lo_z);
69+
70+
// define lower and upper indices of domain
71+
IntVect dom_lo(AMREX_D_DECL( 0, 0, 0));
72+
IntVect dom_hi(AMREX_D_DECL(n_cell_x-1, n_cell_y-1, n_cell_z-1));
73+
74+
// Make a single box that is the entire domain
75+
Box domain(dom_lo, dom_hi);
76+
77+
// number of points in the domain
78+
long npts = domain.numPts();
79+
80+
// Initialize the boxarray "ba" from the single box "domain"
81+
BoxArray ba(domain);
82+
83+
// create IntVect of max_grid_size
84+
IntVect max_grid_size(AMREX_D_DECL(max_grid_size_x,max_grid_size_y,max_grid_size_z));
85+
86+
// Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction
87+
ba.maxSize(max_grid_size);
88+
89+
// How Boxes are distrubuted among MPI processes
90+
DistributionMapping dm(ba);
91+
92+
// This defines the physical box size in each direction
93+
RealBox real_box({ AMREX_D_DECL(prob_lo_x, prob_lo_y, prob_lo_z)},
94+
{ AMREX_D_DECL(prob_hi_x, prob_hi_y, prob_hi_z)} );
95+
96+
// periodic in all direction
97+
Array<int,AMREX_SPACEDIM> is_periodic{AMREX_D_DECL(1,1,1)};
98+
99+
// geometry object for real data
100+
Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic);
101+
102+
// extract dx from the geometry object
103+
GpuArray<Real,AMREX_SPACEDIM> dx = geom.CellSizeArray();
104+
105+
amrex::FFT::R2C my_fft(domain);
106+
107+
// storage for input phi and phi after forward-inverse transformation
108+
MultiFab phi(ba,dm,1,0);
109+
MultiFab phi_after(ba,dm,1,0);
110+
111+
// initialize phi
112+
for (MFIter mfi(phi); mfi.isValid(); ++mfi) {
113+
114+
Array4<Real> const& phi_ptr = phi.array(mfi);
115+
116+
const Box& bx = mfi.fabbox();
117+
118+
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
119+
{
120+
121+
// **********************************
122+
// SET VALUES FOR EACH CELL
123+
// **********************************
124+
125+
Real x = (i+0.5) * dx[0];
126+
Real y = (AMREX_SPACEDIM>=2) ? (j+0.5) * dx[1] : 0.;
127+
Real z = (AMREX_SPACEDIM==3) ? (k+0.5) * dx[2] : 0.;
128+
129+
phi_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)));
130+
131+
});
132+
}
133+
134+
// create storage for the FFT
135+
Box cdomain = geom.Domain();
136+
cdomain.setBig(0,cdomain.length(0)/2);
137+
Geometry cgeom(cdomain, real_box, CoordSys::cartesian, is_periodic);
138+
auto cba = amrex::decompose(cdomain, ParallelContext::NProcsSub(),
139+
{AMREX_D_DECL(true,true,false)});
140+
DistributionMapping cdm(cba);
141+
FabArray<BaseFab<GpuComplex<amrex::Real> > > phi_fft(cba, cdm, 1, 0);
142+
143+
// we will copy the real and imaginary parts of the FFT to this MultiFab
144+
MultiFab phi_fft_realimag(cba,cdm,2,0);
145+
146+
// compute the FFT
147+
my_fft.forward(phi,phi_fft);
148+
149+
// copy FFT into a MultiFab for plotfile visualization
150+
for (MFIter mfi(phi_fft); mfi.isValid(); ++mfi) {
151+
152+
Array4<GpuComplex<Real>> const& phi_fft_ptr = phi_fft.array(mfi);
153+
Array4<Real> phi_fft_realimag_ptr = phi_fft_realimag.array(mfi);
154+
155+
const Box& bx = mfi.fabbox();
156+
157+
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
158+
{
159+
phi_fft_realimag_ptr(i,j,k,0) = phi_fft_ptr(i,j,k).real();
160+
phi_fft_realimag_ptr(i,j,k,1) = phi_fft_ptr(i,j,k).imag();
161+
162+
});
163+
}
164+
165+
// compute the inverse FFT and store result in phi_after
166+
my_fft.backward(phi_after);
167+
168+
// scale phi_after by 1/n_cells so it matches the original phi
169+
long n_cells = n_cell_x;
170+
if (AMREX_SPACEDIM >= 2) n_cells *= n_cell_y;
171+
if (AMREX_SPACEDIM >= 3) n_cells *= n_cell_z;
172+
phi_after.mult(1./n_cells);
173+
174+
// time and step are dummy variables required to WriteSingleLevelPlotfile
175+
Real time = 0.;
176+
int step = 0;
177+
178+
// arguments
179+
// 1: name of plotfile
180+
// 2: MultiFab containing data to plot
181+
// 3: variables names
182+
// 4: geometry object
183+
// 5: "time" of plotfile; not relevant in this example
184+
// 6: "time step" of plotfile; not relevant in this example
185+
WriteSingleLevelPlotfile("plt", phi, {"phi"}, geom, time, step);
186+
WriteSingleLevelPlotfile("plt_after", phi_after, {"phi_after"}, geom, time, step);
187+
WriteSingleLevelPlotfile("plt_fft", phi_fft_realimag, {"phi_fft_real", "phi_fft_imag"}, cgeom, time, step);
188+
189+
} amrex::Finalize();
190+
}

‎ExampleCodes/FFT/Poisson/GNUmakefile‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
AMREX_HOME ?= ../../../../amrex
2+
3+
DEBUG = FALSE
4+
DIM = 3
5+
COMP = gcc
6+
TINY_PROFILE = TRUE
7+
USE_MPI = TRUE
8+
USE_CUDA = FALSE
9+
USE_HIP = FALSE
10+
USE_FFT = TRUE
11+
12+
BL_NO_FORT = TRUE
13+
14+
include $(AMREX_HOME)/Tools/GNUMake/Make.defs
15+
include $(AMREX_HOME)/Src/Base/Make.package
16+
include $(AMREX_HOME)/Src/FFT/Make.package
17+
include Make.package
18+
19+
20+
include $(AMREX_HOME)/Tools/GNUMake/Make.rules

‎ExampleCodes/FFT/Poisson/Make.package‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
CEXE_sources += main.cpp

‎ExampleCodes/FFT/Poisson/inputs‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
n_cell_x = 64
2+
n_cell_y = 64
3+
n_cell_z = 64
4+
5+
max_grid_size_x = 32
6+
max_grid_size_y = 32
7+
max_grid_size_z = 64
8+
9+
prob_lo_x = 0.
10+
prob_lo_y = 0.
11+
prob_lo_z = 0.
12+
13+
prob_hi_x = 1.
14+
prob_hi_y = 1.
15+
prob_hi_z = 1.

0 commit comments

Comments
 (0)