moebiusband/NuSiF-Solver
2
0
Fork 2
Code Issues Pull Requests 1 Packages Projects Releases Wiki Activity

Porting working changes from 3D-seq to 3D-MPI

Browse Source
This commit is contained in:
Aditya Ujeniya
2023-12-04 14:40:02 +01:00
parent 1987056dce
commit 1b3129d3cb
25 changed files with 5162240 additions and 362320 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
BasicSolver/3D-seq/src/main.c
View File

@@ -101,7 +101,7 @@ int main(int argc, char** argv)
double t = 0.0;
// printGrid(&s, s.seg);
// exit(0);
// exit(0)
timeStart = getTimeStamp();

44
BasicSolver/3D-seq/src/particletracing.c
View File

@@ -34,22 +34,18 @@ void injectParticles(ParticleTracer* particletracer, int* seg)
int imax = particletracer->imax;
int jmax = particletracer->jmax;
int kmax = particletracer->kmax;
for (int i = 0; i < particletracer->numberOfParticles; ++i) {
// particletracer->particlePool[particletracer->pointer].x =
// particletracer->linSpaceLine[i].x;
// particletracer->particlePool[particletracer->pointer].y =
// particletracer->linSpaceLine[i].y;
// particletracer->particlePool[particletracer->pointer].z =
// particletracer->linSpaceLine[i].z;
particletracer->particlePool[particletracer->pointer].x = particletracer->x1;
particletracer->particlePool[particletracer->pointer].y = (double)rand() /
RAND_MAX *
particletracer->ylength;
particletracer->particlePool[particletracer->pointer].z = (double)rand() /
RAND_MAX *
particletracer->zlength;
particletracer->particlePool[particletracer->pointer].y = (((double)rand() /
RAND_MAX) *
(particletracer->y2 - particletracer->y1)) +
particletracer->y1;
particletracer->particlePool[particletracer->pointer].z = (((double)rand() /
RAND_MAX) *
(particletracer->z2 - particletracer->z1)) +
particletracer->z1;
int i = particletracer->particlePool[particletracer->pointer].x /
particletracer->dx;
@@ -58,11 +54,13 @@ void injectParticles(ParticleTracer* particletracer, int* seg)
int k = particletracer->particlePool[particletracer->pointer].z /
particletracer->dz;
if (S(i, j, k) == NONE) {
if(S(i+1, j, k) == NONE)
{
particletracer->particlePool[particletracer->pointer].flag = true;
++(particletracer->pointer);
++(particletracer->totalParticles);
}
else particletracer->particlePool[particletracer->pointer].flag = false;
}
}
@@ -254,24 +252,6 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
sizeof(Particle) * particletracer->estimatedNumParticles);
particletracer->linSpaceLine = malloc(
sizeof(Particle) * particletracer->numberOfParticles);
for (int i = 0; i < particletracer->numberOfParticles; ++i) {
// double spacing = (double)i /
// (double)(particletracer->numberOfParticles - 1);
// particletracer->linSpaceLine[i].x = spacing * particletracer->x1 + (1.0 -
// spacing) * particletracer->x2; particletracer->linSpaceLine[i].y = spacing *
// particletracer->y1 + (1.0 - spacing) * particletracer->y2;
// particletracer->linSpaceLine[i].z = spacing * particletracer->z1 + (1.0 -
// spacing) * particletracer->z2;
particletracer->linSpaceLine[i].x = particletracer->x1;
particletracer->linSpaceLine[i].y = (double)rand() / RAND_MAX *
particletracer->ylength;
particletracer->linSpaceLine[i].z = (double)rand() / RAND_MAX *
particletracer->zlength;
particletracer->linSpaceLine[i].flag = true;
}
}
void printParticleTracerParameters(ParticleTracer* particletracer)

121
BasicSolver/3D-seq/src/solver.c
View File

@@ -443,12 +443,10 @@ void computeRHS(Solver* s)
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
for (int i = 1; i < imax + 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) *
idt;
}
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) *
idt;
}
}
}
@@ -482,18 +480,18 @@ void solve(Solver* s)
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
for (int i = 1; i < imax + 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
// 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
P(i, j, k) -= (factor * r);
res += (r * r);
}
P(i, j, k) -= (factor * r);
res += (r * r);
// }
}
}
}
@@ -518,7 +516,7 @@ void solve(Solver* s)
P(imax + 1, j, k) = P(imax, j, k);
}
}
setObjectPBoundaryCondition(s);
// setObjectPBoundaryCondition(s);
res = res / (double)(imax * jmax * kmax);
#ifdef DEBUG
@@ -580,7 +578,7 @@ void solveRB(Solver* s)
}
}
setObjectPBoundaryCondition(s);
// setObjectPBoundaryCondition(s);
for (pass = 0; pass < 2; pass++) {
jsw = ksw;
@@ -589,19 +587,18 @@ void solveRB(Solver* s)
isw = jsw;
for (int j = 1; j < jmax + 1; j++) {
for (int i = isw; i < imax + 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
// 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
P(i, j, k) -= (factor * r);
res += (r * r);
}
P(i, j, k) -= (factor * r);
res += (r * r);
// }
}
isw = 3 - isw;
}
@@ -657,19 +654,18 @@ void solveRBA(Solver* s)
isw = jsw;
for (int j = 1; j < jmax + 1; j++) {
for (int i = isw; i < imax + 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
// 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 +
(P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
idy2 +
(P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
idz2);
P(i, j, k) -= (omega * factor * r);
res += (r * r);
}
P(i, j, k) -= (omega * factor * r);
res += (r * r);
// }
}
isw = 3 - isw;
}
@@ -700,7 +696,7 @@ void solveRBA(Solver* s)
P(imax + 1, j, k) = P(imax, j, k);
}
}
setObjectPBoundaryCondition(s);
// setObjectPBoundaryCondition(s);
res = res / (double)(imax * jmax * kmax);
#ifdef DEBUG
@@ -990,13 +986,13 @@ void setSpecialBoundaryCondition(Solver* s)
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
y = mDy * (j - 0.5);
U(0, j, k) = y * (ylength - y) * 4.0 / (ylength * ylength);
U(1, j, k) = y * (ylength - y) * 4.0 / (ylength * ylength);
}
}
} else if (strcmp(s->problem, "karman") == 0) {
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
U(0, j, k) = 1.0;
U(1, j, k) = 1.0;
}
}
} else if (strcmp(s->problem, "backstep") == 0) {
@@ -1004,7 +1000,11 @@ void setSpecialBoundaryCondition(Solver* s)
int* seg = s->seg;
for (int k = 1; k < kmax + 1; k++) {
for (int j = 1; j < jmax + 1; j++) {
if (S(0, j, k) == NONE) U(0, j, k) = 1.0;
if (S(1, j, k) == NONE) U(0, j, k) = 1.0;
else {
U(0, j, k) = 0.0;
U(1, j, k) = 0.0;
}
}
}
}
@@ -1176,22 +1176,22 @@ void computeFG(Solver* s)
} else {
switch (S(i, j, k)) {
case TOPFACE:
G(i, j, k) = V(i, j, k);
G(i, j, k) = 0.0;
break;
case BOTTOMFACE:
G(i, j, k) = V(i, j, k);
G(i, j, k) = 0.0;
break;
case LEFTFACE:
F(i, j, k) = U(i, j, k);
F(i, j, k) = 0.0;
break;
case RIGHTFACE:
F(i, j, k) = U(i, j, k);
F(i, j, k) = 0.0;
break;
case FRONTFACE:
H(i, j, k) = W(i, j, k);
H(i, j, k) = 0.0;
break;
case BACKFACE:
H(i, j, k) = W(i, j, k);
H(i, j, k) = 0.0;
break;
case FRONTLEFTLINE:
F(i, j, k) = 0.0;
@@ -1214,8 +1214,8 @@ void computeFG(Solver* s)
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);
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
break;
case MIDBOTTOMLEFTLINE:
F(i, j, k) = 0.0;
@@ -1281,6 +1281,11 @@ void computeFG(Solver* s)
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
case LOCAL:
F(i, j, k) = 0.0;
G(i, j, k) = 0.0;
H(i, j, k) = 0.0;
break;
}
}
}
@@ -1361,7 +1366,7 @@ void setObjectBoundaryCondition(Solver* s)
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);
@@ -1490,6 +1495,10 @@ void setObjectBoundaryCondition(Solver* s)
V(i, j, k) = -V(i + 1, j, k);
W(i, j, k) = -W(i, j - 1, k);
break;
case LOCAL:
U(i, j, k) = 0.0;
V(i, j, k) = 0.0;
W(i, j, k) = 0.0;
}
}
}

4
BasicSolver/3D-seq/src/solver.h
View File

@@ -45,7 +45,9 @@ enum OBJECTBOUNDARY {
MIDBOTTOMLEFTLINE,
MIDBOTTOMRIGHTLINE,
/* Local where its an object but not a boundary */
LOCAL
LOCAL,
/*Ghost cells boundary */
OUTSIDEBOUNDARY
};
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
/// @brief

2
BasicSolver/3D-seq/src/vtkWriter.c
View File

@@ -105,7 +105,7 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
for (int k = 0; k < kmax; k++) {
for (int j = 0; j < jmax; j++) {
for (int i = 0; i < imax; i++) {
if (o->fmt == ASCII) {
if (o->fmt == ASCII /*&& k >= 20*/) {
fprintf(o->fh,
"%f %f %f\n",
G(vec.u, i, j, k),