Skip to content

Add scaling for non-cubic rectangular geometries, remove mean velocit… #170

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
drummerdoc wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add scaling for non-cubic rectangular geometries, remove mean velocit… #170

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 27 additions & 14 deletions ExampleCodes/FFT/EnergySpectrum/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,25 +1,39 @@
#include <AMReX.H>
#include <AMReX_FFT.H>
#include <AMReX_PlotFileUtil.H>
#include <AMReX_ParmParse.H>

using namespace amrex;

int main (int argc, char* argv[])
{
amrex::Initialize(argc, argv);
{
std::string infile, outfile;
Vector<std::string> vel_labels = {AMREX_D_DECL("velx","vely","velz")};
ParmParse pp;
pp.get("infile",infile);
Copy link
Copy Markdown
Member

@WeiqunZhang WeiqunZhang Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe we should use query to maintain backward compatibility.

pp.get("outfile",outfile);
pp.queryarr("vel_labels",vel_labels,0,3);
Geometry geom;
MultiFab vx, vy, vz;
Vector<Real> scaleL(3);
{
PlotFileData plot("plot");
geom.define(plot.probDomain(0),
RealBox(plot.probLo(), plot.probHi()),
plot.coordSys(), {1,1,1});
vx = plot.get(0, "velx");
vy = plot.get(0, "vely");
vz = plot.get(0, "velz");
}

PlotFileData plot(infile);
geom.define(plot.probDomain(0),
RealBox(plot.probLo(), plot.probHi()),
plot.coordSys(), {1,1,1});
// Read data, remove mean field
vx = plot.get(0, vel_labels[0]); vx.mult(1./vx.sum(0,0));
vy = plot.get(0, vel_labels[1]); vy.mult(1./vy.sum(0,0));
vz = plot.get(0, vel_labels[2]); vz.mult(1./vz.sum(0,0));
Real L = std::min({geom.ProbLength(0),
geom.ProbLength(1),
geom.ProbLength(2)});
scaleL = {AMREX_D_DECL(1./std::round(geom.ProbLength(0)/L),
1./std::round(geom.ProbLength(1)/L),
1./std::round(geom.ProbLength(2)/L))};
}
cMultiFab cx, cy, cz;
{
// Note that the complex Hermitian output array Y has (nx/2+1,ny,nz) elements.
Expand All @@ -34,8 +48,7 @@ int main (int argc, char* argv[])
fft.forward(vz, cz);
}

// For simplicity, we assume the domain is a cube, and we are not
// going to worry about scaling.
// For simplicity, we are not going to worry about scaling.

int nx = geom.Domain().length(0);
int ny = geom.Domain().length(1);
Expand All @@ -51,7 +64,7 @@ int main (int argc, char* argv[])
int ki = i;
int kj = (j <= ny/2) ? j : ny-j;
int kk = (k <= nz/2) ? k : nz-k;
Real d = std::sqrt(Real(ki*ki+kj*kj+kk*kk));
Real d = std::sqrt(Real(ki*ki*scaleL[0]+kj*kj*scaleL[1]+kk*kk*scaleL[2]));
int di = int(std::round(d));
if (di < nk) {
Real value = amrex::norm(cxa[b](i,j,k))
Expand Down Expand Up @@ -79,8 +92,8 @@ int main (int argc, char* argv[])
ParallelDescriptor::ReduceRealSum(pke_h, nk);

if (ParallelDescriptor::IOProcessor()) {
std::ofstream ofs("spectrum.txt");
for (int i = 0; i < nk; ++i) {
std::ofstream ofs(outfile.c_str());
for (int i = 0; i < nk; ++i) {
ofs << i << " " << pke_h[i] << "\n";
}
}
Expand Down
Loading