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.

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

WIP: setting up gfloat #70

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
1 change: 1 addition & 0 deletions include/mesh.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,7 @@ class mesh_t {
occa::memory o_sM; // Surface mass

// volume, surface, and second order geometric factors
string gfloatString;
occa::memory o_vgeo, o_sgeo, o_ggeo;

//face node mappings
Expand Down
1 change: 1 addition & 0 deletions libs/mesh/meshOccaSetupHex3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ void meshHex3D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*Np*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nfp*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Np*Nggeo*sizeof(dfloat), ggeo);
gfloatString = dfloatString; // define type of ggeo

/* NC: disabling until we re-add treatment of affine elements

Expand Down
1 change: 1 addition & 0 deletions libs/mesh/meshOccaSetupQuad2D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,4 +41,5 @@ void meshQuad2D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*Np*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nfp*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Np*Nggeo*sizeof(dfloat), ggeo);
gfloatString = dfloatString; // define type of ggeo
}
1 change: 1 addition & 0 deletions libs/mesh/meshOccaSetupQuad3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,4 +41,5 @@ void meshQuad3D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*Np*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nfp*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Np*Nggeo*sizeof(dfloat), ggeo);
gfloatString = dfloatString; // define type of ggeo
}
1 change: 1 addition & 0 deletions libs/mesh/meshOccaSetupTet3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,7 @@ void meshTet3D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Nggeo*sizeof(dfloat), ggeo);
gfloatString = dfloatString; // define type of ggeo

free(DT);
free(LIFTT);
Expand Down
3 changes: 2 additions & 1 deletion libs/mesh/meshOccaSetupTri2D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ void meshTri2D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Nggeo*sizeof(dfloat), ggeo);

gfloatString = dfloatString; // define type of ggeo

free(DT);
free(LIFTT);
free(sMT);
Expand Down
3 changes: 2 additions & 1 deletion libs/mesh/meshOccaSetupTri3D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,8 @@ void meshTri3D::OccaSetup(){
o_vgeo = platform.malloc((Nelements+totalHaloPairs)*Nvgeo*sizeof(dfloat), vgeo);
o_sgeo = platform.malloc(Nelements*Nfaces*Nsgeo*sizeof(dfloat), sgeo);
o_ggeo = platform.malloc(Nelements*Nggeo*sizeof(dfloat), ggeo);

gfloatString = dfloatString; // define type of ggeo

free(DT);
free(LIFTT);
free(sMT);
Expand Down
16 changes: 8 additions & 8 deletions solvers/elliptic/okl/ellipticAxHex3D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ SOFTWARE.


@kernel void ellipticAxHex3D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -163,7 +163,7 @@ SOFTWARE.
@kernel void ellipticPartialAxHex3D_v0(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -308,7 +308,7 @@ SOFTWARE.

@kernel void ellipticPartialAxHex3D_v1(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -1066,7 +1066,7 @@ SOFTWARE.
// SPAM KERNELS
@kernel void ellipticPartialAxHex3D_v2(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -1144,7 +1144,7 @@ SOFTWARE.

@kernel void ellipticPartialAxHex3D_v3(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -1364,7 +1364,7 @@ SOFTWARE.
#if 0
@kernel void ellipticPartialAxHex3D_v4(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -1589,7 +1589,7 @@ SOFTWARE.

@kernel void ellipticPartialAxHex3D_v5(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -1735,7 +1735,7 @@ SOFTWARE.

@kernel void ellipticPartialAxHex3D_v6(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down
4 changes: 2 additions & 2 deletions solvers/elliptic/okl/ellipticAxQuad2D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ SOFTWARE.

// square thread version
@kernel void ellipticAxQuad2D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -133,7 +133,7 @@ SOFTWARE.
@kernel void ellipticPartialAxQuad2D(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * DT,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down
4 changes: 2 additions & 2 deletions solvers/elliptic/okl/ellipticAxQuad3D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@

// square thread version
@kernel void ellipticAxQuad3D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -146,7 +146,7 @@
@kernel void ellipticPartialAxQuad3D(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down
6 changes: 3 additions & 3 deletions solvers/elliptic/okl/ellipticAxTet3D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ SOFTWARE.


@kernel void ellipticAxTet3D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -89,7 +89,7 @@ SOFTWARE.

@kernel void ellipticPartialAxTet3D_v0(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -187,7 +187,7 @@ SOFTWARE.
@kernel void ellipticPartialAxTet3D(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down
4 changes: 2 additions & 2 deletions solvers/elliptic/okl/ellipticAxTri2D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ SOFTWARE.


@kernel void ellipticAxTri2D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -97,7 +97,7 @@ SOFTWARE.
@kernel void ellipticPartialAxTri2D(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * D,
@restrict const dfloat * S,
@restrict const dfloat * MM,
Expand Down
4 changes: 2 additions & 2 deletions solvers/elliptic/okl/ellipticAxTri3D.okl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@


@kernel void ellipticAxTri3D(const dlong Nelements,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * Dmatrices,
@restrict const dfloat * Smatrices,
@restrict const dfloat * MM,
Expand Down Expand Up @@ -97,7 +97,7 @@
@kernel void ellipticPartialAxTri3D(const dlong Nelements,
@restrict const dlong * elementList,
@restrict const dlong * GlobalToLocal,
@restrict const dfloat * ggeo,
@restrict const gfloat * ggeo,
@restrict const dfloat * Dmatrices,
@restrict const dfloat * Smatrices,
@restrict const dfloat * MM,
Expand Down
13 changes: 10 additions & 3 deletions solvers/elliptic/setups/setupHex3D.rc
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ ISOPARAMETRIC
1

[POLYNOMIAL DEGREE]
4
6

[THREAD MODEL]
CUDA
Expand Down Expand Up @@ -68,7 +68,8 @@ MULTIGRID

# can be ALLDEGREES, HALFDEGREES, HALFDOFS
[MULTIGRID COARSENING]
HALFDEGREES
ALLDEGREES
#HALFDEGREES

# can be LOCALPATCH, or DAMPEDJACOBI
# LOCALPATCH smoother can include EXACT
Expand All @@ -80,6 +81,12 @@ CHEBYSHEV
[MULTIGRID CHEBYSHEV DEGREE]
2

# precision of geofacs for p-coarsened multigrid levels
# can be float or double
[MULTIGRID GEOFAC TYPE]
float
#double

###########################################

########## ParAlmond Options ##############
Expand All @@ -104,7 +111,7 @@ CHEBYSHEV
###########################################

[OUTPUT TO FILE]
TRUE
FALSE

[OUTPUT FILE NAME]
elliptic
Expand Down
2 changes: 2 additions & 0 deletions solvers/elliptic/src/ellipticOperator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@ void elliptic_t::Operator(occa::memory &o_q, occa::memory &o_Aq){

ogsMasked->GatheredHaloExchangeStart(o_q, 1, ogs_dfloat);

// printf("ellipticOperator using precision gfloatString=%s \n", mesh.gfloatString.c_str());

if(mesh.NlocalGatherElements){
// if(integrationType==0) { // GLL or non-hex
// if(mapType==0)
Expand Down
49 changes: 49 additions & 0 deletions solvers/elliptic/src/ellipticPreconMultiGrid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,40 @@ void MultiGridPrecon::Operator(occa::memory& o_r, occa::memory& o_Mr) {
if(elliptic.allNeumann) elliptic.ZeroMean(o_Mr);
}

void occaConvertType(platform_t &platform, dlong N, occa::memory &o_x, string &gfloatString){


if(gfloatString==dfloatString) return;

// bring o_x data back to host
dfloat *h_x = (dfloat*) calloc(N, sizeof(dfloat));
o_x.copyTo(h_x);

// free device storage
o_x.free();

// type convert on host and copy back to device
if(gfloatString=="float"){
float *x = (float*) calloc(N, sizeof(float));
for(int n=0;n<N;++n){
x[n] = (float)h_x[n];
}
o_x = platform.device.malloc(N*sizeof(float), x);
free(x);
}

if(gfloatString=="double"){
double *x = (double*) calloc(N, sizeof(double));
for(int n=0;n<N;++n){
x[n] = (double)h_x[n];
}
o_x = platform.device.malloc(N*sizeof(double), x);
free(x);
}

free(h_x);
}

MultiGridPrecon::MultiGridPrecon(elliptic_t& _elliptic):
elliptic(_elliptic), mesh(_elliptic.mesh), settings(_elliptic.settings),
parAlmond(elliptic.platform, settings, mesh.comm) {
Expand All @@ -52,6 +86,21 @@ MultiGridPrecon::MultiGridPrecon(elliptic_t& _elliptic):
while(Nc>1) {
//build mesh and elliptic objects for this degree
mesh_t &meshF = mesh.SetupNewDegree(Nf);

// TW: rewrite o_ggeo here ? and set a variable o_ggeoType in mesh ?
if(Nf<mesh.N){ //
meshF.gfloatString = dfloatString;
settings.getSetting("MULTIGRID GEOFAC TYPE", meshF.gfloatString);

dlong M = 0;
if(meshF.elementType==TRIANGLES || meshF.elementType==TETRAHEDRA)
M = meshF.Nelements*meshF.Nggeo;
else
M = meshF.Nelements*meshF.Np*meshF.Nggeo;

occaConvertType(meshF.platform, M, meshF.o_ggeo, meshF.gfloatString);
}

elliptic_t &ellipticF = elliptic.SetupNewDegree(meshF);

//share masking data with previous MG level
Expand Down
6 changes: 6 additions & 0 deletions solvers/elliptic/src/ellipticSettings.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,11 @@ void ellipticAddSettings(settings_t& settings,
"2",
"Smoothing iterations in Chebyshev smoother");

settings.newSetting(prefix+"MULTIGRID GEOFAC TYPE",
"double",
"Precision of geometric factors to use for p coarsened p-multigrid levels",
{"float", "double"});

settings.newSetting(prefix+"VERBOSE",
"FALSE",
"Enable verbose output",
Expand All @@ -116,6 +121,7 @@ void ellipticSettings_t::report() {
reportSetting("MULTIGRID SMOOTHER");
if (compareSetting("MULTIGRID SMOOTHER","CHEBYSHEV"))
reportSetting("MULTIGRID CHEBYSHEV DEGREE");
reportSetting("MULTIGRID GEOFAC TYPE");
}

if (compareSetting("PRECONDITIONER","MULTIGRID")
Expand Down
3 changes: 2 additions & 1 deletion solvers/elliptic/src/ellipticSetup.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,8 @@ elliptic_t& elliptic_t::Setup(platform_t& platform, mesh_t& mesh,

int NblockV = mymax(1,blockMax/mesh.Np);
kernelInfo["defines/" "p_NblockV"]= NblockV;

kernelInfo["defines/" "gfloat"] = mesh.gfloatString;

// Ax kernel
if (settings.compareSetting("DISCRETIZATION","CONTINUOUS")) {
sprintf(fileName, DELLIPTIC "/okl/ellipticAx%s.okl", suffix);
Expand Down
1 change: 1 addition & 0 deletions solvers/elliptic/src/ellipticSetupNewDegree.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,7 @@ elliptic_t& elliptic_t::SetupNewDegree(mesh_t& meshC){

int NblockV = mymax(1,blockMax/meshC.Np);
kernelInfo["defines/" "p_NblockV"]= NblockV;
kernelInfo["defines/" "gfloat"] = meshC.gfloatString;

// Ax kernel
if (settings.compareSetting("DISCRETIZATION","CONTINUOUS")) {
Expand Down
2 changes: 2 additions & 0 deletions test/testElliptic.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ def ellipticSettings(rcformat="2.0", data_file=ellipticData2D,
linear_solver="PCG",
precon="MULTIGRID",
multigrid_smoother="CHEBYSHEV",
multigrid_geofac_type="float",
paralmond_cycle="VCYCLE",
paralmond_strength="SYMMETRIC",
paralmond_aggregation="UNSMOOTHED",
Expand All @@ -64,6 +65,7 @@ def ellipticSettings(rcformat="2.0", data_file=ellipticData2D,
setting_t("LINEAR SOLVER", linear_solver),
setting_t("PRECONDITIONER", precon),
setting_t("MULTIGRID SMOOTHER", multigrid_smoother),
setting_t("MULTIGRID GEOFAC TYPE", multigrid_geofac_type),
setting_t("PARALMOND CYCLE", paralmond_cycle),
setting_t("PARALMOND STRENGTH", paralmond_strength),
setting_t("PARALMOND AGGREGATION", paralmond_aggregation),
Expand Down