Migration to 3D Mpi complete for backstep

This commit is contained in:
Aditya Ujeniya 2023-10-28 12:04:37 +02:00
parent 175ec58667
commit 120be69102
22 changed files with 22097 additions and 280116 deletions

View File

@ -39,11 +39,15 @@ $(BUILD_DIR)/%.s: %.c
.PHONY: clean distclean tags info asm format
clean:
clean: vis
$(info ===> CLEAN)
@rm -rf $(BUILD_DIR)
@rm -f tags
vis:
$(info ===> REMOVING VIZUALISATION FILES)
@rm -f vtk_files/particle*.vtk
distclean: clean
$(info ===> DIST CLEAN)
@rm -f $(TARGET)

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,82 @@
#==============================================================================
# Laminar Canal Flow
#==============================================================================
# Problem specific Data:
# ---------------------
name backstep # name of flow setup
bcLeft 3 # flags for boundary conditions
bcRight 3 # 1 = no-slip 3 = outflow
bcBottom 1 # 2 = free-slip 4 = periodic
bcTop 1 #
bcFront 1 #
bcBack 1 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
gz 0.0 #
re 5000.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
w_init 0.0 # initial value for velocity in z-direction
p_init 1.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 7.0 # domain size in x-direction
ylength 1.5 # domain size in y-direction
zlength 1.0 # domain size in z-direction
imax 70 # number of interior cells in x-direction
jmax 15 # number of interior cells in y-direction
kmax 10 # number of interior cells in z-direction
# Time Data:
# ---------
te 60.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
omg 0.52
omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
# -----------------------
numberOfParticles 200
startTime 0
injectTimePeriod 1.0
writeTimePeriod 0.2
x1 0.0
y1 0.0
z1 1.0
x2 0.0
y2 4.0
z2 1.0
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 1
xCenter 0.0
yCenter 0.0
zCenter 0.0
xRectLength 2.0
yRectLength 1.0
zRectLength 2.0
circleRadius 1.0
#===============================================================================

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -18,9 +18,9 @@ gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
gz 0.0 #
re 100.0 # Reynolds number
re 5000.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction
u_init 0.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
w_init 0.0 # initial value for velocity in z-direction
p_init 0.0 # initial value for pressure
@ -31,14 +31,14 @@ p_init 0.0 # initial value for pressure
xlength 30.0 # domain size in x-direction
ylength 4.0 # domain size in y-direction
zlength 4.0 # domain size in z-direction
imax 40 # number of interior cells in x-direction
jmax 10 # number of interior cells in y-direction
kmax 10 # number of interior cells in z-direction
imax 100 # number of interior cells in x-direction
jmax 40 # number of interior cells in y-direction
kmax 40 # number of interior cells in z-direction
# Time Data:
# ---------
te 60.0 # final time
te 200.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
@ -47,17 +47,17 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt)
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
rho 0.52
omg 1.3 # relaxation parameter for SOR iteration
omg 0.52
omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
# -----------------------
numberOfParticles 30
startTime 10
numberOfParticles 200
startTime 0
injectTimePeriod 2.0
writeTimePeriod 2.0
writeTimePeriod 1.0
x1 1.0
y1 0.0
@ -65,4 +65,17 @@ z1 1.0
x2 1.0
y2 4.0
z2 1.0
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 1
xCenter 10.0
yCenter 2.0
zCenter 2.0
xRectLength 8.0
yRectLength 2.0
zRectLength 2.0
circleRadius 1.0
#===============================================================================

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -48,7 +48,7 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt)
itermax 1000 # maximal number of pressure iteration in one time step
eps 0.001 # stopping tolerance for pressure iteration
rho 0.5
omg 1.8 # relaxation parameter for SOR iteration
omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
@ -65,4 +65,17 @@ z1 1.0
x2 1.0
y2 4.0
z2 1.0
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 1
xCenter 10.0
yCenter 2.0
zCenter 2.0
xRectLength 8.0
yRectLength 2.0
zRectLength 2.0
circleRadius 1.0
#===============================================================================

View File

@ -0,0 +1,82 @@
#==============================================================================
# Laminar Canal Flow
#==============================================================================
# Problem specific Data:
# ---------------------
name karman # name of flow setup
bcLeft 3 # flags for boundary conditions
bcRight 3 # 1 = no-slip 3 = outflow
bcBottom 1 # 2 = free-slip 4 = periodic
bcTop 1 #
bcFront 1 #
bcBack 1 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
gz 0.0 #
re 100.0 # Reynolds number
u_init 0.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
w_init 0.0 # initial value for velocity in z-direction
p_init 0.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 30.0 # domain size in x-direction
ylength 8.0 # domain size in y-direction
zlength 8.0 # domain size in z-direction
imax 100 # number of interior cells in x-direction
jmax 40 # number of interior cells in y-direction
kmax 40 # number of interior cells in z-direction
# Time Data:
# ---------
te 100.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
omg 0.52
omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
# -----------------------
numberOfParticles 200
startTime 0
injectTimePeriod 2.0
writeTimePeriod 1.0
x1 1.0
y1 0.0
z1 1.0
x2 1.0
y2 4.0
z2 1.0
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 2
xCenter 5.0
yCenter 4.0
zCenter 4.0
xRectLength 8.0
yRectLength 2.0
zRectLength 2.0
circleRadius 1.0
#===============================================================================

View File

@ -101,6 +101,14 @@ static void setupCommunication(Comm* c, Direction direction, int layer)
MPI_DOUBLE,
&c->rbufferTypes[direction]);
MPI_Type_commit(&c->rbufferTypes[direction]);
MPI_Type_create_subarray(NDIMS,
sizes,
subSizes,
starts,
MPI_ORDER_C,
MPI_INT,
&c->rbufferTypesInt[direction]);
MPI_Type_commit(&c->rbufferTypesInt[direction]);
} else if (layer == BULK) {
MPI_Type_create_subarray(NDIMS,
sizes,
@ -110,6 +118,14 @@ static void setupCommunication(Comm* c, Direction direction, int layer)
MPI_DOUBLE,
&c->sbufferTypes[direction]);
MPI_Type_commit(&c->sbufferTypes[direction]);
MPI_Type_create_subarray(NDIMS,
sizes,
subSizes,
starts,
MPI_ORDER_C,
MPI_INT,
&c->sbufferTypesInt[direction]);
MPI_Type_commit(&c->sbufferTypesInt[direction]);
}
}
@ -233,6 +249,22 @@ void commExchange(Comm* c, double* grid)
c->comm);
}
void commExchangeInt(Comm* c, int* grid)
{
int counts[6] = { 1, 1, 1, 1, 1, 1 };
MPI_Aint displs[6] = { 0, 0, 0, 0, 0, 0 };
MPI_Neighbor_alltoallw(grid,
counts,
displs,
c->sbufferTypesInt,
grid,
counts,
displs,
c->rbufferTypesInt,
c->comm);
}
void commShift(Comm* c, double* f, double* g, double* h)
{
MPI_Request requests[6] = { MPI_REQUEST_NULL,

View File

@ -27,6 +27,8 @@ typedef struct {
MPI_Comm comm;
MPI_Datatype sbufferTypes[NDIRS];
MPI_Datatype rbufferTypes[NDIRS];
MPI_Datatype sbufferTypesInt[NDIRS];
MPI_Datatype rbufferTypesInt[NDIRS];
int neighbours[NDIRS];
int coords[NDIMS], dims[NDIMS];
int imaxLocal, jmaxLocal, kmaxLocal;
@ -37,6 +39,8 @@ extern void commInit(Comm* comm, int kmax, int jmax, int imax);
extern void commFree(Comm* comm);
extern void commPrintConfig(Comm*);
extern void commExchange(Comm*, double*);
extern void commExchangeInt(Comm*, int*);
extern void commShift(Comm* c, double* f, double* g, double* h);
extern void commReduction(double* v, int op);
extern int commIsBoundary(Comm* c, Direction direction);

View File

@ -71,13 +71,16 @@ int main(int argc, char** argv)
if (tau > 0.0) computeTimestep(&solver);
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
setObjectBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
// if (nt % 100 == 0) normalizePressure(&solver);
solver_generic[variant - 1](&solver);
adaptUV(&solver);
trace(&particletracer, solver.u, solver.v, solver.w, t);
trace(&particletracer, solver.u, solver.v, solver.w, solver.seg, t);
t += solver.dt;
nt++;

View File

@ -102,6 +102,16 @@ void readParameter(Parameter* param, const char* filename)
PARSE_REAL(y2);
PARSE_REAL(z2);
/* Added obstacle geometry parameters */
PARSE_INT(shape);
PARSE_REAL(xCenter);
PARSE_REAL(yCenter);
PARSE_REAL(zCenter);
PARSE_REAL(xRectLength);
PARSE_REAL(yRectLength);
PARSE_REAL(zRectLength);
PARSE_REAL(circleRadius);
}
}

View File

@ -23,7 +23,9 @@ typedef struct {
double startTime, injectTimePeriod, writeTimePeriod;
double x1, y1, z1, x2, y2, z2;
} Parameter;
int shape;
double xCenter, yCenter, zCenter, xRectLength, yRectLength, zRectLength, circleRadius;} Parameter;
void initParameter(Parameter*);
void readParameter(Parameter*, const char*);

View File

@ -16,6 +16,7 @@
#define U(i, j, k) u[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
#define V(i, j, k) v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
#define W(i, j, k) w[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
#define S(i, j, k) seg[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
static int ts = 0;
@ -66,9 +67,12 @@ void printParticles(ParticleTracer* particletracer)
particletracer->xOffsetEnd, particletracer->yOffsetEnd);
}
}
void injectParticles(ParticleTracer* particletracer)
void injectParticles(ParticleTracer* particletracer,int* restrict seg)
{
double x, y, z;
int imaxLocal = particletracer->imaxLocal;
int jmaxLocal = particletracer->jmaxLocal;
int kmaxLocal = particletracer->kmaxLocal;
for(int i = 0; i < particletracer->numberOfParticles; ++i)
{
x = particletracer->linSpaceLine[i].x;
@ -82,14 +86,21 @@ void injectParticles(ParticleTracer* particletracer)
particletracer->particlePool[particletracer->pointer].x = x;
particletracer->particlePool[particletracer->pointer].y = y;
particletracer->particlePool[particletracer->pointer].z = z;
int i = particletracer->particlePool[particletracer->pointer].x / particletracer->dx;
int j = particletracer->particlePool[particletracer->pointer].y / particletracer->dy;
int k = particletracer->particlePool[particletracer->pointer].z / particletracer->dz;
if(S(i,j,k) == NONE)
{
particletracer->particlePool[particletracer->pointer].flag = true;
++(particletracer->pointer);
++(particletracer->totalParticles);
}
}
}
}
void advanceParticles(ParticleTracer* particletracer, double* restrict u, double* restrict v, double* restrict w,double time)
void advanceParticles(ParticleTracer* particletracer, double* restrict u, double* restrict v, double* restrict w, int* restrict seg, double time)
{
int imax = particletracer->imax;
int jmax = particletracer->jmax;
@ -108,11 +119,11 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
double zlength = particletracer->zlength;
Particle buff[particletracer->size][30];
Particle buff[particletracer->size][100];
for(int i = 0; i < particletracer->size; ++i)
{
for(int j = 0; j < 30; ++j)
for(int j = 0; j < 100; ++j)
{
buff[i][j].x = 0.0;
buff[i][j].y = 0.0;
@ -251,6 +262,16 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
}
}
particletracer->particlePool[i].flag = false;
int i_new = new_x/dx, j_new = new_y/dy, k_new = new_z/dz;
int iOffset = particletracer->xOffset / dx,
jOffset = particletracer->yOffset / dy,
kOffset = particletracer->zOffset / dz;
if(S(i_new - iOffset, j_new - jOffset, k_new - kOffset) != NONE)
{
particletracer->particlePool[i].flag = false;
}
}
}
}
@ -259,21 +280,21 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
{
if(i != particletracer->rank)
{
MPI_Send(buff[i], 30, particletracer->mpi_particle, i, 0, particletracer->comm);
MPI_Send(buff[i], 100, particletracer->mpi_particle, i, 0, particletracer->comm);
}
}
for(int i = 0; i < particletracer->size; ++i)
{
if(i != particletracer->rank)
{
MPI_Recv(buff[i], 30, particletracer->mpi_particle, i, 0, particletracer->comm, MPI_STATUS_IGNORE);
MPI_Recv(buff[i], 100, particletracer->mpi_particle, i, 0, particletracer->comm, MPI_STATUS_IGNORE);
}
}
for(int i = 0; i < particletracer->size; ++i)
{
if(i != particletracer->rank)
{
for(int j = 0; j < 30; ++j)
for(int j = 0; j < 100; ++j)
{
if(buff[i][j].flag == true)
{
@ -458,27 +479,16 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
memcpy(particletracer->coords, solver->comm.coords, sizeof(solver->comm.coords));
particletracer->xLocal = particletracer->imaxLocal * particletracer->dx;
particletracer->yLocal = particletracer->jmaxLocal * particletracer->dy;
particletracer->zLocal = particletracer->kmaxLocal * particletracer->dz;
double xLocal[particletracer->size];
double yLocal[particletracer->size];
double zLocal[particletracer->size];
double offset[6][particletracer->size];
MPI_Allgather(&particletracer->xLocal, 1, MPI_DOUBLE, xLocal, 1, MPI_DOUBLE, particletracer->comm);
MPI_Allgather(&particletracer->yLocal, 1, MPI_DOUBLE, yLocal, 1, MPI_DOUBLE, particletracer->comm);
MPI_Allgather(&particletracer->zLocal, 1, MPI_DOUBLE, zLocal, 1, MPI_DOUBLE, particletracer->comm);
particletracer->xOffset = sumOffset(xLocal, particletracer->rank, (particletracer->dims[1] * particletracer->dims[2]), particletracer->coords[0]);
particletracer->yOffset = sumOffset(yLocal, particletracer->rank, particletracer->dims[2], particletracer->coords[1]);
particletracer->zOffset = sumOffset(zLocal, particletracer->rank, 1, particletracer->coords[2]);
particletracer->xOffsetEnd = particletracer->xOffset + particletracer->xLocal;
particletracer->yOffsetEnd = particletracer->yOffset + particletracer->yLocal;
particletracer->zOffsetEnd = particletracer->zOffset + particletracer->zLocal;
particletracer->xOffset = solver->xOffset;
particletracer->yOffset = solver->yOffset;
particletracer->zOffset = solver->zOffset;
particletracer->xOffsetEnd = solver->xOffsetEnd;
particletracer->yOffsetEnd = solver->yOffsetEnd;
particletracer->zOffsetEnd = solver->zOffsetEnd;
//printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, zOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n", particletracer->rank, particletracer->xOffset, particletracer->yOffset, particletracer->zOffset, particletracer->xOffsetEnd, particletracer->yOffsetEnd, particletracer->zOffsetEnd);
@ -560,13 +570,13 @@ void printParticleTracerParameters(ParticleTracer* particletracer)
printf("\txLocal : %.2f, yLocal : %.2f, zLocal : %.2f\n", particletracer->xLocal, particletracer->yLocal, particletracer->zLocal);
}
void trace(ParticleTracer* particletracer, double* restrict u, double* restrict v, double* restrict w, double time)
void trace(ParticleTracer* particletracer, double* restrict u, double* restrict v, double* restrict w, int* restrict seg,double time)
{
if (time >= particletracer->startTime)
{
if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
{
injectParticles(particletracer);
injectParticles(particletracer, seg);
particletracer->lastInjectTime = time;
}
if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
@ -574,7 +584,7 @@ void trace(ParticleTracer* particletracer, double* restrict u, double* restrict
writeParticles(particletracer);
particletracer->lastWriteTime = time;
}
advanceParticles(particletracer, u, v, w, time);
advanceParticles(particletracer, u, v, w, seg, time);
compress(particletracer);
particletracer->lastUpdateTime = time;
}

View File

@ -50,12 +50,12 @@ typedef struct {
} ParticleTracer;
void initParticleTracer(ParticleTracer*, Parameter*, Solver* );
void injectParticles(ParticleTracer*);
void advanceParticles(ParticleTracer*, double* , double* , double* , double);
void injectParticles(ParticleTracer*, int*);
void advanceParticles(ParticleTracer*, double* , double* , double* ,int*, double);
void freeParticles(ParticleTracer*);
void writeParticles(ParticleTracer*);
void printParticleTracerParameters(ParticleTracer*);
void printParticles(ParticleTracer*);
void trace(ParticleTracer*, double* , double* , double* , double );
void trace(ParticleTracer*, double* , double* , double* , int*, double );
void compress(ParticleTracer* );
#endif

View File

@ -32,6 +32,13 @@
w[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
#define RHS(i, j, k) \
rhs[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
#define S(i, j, k) seg[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
static double distance(
double i, double j, double k, double iCenter, double jCenter, double kCenter)
{
return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) + pow(kCenter - k, 2) * 1.0);
}
static void printConfig(Solver* s)
{
@ -72,6 +79,125 @@ static void printConfig(Solver* s)
commPrintConfig(&s->comm);
}
static double sumOffset(double* sizes, int init, int offset, int coord)
{
double sum = 0;
for (int i = init - offset; coord > 0; i-=offset, --coord)
{
sum += sizes[i];
}
return sum;
}
void printGrid(Solver* solver, int* grid)
{
int imax = solver->grid.imax;
int jmax = solver->grid.jmax;
for (int k = 0; k < solver->grid.kmax + 2; k++) {
printf("K : %02d:\n", k);
for (int j = 0; j < solver->grid.jmax + 2; j++) {
printf("J : %02d: ", j);
for (int i = 0; i < solver->grid.imax + 2; i++) {
switch (grid[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]) {
case FRONTFACE:
printf("FF ");
break;
case BACKFACE:
printf("BF ");
break;
case LEFTFACE:
printf("LF ");
break;
case RIGHTFACE:
printf("RF ");
break;
case TOPFACE:
printf("TF ");
break;
case BOTTOMFACE:
printf("BMF ");
break;
case FRONTTOPLEFTCORNER:
printf("FTLC ");
break;
case FRONTTOPRIGHTCORNER:
printf("FTRC ");
break;
case FRONTBOTTOMLEFTCORNER:
printf("FBLC ");
break;
case FRONTBOTTOMRIGHTCORNER:
printf("FBRC ");
break;
case BACKTOPLEFTCORNER:
printf("BTLC ");
break;
case BACKTOPRIGHTCORNER:
printf("BTRC ");
break;
case BACKBOTTOMLEFTCORNER:
printf("BBLC ");
break;
case BACKBOTTOMRIGHTCORNER:
printf("BBRC ");
break;
case FRONTTOPLINE:
printf("FTL ");
break;
case FRONTBOTTOMLINE:
printf("FBL ");
break;
case FRONTLEFTLINE:
printf("FLL ");
break;
case FRONTRIGHTLINE:
printf("FRL ");
break;
case MIDTOPLEFTLINE:
printf("MTLL ");
break;
case MIDTOPRIGHTLINE:
printf("MTRL ");
break;
case MIDBOTTOMLEFTLINE:
printf("MBTL ");
break;
case MIDBOTTOMRIGHTLINE:
printf("MBRL ");
break;
case BACKTOPLINE:
printf("BTL ");
break;
case BACKBOTTOMLINE:
printf("BBL ");
break;
case BACKLEFTLINE:
printf("BLL ");
break;
case BACKRIGHTLINE:
printf("BRL ");
break;
case LOCAL:
printf("L ");
break;
case NONE:
printf("N ");
break;
}
}
printf("\n");
}
printf("\n\n");
}
fflush(stdout);
}
void initSolver(Solver* s, Parameter* params)
{
s->problem = params->name;
@ -121,6 +247,8 @@ void initSolver(Solver* s, Parameter* params)
s->f = allocate(64, size * sizeof(double));
s->g = allocate(64, size * sizeof(double));
s->h = allocate(64, size * sizeof(double));
s->seg = allocate(64, size * sizeof(int));
for (int i = 0; i < size; i++) {
s->u[i] = params->u_init;
@ -131,6 +259,8 @@ void initSolver(Solver* s, Parameter* params)
s->f[i] = 0.0;
s->g[i] = 0.0;
s->h[i] = 0.0;
s->seg[i] = NONE;
}
double dx = s->grid.dx;
@ -140,6 +270,225 @@ void initSolver(Solver* s, Parameter* params)
double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy) + 1.0 / (dz * dz);
s->dtBound = 0.5 * s->re * 1.0 / invSqrSum;
double xCenter = 0, yCenter = 0, zCenter = 0, radius = 0;
double x1 = 0, x2 = 0, y1 = 0, y2 = 0, z1 = 0, z2 = 0;
int iOffset = 0, jOffset = 0, kOffset = 0;
s->xLocal = s->comm.imaxLocal * s->grid.dx;
s->yLocal = s->comm.jmaxLocal * s->grid.dy;
s->zLocal = s->comm.kmaxLocal * s->grid.dz;
double xLocal[s->comm.size];
double yLocal[s->comm.size];
double zLocal[s->comm.size];
MPI_Allgather(&s->xLocal, 1, MPI_DOUBLE, xLocal, 1, MPI_DOUBLE, s->comm.comm);
MPI_Allgather(&s->yLocal, 1, MPI_DOUBLE, yLocal, 1, MPI_DOUBLE, s->comm.comm);
MPI_Allgather(&s->zLocal, 1, MPI_DOUBLE, zLocal, 1, MPI_DOUBLE, s->comm.comm);
s->xOffset = sumOffset(xLocal, s->comm.rank, (s->comm.dims[1] * s->comm.dims[2]), s->comm.coords[0]);
s->yOffset = sumOffset(yLocal, s->comm.rank, s->comm.dims[2], s->comm.coords[1]);
s->zOffset = sumOffset(zLocal, s->comm.rank, 1, s->comm.coords[2]);
s->xOffsetEnd = s->xOffset + s->xLocal;
s->yOffsetEnd = s->yOffset + s->yLocal;
s->zOffsetEnd = s->zOffset + s->zLocal;
// printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, zOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n", s->comm.rank, s->xOffset, s->yOffset, s->zOffset, s->xOffsetEnd, s->yOffsetEnd, s->zOffsetEnd);
// exit(0);
int* seg = s->seg;
iOffset = s->xOffset / dx;
jOffset = s->yOffset / dy;
kOffset = s->zOffset / dz;
switch (params->shape) {
case NOSHAPE:
break;
case RECT:
x1 = params->xCenter - params->xRectLength / 2;
x2 = params->xCenter + params->xRectLength / 2;
y1 = params->yCenter - params->yRectLength / 2;
y2 = params->yCenter + params->yRectLength / 2;
z1 = params->zCenter - params->zRectLength / 2;
z2 = params->zCenter + params->zRectLength / 2;
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
if ((x1 <= ((i + iOffset) * dx)) && (((i + iOffset) * dx) <= x2) && (y1 <= ((j +jOffset) * dy)) &&
(((j+ jOffset) * dy) <= y2) && ((z1 <= ((k + kOffset) * dz)) && (((k + kOffset)* dz) <= z2))) {
S(i, j, k) = LOCAL;
}
}
}
}
break;
case CIRCLE:
xCenter = params->xCenter;
yCenter = params->yCenter;
zCenter = params->zCenter;
radius = params->circleRadius;
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
if (distance(((i + iOffset) * dx),
((j +jOffset) * dy),
((k + kOffset) * dz),
xCenter,
yCenter,
zCenter) <= radius) {
S(i, j, k) = LOCAL;
}
}
}
}
break;
default:
break;
}
commExchangeInt(&s->comm, seg);
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
/* Assigning enum values to Corners */
if (S(i - 1, j + 1, k - 1) == NONE && S(i - 1, j, k) == NONE &&
S(i, j + 1, k) == NONE && S(i, j, k - 1) == NONE &&
S(i + 1, j - 1, k + 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTTOPLEFTCORNER; //
}
if (S(i + 1, j + 1, k - 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j + 1, k) == NONE && S(i, j, k - 1) == NONE &&
S(i - 1, j - 1, k + 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTTOPRIGHTCORNER; //
}
if (S(i - 1, j - 1, k - 1) == NONE && S(i - 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i, j, k - 1) == NONE &&
S(i + 1, j + 1, k + 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTBOTTOMLEFTCORNER; //
}
if (S(i + 1, j - 1, k - 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i, j, k - 1) == NONE &&
S(i - 1, j + 1, k + 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTBOTTOMRIGHTCORNER; //
}
if (S(i - 1, j + 1, k + 1) == NONE && S(i - 1, j, k) == NONE &&
S(i, j + 1, k) == NONE && S(i, j, k + 1) == NONE &&
S(i + 1, j - 1, k - 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = BACKTOPLEFTCORNER; //
}
if (S(i + 1, j + 1, k + 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j + 1, k) == NONE && S(i, j, k + 1) == NONE &&
S(i - 1, j - 1, k - 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = BACKTOPRIGHTCORNER;
}
if (S(i - 1, j - 1, k + 1) == NONE && S(i - 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i, j, k + 1) == NONE &&
S(i + 1, j + 1, k - 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = BACKBOTTOMLEFTCORNER;
}
if (S(i + 1, j - 1, k + 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i, j, k + 1) == NONE &&
S(i - 1, j + 1, k - 1) == LOCAL && S(i, j, k) == LOCAL) {
S(i, j, k) = BACKBOTTOMRIGHTCORNER;
}
/* Assigning enum values to Lines */
if (S(i - 1, j, k - 1) == NONE && S(i, j, k - 1) == NONE &&
S(i - 1, j, k) == NONE && S(i + 1, j, k + 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTLEFTLINE;
}
if (S(i + 1, j, k - 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j, k - 1) == NONE && S(i - 1, j, k + 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTRIGHTLINE;
}
if (S(i, j + 1, k - 1) == NONE && S(i, j + 1, k) == NONE &&
S(i, j, k - 1) == NONE && S(i, j - 1, k + 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTTOPLINE;
}
if (S(i, j - 1, k - 1) == NONE && S(i, j, k - 1) == NONE &&
S(i, j - 1, k) == NONE && S(i, j + 1, k + 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTBOTTOMLINE;
}
if (S(i - 1, j + 1, k) == NONE && S(i, j + 1, k) == NONE &&
S(i - 1, j, k) == NONE && S(i + 1, j - 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = MIDTOPLEFTLINE;
}
if (S(i + 1, j + 1, k) == NONE && S(i + 1, j, k) == NONE &&
S(i, j + 1, k) == NONE && S(i - 1, j - 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = MIDTOPRIGHTLINE;
}
if (S(i - 1, j - 1, k) == NONE && S(i - 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i + 1, j + 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = MIDBOTTOMLEFTLINE;
}
if (S(i + 1, j - 1, k) == NONE && S(i + 1, j, k) == NONE &&
S(i, j - 1, k) == NONE && S(i - 1, j + 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = MIDBOTTOMRIGHTLINE;
}
if (S(i - 1, j, k + 1) == NONE && S(i - 1, j, k) == NONE &&
S(i, j, k + 1) == NONE && S(i + 1, j, k - 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BACKLEFTLINE;
}
if (S(i + 1, j, k + 1) == NONE && S(i + 1, j, k) == NONE &&
S(i, j, k + 1) == NONE && S(i - 1, j, k - 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BACKRIGHTLINE;
}
if (S(i, j + 1, k + 1) == NONE && S(i, j + 1, k) == NONE &&
S(i, j, k + 1) == NONE && S(i, j - 1, k - 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BACKTOPLINE;
}
if (S(i, j - 1, k + 1) == NONE && S(i, j - 1, k) == NONE &&
S(i, j, k + 1) == NONE && S(i, j + 1, k - 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BACKBOTTOMLINE;
}
/* Assigning enum values to Faces */
if (S(i, j, k - 1) == NONE && S(i, j, k + 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = FRONTFACE; //
}
if (S(i, j, k + 1) == NONE && S(i, j, k - 1) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BACKFACE; //
}
if (S(i, j - 1, k) == NONE && S(i, j + 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = BOTTOMFACE; //
}
if (S(i, j + 1, k) == NONE && S(i, j - 1, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = TOPFACE; //
}
if (S(i - 1, j, k) == NONE && S(i + 1, j, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = LEFTFACE; //
}
if (S(i + 1, j, k) == NONE && S(i - 1, j, k) == LOCAL &&
S(i, j, k) == LOCAL) {
S(i, j, k) = RIGHTFACE; //
}
}
}
}
#ifdef VERBOSE
//printConfig(s);
#endif /* VERBOSE */
@ -160,12 +509,15 @@ void computeRHS(Solver* s)
double* f = s->f;
double* g = s->g;
double* h = s->h;
int* seg = s->seg;
commShift(&s->comm, f, g, h);
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
if(S(i,j,k) == NONE)
{
RHS(i, j, k) = ((F(i, j, k) - F(i - 1, j, k)) * idx +
(G(i, j, k) - G(i, j - 1, k)) * idy +
(H(i, j, k) - H(i, j, k - 1)) * idz) *
@ -173,6 +525,7 @@ void computeRHS(Solver* s)
}
}
}
}
}
void solveRB(Solver* s)
@ -202,6 +555,7 @@ void solveRB(Solver* s)
int it = 0;
double res = 1.0;
int pass, ksw, jsw, isw;
int* seg = s->seg;
while ((res >= epssq) && (it < itermax)) {
@ -217,7 +571,8 @@ void solveRB(Solver* s)
isw = jsw;
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imaxLocal + 1; i += 2) {
if(S(i,j,k) == NONE)
{
double r =
RHS(i, j, k) -
((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
@ -228,6 +583,7 @@ void solveRB(Solver* s)
P(i, j, k) -= (factor * r);
res += (r * r);
}
}
isw = 3 - isw;
}
@ -283,6 +639,7 @@ void solveRB(Solver* s)
}
}
}
setObjectPBoundaryCondition(s);
commReduction(&res, SUM);
res = res / (double)(imax * jmax * kmax);
@ -329,6 +686,7 @@ void solve(Solver* s)
double epssq = eps * eps;
int it = 0;
double res = 1.0;
int* seg = s->seg;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
@ -338,7 +696,8 @@ void solve(Solver* s)
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
if(S(i,j,k) == NONE)
{
double r = RHS(i, j, k) -
((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) *
idx2 +
@ -352,6 +711,7 @@ void solve(Solver* s)
}
}
}
}
if (commIsBoundary(&s->comm, FRONT)) {
for (int j = 1; j < jmaxLocal + 1; j++) {
@ -400,6 +760,7 @@ void solve(Solver* s)
}
}
}
setObjectPBoundaryCondition(s);
commReduction(&res, SUM);
res = res / (double)(imax * jmax * kmax);
@ -449,6 +810,7 @@ void solveRBA(Solver* s)
int it = 0;
double res = 1.0;
int pass, ksw, jsw, isw;
int* seg = s->seg;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
@ -463,7 +825,8 @@ void solveRBA(Solver* s)
isw = jsw;
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imaxLocal + 1; i += 2) {
if(S(i,j,k) == NONE)
{
double r =
RHS(i, j, k) -
((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
@ -475,6 +838,7 @@ void solveRBA(Solver* s)
P(i, j, k) -= (omega * factor * r);
res += (r * r);
}
}
isw = 3 - isw;
}
jsw = 3 - jsw;
@ -531,6 +895,7 @@ void solveRBA(Solver* s)
}
}
}
setObjectPBoundaryCondition(s);
commReduction(&res, SUM);
res = res / (double)(imax * jmax * kmax);
@ -870,6 +1235,7 @@ void computeFG(Solver* s)
double* f = s->f;
double* g = s->g;
double* h = s->h;
int* seg = s->seg;
double gx = s->gx;
double gy = s->gy;
@ -894,6 +1260,8 @@ void computeFG(Solver* s)
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
if(S(i,j,k) == NONE)
{
du2dx = inverseDx * 0.25 *
((U(i, j, k) + U(i + 1, j, k)) *
(U(i, j, k) + U(i + 1, j, k)) -
@ -1020,6 +1388,117 @@ void computeFG(Solver* s)
H(i, j, k) = W(i, j, k) + dt * (inverseRe * (dw2dx2 + dw2dy2 + dw2dz2) -
duwdx - dvwdy - dw2dz + gz);
}
else{
switch (S(i, j, k)) {
case TOPFACE:
G(i, j, k) = V(i, j, k);
break;
case BOTTOMFACE:
G(i, j, k) = V(i, j, k);
break;
case LEFTFACE:
F(i, j, k) = U(i, j, k);
break;
case RIGHTFACE:
F(i, j, k) = U(i, j, k);
break;
case FRONTFACE:
H(i, j, k) = W(i, j, k);
break;
case BACKFACE:
H(i, j, k) = W(i, j, k);
break;
case FRONTLEFTLINE:
F(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTRIGHTLINE:
F(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTTOPLINE:
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTBOTTOMLINE:
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case MIDTOPLEFTLINE:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
break;
case MIDTOPRIGHTLINE:
F(i, j, k) = U(i, j, k);
G(i, j, k) = V(i, j, k);
break;
case MIDBOTTOMLEFTLINE:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
break;
case MIDBOTTOMRIGHTLINE:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
break;
case BACKLEFTLINE:
F(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKRIGHTLINE:
F(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKTOPLINE:
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKBOTTOMLINE:
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTTOPLEFTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTTOPRIGHTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTBOTTOMLEFTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case FRONTBOTTOMRIGHTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKTOPLEFTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKTOPRIGHTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKBOTTOMLEFTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case BACKBOTTOMRIGHTCORNER:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
}
}
}
}
}
@ -1106,3 +1585,249 @@ void adaptUV(Solver* s)
}
}
}
void setObjectBoundaryCondition(Solver* s)
{
int imaxLocal = s->comm.imaxLocal;
int jmaxLocal = s->comm.jmaxLocal;
int kmaxLocal = s->comm.kmaxLocal;
double* u = s->u;
double* v = s->v;
double* w = s->w;
int* seg = s->seg;
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
switch (S(i, j, k)) {
case TOPFACE:
U(i, j, k) = -U(i, j+1, k);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j+1, k);;
break;
case BOTTOMFACE:
U(i, j, k) = -U(i, j-1, k);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j-1, k);
break;
case LEFTFACE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = -W(i-1, j, k);
break;
case RIGHTFACE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = -W(i+1, j, k);
break;
case FRONTFACE:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = -V(i, j, k-1);
W(i, j, k) = 0.0;
break;
case BACKFACE:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = -V(i, j, k+1);
W(i, j, k) = 0.0;
break;
case FRONTLEFTLINE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i, j, k-1);
W(i, j, k) = -W(i-1, j, k);
break;
case FRONTRIGHTLINE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i, j, k-1);
W(i, j, k) = -W(i+1, j, k);
break;
case FRONTTOPLINE:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j+1, k);
break;
case FRONTBOTTOMLINE:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j-1, k);
break;
case MIDTOPLEFTLINE:
U(i, j, k) = -U(i, j+1, k);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = 0.0;
break;
case MIDTOPRIGHTLINE:
U(i, j, k) = 0.0;
V(i, j, k) = 0.0;
U(i-1, j, k) = -U(i-1, j+1, k);
V(i, j-1, k) = -V(i+1, j-1, k);
W(i, j, k) = 0.0;
break;
case MIDBOTTOMLEFTLINE:
U(i, j, k) = -U(i, j-1, k);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = 0.0;
break;
case MIDBOTTOMRIGHTLINE:
U(i, j, k) = -U(i, j-1, k);
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = 0.0;
break;
case BACKLEFTLINE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i, j, k+1);
W(i, j, k) = -W(i-1, j, k);
break;
case BACKRIGHTLINE:
U(i, j, k) = 0.0;
V(i, j, k) = -V(i, j, k+1);
W(i, j, k) = -W(i+1, j, k);
break;
case BACKTOPLINE:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j+1, k);
break;
case BACKBOTTOMLINE:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = 0.0;
W(i, j, k) = -W(i, j-1, k);
break;
case FRONTTOPLEFTCORNER:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = -W(i, j+1, k);
break;
case FRONTTOPRIGHTCORNER:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = -W(i, j+1, k);
break;
case FRONTBOTTOMLEFTCORNER:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = -W(i, j-1, k);
break;
case FRONTBOTTOMRIGHTCORNER:
U(i, j, k) = -U(i, j, k-1);
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = -W(i, j-1, k);
break;
case BACKTOPLEFTCORNER:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = -W(i, j+1, k);
break;
case BACKTOPRIGHTCORNER:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = -W(i, j+1, k);
break;
case BACKBOTTOMLEFTCORNER:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = -V(i-1, j, k);
W(i, j, k) = -W(i, j-1, k);
break;
case BACKBOTTOMRIGHTCORNER:
U(i, j, k) = -U(i, j, k+1);
V(i, j, k) = -V(i+1, j, k);
W(i, j, k) = -W(i, j-1, k);
break;
}
}
}
}
}
void setObjectPBoundaryCondition(Solver* s)
{
int imaxLocal = s->comm.imaxLocal;
int jmaxLocal = s->comm.jmaxLocal;
int kmaxLocal = s->comm.kmaxLocal;
double* p = s->p;
int* seg = s->seg;
for (int k = 1; k < kmaxLocal + 1; k++) {
for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imaxLocal + 1; i++) {
switch (S(i, j, k)) {
case TOPFACE:
P(i,j,k) = P(i,j+1,k);
break;
case BOTTOMFACE:
P(i,j,k) = P(i,j-1,k);
break;
case LEFTFACE:
P(i,j,k) = P(i-1,j,k);
break;
case RIGHTFACE:
P(i,j,k) = P(i+1,j,k);
break;
case FRONTFACE:
P(i,j,k) = P(i,j,k-1);
break;
case BACKFACE:
P(i,j,k) = P(i,j,k+1);
break;
case FRONTLEFTLINE:
P(i,j,k) = (P(i,j,k-1) + P(i-1,j,k)) / 2;
break;
case FRONTRIGHTLINE:
P(i,j,k) = (P(i,j,k-1) + P(i+1,j,k)) / 2;
break;
case FRONTTOPLINE:
P(i,j,k) = (P(i,j,k-1) + P(i,j+1,k)) / 2;
break;
case FRONTBOTTOMLINE:
P(i,j,k) = (P(i,j,k-1) + P(i,j-1,k)) / 2;
break;
case MIDTOPLEFTLINE:
P(i,j,k) = (P(i-1,j,k) + P(i,j+1,k)) / 2;
break;
case MIDTOPRIGHTLINE:
P(i,j,k) = (P(i+1,j,k) + P(i,j+1,k)) / 2;
break;
case MIDBOTTOMLEFTLINE:
P(i,j,k) = (P(i-1,j,k) + P(i,j-1,k)) / 2;
break;
case MIDBOTTOMRIGHTLINE:
P(i,j,k) = (P(i+1,j,k) + P(i,j-1,k)) / 2;
break;
case BACKLEFTLINE:
P(i,j,k) = (P(i,j,k+1) + P(i-1,j,k)) / 2;
break;
case BACKRIGHTLINE:
P(i,j,k) = (P(i,j,k+1) + P(i+1,j,k)) / 2;
break;
case BACKTOPLINE:
P(i,j,k) = (P(i,j,k+1) + P(i,j+1,k)) / 2;
break;
case BACKBOTTOMLINE:
P(i,j,k) = (P(i,j,k+1) + P(i,j-1,k)) / 2;
break;
case FRONTTOPLEFTCORNER:
P(i,j,k) = (P(i,j,k-1) + P(i-1,j,k) + P(i,j+1,k)) / 3;
break;
case FRONTTOPRIGHTCORNER:
P(i,j,k) = (P(i,j,k-1) + P(i+1,j,k) + P(i,j+1,k)) / 3;
break;
case FRONTBOTTOMLEFTCORNER:
P(i,j,k) = (P(i,j,k-1) + P(i-1,j,k) + P(i,j-1,k)) / 3;
break;
case FRONTBOTTOMRIGHTCORNER:
P(i,j,k) = (P(i,j,k-1) + P(i+1,j,k) + P(i,j-1,k)) / 3;
break;
case BACKTOPLEFTCORNER:
P(i,j,k) = (P(i,j,k+1) + P(i-1,j,k) + P(i,j+1,k)) / 3;
break;
case BACKTOPRIGHTCORNER:
P(i,j,k) = (P(i,j,k+1) + P(i+1,j,k) + P(i,j+1,k)) / 3;
break;
case BACKBOTTOMLEFTCORNER:
P(i,j,k) = (P(i,j,k+1) + P(i-1,j,k) + P(i,j-1,k)) / 3;
break;
case BACKBOTTOMRIGHTCORNER:
P(i,j,k) = (P(i,j,k+1) + P(i+1,j,k) + P(i,j-1,k)) / 3;
break;
}
}
}
}
}

View File

@ -10,7 +10,46 @@
#include "grid.h"
#include "parameter.h"
enum OBJECTBOUNDARY {
NONE = 0,
/* Front Corners */
FRONTTOPLEFTCORNER,
FRONTTOPRIGHTCORNER,
FRONTBOTTOMLEFTCORNER,
FRONTBOTTOMRIGHTCORNER,
/* Back Corners */
BACKTOPLEFTCORNER,
BACKTOPRIGHTCORNER,
BACKBOTTOMLEFTCORNER,
BACKBOTTOMRIGHTCORNER,
/* Faces */
FRONTFACE,
BACKFACE,
LEFTFACE,
RIGHTFACE,
TOPFACE,
BOTTOMFACE,
/* Front Lines remaining after Corners and Faces */
FRONTLEFTLINE,
FRONTRIGHTLINE,
FRONTTOPLINE,
FRONTBOTTOMLINE,
/* Bottom Lines remaining after Corners and Faces */
BACKLEFTLINE,
BACKRIGHTLINE,
BACKTOPLINE,
BACKBOTTOMLINE,
/* Mid Lines remaining after Corners and Faces */
MIDTOPLEFTLINE,
MIDTOPRIGHTLINE,
MIDBOTTOMLEFTLINE,
MIDBOTTOMRIGHTLINE,
/* Local where its an object but not a boundary */
LOCAL
};
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
typedef struct {
/* geometry and grid information */
@ -19,6 +58,8 @@ typedef struct {
double *p, *rhs;
double *f, *g, *h;
double *u, *v, *w;
int* seg;
double xLocal, yLocal, zLocal, xOffset, yOffset, zOffset, xOffsetEnd, yOffsetEnd, zOffsetEnd;
/* parameters */
double eps, omega, rho;
double re, tau, gamma;
@ -33,15 +74,18 @@ typedef struct {
Comm comm;
} Solver;
void initSolver(Solver*, Parameter*);
void computeRHS(Solver*);
void solve(Solver*);
void solveRB(Solver*);
void solveRBA(Solver*);
void normalizePressure(Solver*);
void computeTimestep(Solver*);
void setBoundaryConditions(Solver*);
void setSpecialBoundaryCondition(Solver*);
void computeFG(Solver*);
void adaptUV(Solver*);
extern void initSolver(Solver*, Parameter*);
extern void computeRHS(Solver*);
extern void solve(Solver*);
extern void solveRB(Solver*);
extern void solveRBA(Solver*);
extern void normalizePressure(Solver*);
extern void computeTimestep(Solver*);
extern void setBoundaryConditions(Solver*);
extern void setSpecialBoundaryCondition(Solver*);
extern void setObjectBoundaryCondition(Solver*);
extern void setObjectPBoundaryCondition(Solver*);
extern void computeFG(Solver*);
extern void adaptUV(Solver*);
#endif

View File

@ -179,6 +179,7 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
particletracer->particlePool[i].flag = false;
}
int i_new = new_x/dx, j_new = new_y/dy, k_new = new_z/dz;
if(S(i_new, j_new, k_new) != NONE)
{
particletracer->particlePool[i].flag = false;