Skip to content
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

Jump to bottom

add option to solve discrete rather than exact laplacian #137

Merged
merged 2 commits into from
Oct 2, 2024
Merged

add option to solve discrete rather than exact laplacian #137

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b9e1740
add option to solve discrete rather than exact laplacian
asalmgren Sep 28, 2024
6fb49b5
clean up enum
asalmgren Oct 2, 2024
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
3 changes: 3 additions & 0 deletions ExampleCodes/heFFTe/Poisson/inputs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,3 +13,6 @@ prob_lo_z = 0.
prob_hi_x = 1.
prob_hi_y = 1.
prob_hi_z = 1.

laplacian_type = exact
laplacian_type = discrete
56 changes: 50 additions & 6 deletions ExampleCodes/heFFTe/Poisson/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,10 @@

using namespace amrex;

enum struct LaplacianType {
Unknown, Exact, Discrete
};

Copy link
Member

@WeiqunZhang WeiqunZhang Oct 2, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Suggested change
enum struct LaplacianType {
Unknown, Exact, Discrete
};
AMREX_ENUM(LaplacianType,
Unknown, Exact, Discrete
);

Then we can query LaplacianType directly without using string.

pp.query_enum_case_insensitive("laplacian_type", laplacian_type)

int main (int argc, char* argv[])
{
amrex::Initialize(argc, argv); {
Expand Down Expand Up @@ -36,6 +40,9 @@ int main (int argc, char* argv[])
Real prob_hi_y = 1.;
Real prob_hi_z = 1.;

std::string laplacian_type_string;
LaplacianType laplacian_type = LaplacianType::Unknown;

// **********************************
// READ PARAMETER VALUES FROM INPUTS FILE
// **********************************
Expand All @@ -61,8 +68,26 @@ int main (int argc, char* argv[])
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("laplacian_type",laplacian_type_string);
}

if ( (laplacian_type_string == "Exact") ||
(laplacian_type_string == "exact") ) {
laplacian_type = LaplacianType::Exact;
}
else
if ( (laplacian_type_string == "Discrete") ||
(laplacian_type_string == "discrete") ) {
laplacian_type = LaplacianType::Discrete;
}


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);
Expand All @@ -77,7 +102,6 @@ int main (int argc, char* argv[])

// 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);
Expand Down Expand Up @@ -172,7 +196,6 @@ int main (int argc, char* argv[])
// 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<heffte::backend::cufft> fft
#elif AMREX_USE_HIP
Expand Down Expand Up @@ -203,6 +226,7 @@ int main (int argc, char* argv[])

#endif

Real start_step = static_cast<Real>(ParallelDescriptor::second());
using heffte_complex = typename heffte::fft_output<Real>::type;
heffte_complex* spectral_data = (heffte_complex*) spectral_field.dataPtr();

Expand All @@ -215,6 +239,11 @@ int main (int argc, char* argv[])
// Now we take the standard FFT and scale it by 1/k^2
Array4< GpuComplex<Real> > 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;
Expand All @@ -225,13 +254,25 @@ int main (int argc, char* argv[])
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)
Real k2 = -a*a;
k2 = 2*(std::cos(a*dx[0])-1.)/(dx[0]*dx[0]);
#elif (AMREX_SPACEDIM == 2)
Real k2 = -(a*a + b*b);
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)
Real k2 = -(a*a + b*b + c*c);
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;
Expand All @@ -248,9 +289,12 @@ int main (int argc, char* argv[])
}
}

// scale by 1/npts (both forward and inverse need 1/sqrtnpts scaling so I am doing it all here)
// 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<Real>(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);

Expand Down
Loading