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

Minor improvements to smoother #267

Merged
merged 4 commits into from
Mar 31, 2021
Merged

Minor improvements to smoother #267

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5658a88
Remove un-needed pfloat -> dfloat operation at the end of MGLevel::sm…
Mar 30, 2021
7ff7562
Remove Richardson smoother, prefer enum class to enum
Mar 30, 2021
f657c2b
Fuse conversion of two vectors from dfloat type to pfloat type
Mar 30, 2021
a12efc2
Remove unused method definitions in MGLevel
Mar 30, 2021
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
12 changes: 12 additions & 0 deletions okl/elliptic/ellipticLinAlg.okl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,18 @@ SOFTWARE.
}
}

@kernel void fusedCopyDfloatToPfloat(const dlong N,
@restrict const dfloat * x_dfloat,
@restrict const dfloat * y_dfloat,
@restrict pfloat * x_pfloat,
@restrict pfloat * y_pfloat){
for(dlong n=0;n<N;++n;@tile(p_blockSize,@outer,@inner)){
if(n<N){
x_pfloat[n] = x_dfloat[n];
y_pfloat[n] = y_dfloat[n];
}
}
}
@kernel void copyDfloatToPfloat(const dlong N,
@restrict pfloat * y,
@restrict const dfloat * x){
Expand Down
1 change: 1 addition & 0 deletions src/elliptic/elliptic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,7 @@ struct elliptic_t
occa::kernel scaledAddPfloatKernel;
occa::kernel dotMultiplyPfloatKernel;
occa::kernel copyDfloatToPfloatKernel;
occa::kernel fusedCopyDfloatToPfloatKernel;
occa::kernel copyPfloatToDPfloatKernel;
occa::kernel updateSmoothedSolutionVecKernel;
occa::kernel updateChebyshevSolutionVecKernel;
Expand Down
23 changes: 13 additions & 10 deletions src/elliptic/ellipticMultiGrid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,16 @@
#ifndef ELLIPTIC_MGLEVEL_HPP
#define ELLIPTIC_MGLEVEL_HPP

typedef enum {RICHARDSON = 1,
CHEBYSHEV = 2,
SCHWARZ = 3} SmoothType;
typedef enum {JACOBI = 1,
SCHWARZ_SMOOTH = 3} SmootherType;
enum class SmootherType
{
CHEBYSHEV,
SCHWARZ,
};
enum class SecondarySmootherType
{
JACOBI,
SCHWARZ,
};

class MGLevel : public parAlmond::multigridLevel
{
Expand All @@ -50,9 +55,9 @@ class MGLevel : public parAlmond::multigridLevel
occa::memory o_invDegree;

//smoothing params
SmoothType stype;
SmootherType smtypeUp;
SmootherType smtypeDown;
SmootherType stype;
SecondarySmootherType smtypeUp;
SecondarySmootherType smtypeDown;

dfloat lambda1, lambda0;
int ChebyshevIterations;
Expand Down Expand Up @@ -128,7 +133,6 @@ class MGLevel : public parAlmond::multigridLevel

void smoother(occa::memory o_x, occa::memory o_Sx, bool xIsZero);

void smoothRichardson(occa::memory &o_r, occa::memory &o_x, bool xIsZero);
void smoothChebyshev (occa::memory &o_r, occa::memory &o_x, bool xIsZero);
void smoothSchwarz (occa::memory &o_r, occa::memory &o_x, bool xIsZero);

Expand All @@ -139,7 +143,6 @@ class MGLevel : public parAlmond::multigridLevel
void setupSmoother(elliptic_t* base);
dfloat maxEigSmoothAx();

void buildCoarsenerTriTet(mesh_t** meshLevels, int Nf, int Nc);
void buildCoarsenerQuadHex(mesh_t** meshLevels, int Nf, int Nc);
private:
void smoothChebyshevOneIteration (occa::memory &o_r, occa::memory &o_x, bool xIsZero);
Expand Down
51 changes: 11 additions & 40 deletions src/elliptic/ellipticMultiGridLevel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ void MGLevel::Ax(occa::memory o_x, occa::memory o_Ax)

void MGLevel::residual(occa::memory o_rhs, occa::memory o_x, occa::memory o_res)
{
if(stype != SCHWARZ) {
if(stype != SmootherType::SCHWARZ) {
ellipticOperator(elliptic,o_x,o_res, dfloatString);
// subtract r = b - A*x
const dlong Nlocal = mesh->Np * mesh->Nelements;
Expand Down Expand Up @@ -73,24 +73,18 @@ void MGLevel::prolongate(occa::memory o_x, occa::memory o_Px)

void MGLevel::smooth(occa::memory o_rhs, occa::memory o_x, bool x_is_zero)
{
if(!x_is_zero && stype == SCHWARZ) return;
if(!x_is_zero && stype == SmootherType::SCHWARZ) return;
if(!strstr(pfloatString,dfloatString)) {
elliptic->copyDfloatToPfloatKernel(Nrows, o_xPfloat, o_x);
elliptic->copyDfloatToPfloatKernel(Nrows, o_rhsPfloat, o_rhs);
if (stype == RICHARDSON)
this->smoothRichardson(o_rhsPfloat, o_xPfloat, x_is_zero);
else if (stype == CHEBYSHEV)
elliptic->fusedCopyDfloatToPfloatKernel(Nrows, o_x, o_rhs, o_xPfloat, o_rhsPfloat);
if (stype == SmootherType::CHEBYSHEV)
this->smoothChebyshev(o_rhsPfloat, o_xPfloat, x_is_zero);
else if (stype == SCHWARZ)
else
this->smoothSchwarz(o_rhsPfloat, o_xPfloat, x_is_zero);
elliptic->copyPfloatToDPfloatKernel(Nrows, o_xPfloat, o_x);
elliptic->copyPfloatToDPfloatKernel(Nrows, o_rhsPfloat, o_rhs);
} else {
if (stype == RICHARDSON)
this->smoothRichardson(o_rhs, o_x, x_is_zero);
else if (stype == CHEBYSHEV)
if (stype == SmootherType::CHEBYSHEV)
this->smoothChebyshev(o_rhs, o_x, x_is_zero);
else if (stype == SCHWARZ)
else
this->smoothSchwarz(o_rhs, o_x, x_is_zero);
}
}
Expand All @@ -99,41 +93,18 @@ void MGLevel::smoother(occa::memory o_x, occa::memory o_Sx, bool x_is_zero)
{
// x_is_zero = true <-> downward leg
if(x_is_zero) {
if (smtypeDown == JACOBI)
if (smtypeDown == SecondarySmootherType::JACOBI)
this->smootherJacobi(o_x, o_Sx);
else if (smtypeDown == SCHWARZ_SMOOTH)
//this->smootherSchwarz(o_x, o_Sx);
else
this->smoothSchwarz(o_x, o_Sx, true); // no-op if false
} else {
if (smtypeUp == JACOBI)
if (smtypeUp == SecondarySmootherType::JACOBI)
this->smootherJacobi(o_x, o_Sx);
else if (smtypeUp == SCHWARZ_SMOOTH)
//this->smootherSchwarz(o_x, o_Sx);
else
this->smoothSchwarz(o_x, o_Sx, true); // no-op if false
}
}

void MGLevel::smoothRichardson(occa::memory &o_r, occa::memory &o_x, bool xIsZero)
{
occa::memory o_res = o_smootherResidual;

if (xIsZero) {
this->smoother(o_r, o_x, xIsZero);
return;
}

pfloat one = 1.;
pfloat mone = -1.;

//res = r-Ax
this->Ax(o_x,o_res);
elliptic->scaledAddPfloatKernel(Nrows, one, o_r, mone, o_res);

//smooth the fine problem x = x + S(r-Ax)
this->smoother(o_res, o_res, xIsZero);
elliptic->scaledAddPfloatKernel(Nrows, one, o_res, one, o_x);
}

void MGLevel::smoothChebyshevOneIteration (occa::memory &o_r, occa::memory &o_x, bool xIsZero)
{
const pfloat theta = 0.5 * (lambda1 + lambda0);
Expand Down
46 changes: 18 additions & 28 deletions src/elliptic/ellipticMultiGridLevelSetup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,15 +113,15 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)

if (options.compareArgs("MULTIGRID SMOOTHER","ASM") ||
options.compareArgs("MULTIGRID SMOOTHER","RAS")) {
stype = SCHWARZ;
smtypeUp = JACOBI;
smtypeDown = JACOBI;
stype = SmootherType::SCHWARZ;
smtypeUp = SecondarySmootherType::JACOBI;
smtypeDown = SecondarySmootherType::JACOBI;
build(ellipticBase);

if(options.compareArgs("MULTIGRID SMOOTHER","CHEBYSHEV")) {
smtypeUp = SCHWARZ_SMOOTH;
smtypeDown = SCHWARZ_SMOOTH;
stype = CHEBYSHEV;
smtypeUp = SecondarySmootherType::SCHWARZ;
smtypeDown = SecondarySmootherType::SCHWARZ;
stype = SmootherType::CHEBYSHEV;
if (!options.getArgs("MULTIGRID CHEBYSHEV DEGREE", ChebyshevIterations))
ChebyshevIterations = 2; //default to degree 2
//estimate the max eigenvalue of S*A
Expand All @@ -138,13 +138,13 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
casted_invDiagA[i] = static_cast<pfloat>(invDiagA[i]);
o_invDiagA = platform->device.malloc(mesh->Np * mesh->Nelements * sizeof(pfloat), casted_invDiagA.data());
if(options.compareArgs("MULTIGRID UPWARD SMOOTHER","JACOBI"))
smtypeUp = JACOBI;
smtypeUp = SecondarySmootherType::JACOBI;
if(options.compareArgs("MULTIGRID DOWNWARD SMOOTHER","JACOBI"))
smtypeDown = JACOBI;
smtypeDown = SecondarySmootherType::JACOBI;
}
} else if (options.compareArgs("MULTIGRID SMOOTHER","DAMPEDJACOBI")) { //default to damped jacobi
smtypeUp = JACOBI;
smtypeDown = JACOBI;
smtypeUp = SecondarySmootherType::JACOBI;
smtypeDown = SecondarySmootherType::JACOBI;
dfloat* invDiagA;
ellipticBuildJacobi(elliptic,&invDiagA);
std::vector<pfloat> casted_invDiagA(mesh->Np * mesh->Nelements, 0.0);
Expand All @@ -154,7 +154,7 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
o_invDiagA = platform->device.malloc(mesh->Np * mesh->Nelements * sizeof(pfloat), casted_invDiagA.data());

if (options.compareArgs("MULTIGRID SMOOTHER","CHEBYSHEV")) {
stype = CHEBYSHEV;
stype = SmootherType::CHEBYSHEV;

if (!options.getArgs("MULTIGRID CHEBYSHEV DEGREE", ChebyshevIterations))
ChebyshevIterations = 2; //default to degree 2
Expand All @@ -165,19 +165,11 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
lambda1 = 1.1 * rho;
lambda0 = rho / 10.;
}else {
stype = RICHARDSON;

//estimate the max eigenvalue of S*A
dfloat rho = this->maxEigSmoothAx();

//set the stabilty weight (jacobi-type interation)
lambda0 = (4. / 3.) / rho;

for (dlong n = 0; n < mesh->Np * mesh->Nelements; n++)
casted_invDiagA[n] *= static_cast<pfloat>(lambda0);

//update diagonal with weight
o_invDiagA.copyFrom(casted_invDiagA.data());
std::string invalidSmootherName;
options.getArgs("MULTIGRID SMOOTHER", invalidSmootherName);
if(platform->comm.mpiRank == 0) printf("Smoother %s is not supported!\n", invalidSmootherName.c_str());
ABORT(EXIT_FAILURE);
}
free(invDiagA);
}
Expand All @@ -201,13 +193,11 @@ void MGLevel::Report()

char smootherString[BUFSIZ];
if (degree != 1) {
if (stype == RICHARDSON && smtypeDown == JACOBI)
strcpy(smootherString, "Damped Jacobi ");
else if (stype == CHEBYSHEV && smtypeDown == JACOBI)
if (stype == SmootherType::CHEBYSHEV && smtypeDown == SecondarySmootherType::JACOBI)
strcpy(smootherString, "Chebyshev+Jacobi ");
else if (stype == SCHWARZ)
else if (stype == SmootherType::SCHWARZ)
strcpy(smootherString, "Schwarz ");
else if (stype == CHEBYSHEV && smtypeDown == SCHWARZ_SMOOTH)
else if (stype == SmootherType::CHEBYSHEV && smtypeDown == SecondarySmootherType::SCHWARZ)
strcpy(smootherString, "Chebyshev+Schwarz");
}

Expand Down
4 changes: 4 additions & 0 deletions src/elliptic/ellipticSolveSetup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,6 +271,10 @@ void ellipticSolveSetup(elliptic_t* elliptic, occa::properties kernelInfo)
"mask",
pfloatKernelInfo);
filename = install_dir + "/okl/elliptic/ellipticLinAlg.okl";
elliptic->fusedCopyDfloatToPfloatKernel =
platform->device.buildKernel(filename,
"fusedCopyDfloatToPfloat",
kernelInfo);
elliptic->copyDfloatToPfloatKernel =
platform->device.buildKernel(filename,
"copyDfloatToPfloat",
Expand Down