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

Sample comit

Browse Source
This commit is contained in:
Aditya Ujeniya
2023-09-12 14:10:47 +02:00
parent 1ab581f274
commit eb14891126
101 changed files with 351591 additions and 3215 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -57,6 +57,10 @@ int main (int argc, char** argv)
void (*solver_generic[])() = {solve, solveRB, solveRBA};
//print(&solver, solver.s);
//exit(0);
S = getTimeStamp();
while (t <= te)
@@ -71,7 +75,7 @@ int main (int argc, char** argv)
adaptUV(&solver);
/* Added function for particle tracing. Will inject and advance particles as per timePeriod */
trace(&particletracer, solver.u, solver.v, t);
//trace(&particletracer, solver.u, solver.v, t);
t += solver.dt;
nt++;
@@ -84,6 +88,9 @@ int main (int argc, char** argv)
#endif
}
printf("Total particles : %d\n", particletracer.totalParticles);
//print(&solver, solver.p);
E = getTimeStamp();
stopProgress();

9
BasicSolver/2D-seq/src/parameter.c
View File

@@ -90,6 +90,15 @@ void readParameter(Parameter* param, const char* filename)
PARSE_REAL(y1);
PARSE_REAL(x2);
PARSE_REAL(y2);
/* Added obstacle geometry parameters */
PARSE_INT(shape);
PARSE_REAL(xCenter);
PARSE_REAL(yCenter);
PARSE_REAL(xRectLength);
PARSE_REAL(yRectLength);
PARSE_REAL(circleRadius);
}
}

4
BasicSolver/2D-seq/src/parameter.h
View File

@@ -23,6 +23,10 @@ typedef struct {
double startTime, injectTimePeriod, writeTimePeriod;
double x1, y1, x2, y2;
int shape;
double xCenter, yCenter, xRectLength, yRectLength, circleRadius;
} Parameter;
void initParameter(Parameter*);

1
BasicSolver/2D-seq/src/particletracing.c
View File

@@ -205,7 +205,6 @@ void trace(ParticleTracer* particletracer, double* u, double* v, double time)
if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
{
injectParticles(particletracer);
printf("Rank : %d particles after particles injected !", particletracer->rank);
particletracer->lastInjectTime = time;
}
if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)

2
BasicSolver/2D-seq/src/particletracing.h
View File

@@ -43,6 +43,6 @@ void freeParticles(ParticleTracer*);
void writeParticles(ParticleTracer*);
void printParticleTracerParameters(ParticleTracer*);
void printParticles(ParticleTracer*);
void trace(ParticleTracer*, Solver*, double* , double* , double );
void trace(ParticleTracer*, double* , double* , double );
void compress(ParticleTracer* );
#endif

101
BasicSolver/2D-seq/src/solver.c
View File

@@ -20,9 +20,16 @@
#define G(i, j) g[(j) * (imax + 2) + (i)]
#define U(i, j) u[(j) * (imax + 2) + (i)]
#define V(i, j) v[(j) * (imax + 2) + (i)]
#define S(i, j) s[(j) * (imax + 2) + (i)]
#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
static void print(Solver* solver, double* grid)
static double distance(i, j, iCenter, jCenter)
{
return sqrt(pow(iCenter - i, 2)
+ pow(jCenter - j, 2) * 1.0);
}
void print(Solver* solver, double* grid)
{
int imax = solver->imax;
@@ -91,6 +98,7 @@ void initSolver(Solver* solver, Parameter* params)
size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
solver->u = allocate(64, size);
solver->v = allocate(64, size);
solver->s = allocate(64, size);
solver->p = allocate(64, size);
solver->rhs = allocate(64, size);
solver->f = allocate(64, size);
@@ -103,12 +111,71 @@ void initSolver(Solver* solver, Parameter* params)
solver->rhs[i] = 0.0;
solver->f[i] = 0.0;
solver->g[i] = 0.0;
solver->s[i] = 1.0;
}
double dx = solver->dx;
double dy = solver->dy;
double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum;
int iCenter = 0, jCenter = 0, iRectLength = 0, jRectLength = 0, radius = 0;
double* s = solver->s;
switch(params->shape)
{
case NOSHAPE:
break;
case RECT:
iCenter = params->xCenter/dx;
jCenter = params->yCenter/dy;
iRectLength = params->xRectLength/dx;
jRectLength = params->yRectLength/dy;
int i1 = iCenter - iRectLength/2;
int i2 = iCenter + iRectLength/2;
int j1 = jCenter - jRectLength/2;
int j2 = jCenter + jRectLength/2;
printf("xCenter : %f, yCenter : %f, iCenter : %d, jCenter : %d\n", params->xCenter, params->yCenter, iCenter, jCenter);
printf("xRectLength : %f, yRectLength : %f, dx : %f, dy : %f\n", params->xRectLength, params->yRectLength, dx, dy);
printf("iRectLength : %d, jRectLength : %d, i1 : %d, i2 : %d, j1 : %d, j2 : %d\n", iRectLength, jRectLength, i1, i2, j1, j2);
for (int j = 1; j < jmax + 1; j++)
{
for (int i = 1; i < imax + 1; i++)
{
if(i1 <= i && i <= i2 && j1 <= j && j <= j2)
{
// if(i == i1 || i == i2 || j == j1 || j == j2) S(i, j) = 2.0;
S(i, j) = 0.0;
}
}
}
break;
case CIRCLE:
iCenter = params->xCenter/dx;
jCenter = params->yCenter/dy;
radius = params->circleRadius/dx; // have to change the logic to get max out of both the dx and dy
for (int j = 1; j < jmax + 1; j++)
{
for (int i = 1; i < imax + 1; i++)
{
if(distance(i, j, iCenter, jCenter) <= radius)
{
S(i, j) = 0.0;
}
}
}
break;
default:
break;
}
#ifdef VERBOSE
printConfig(solver);
#endif
@@ -199,10 +266,12 @@ void solveRB(Solver* solver)
double idy2 = 1.0 / dy2;
double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double* p = solver->p;
double* s = solver->s;
double* rhs = solver->rhs;
double epssq = eps * eps;
int it = 0;
double res = 1.0;
double val = 1.0;
int pass, jsw, isw;
while ((res >= epssq) && (it < itermax)) {
@@ -214,12 +283,16 @@ void solveRB(Solver* solver)
for (int j = 1; j < jmax + 1; j++) {
for (int i = isw; i < imax + 1; i += 2) {
// val = S(i+1, j) + S(i-1, j) + S(i, j+1) + S(i, j-1);
// if(val == 0 || S(i,j) == 0)
// {
// continue;
// }
double r = RHS(i, j) -
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
(P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
P(i, j) -= (factor * r);
P(i, j) -= (factor * r );
res += (r * r);
}
isw = 3 - isw;
@@ -568,13 +641,31 @@ void adaptUV(Solver* solver)
double* p = solver->p;
double* u = solver->u;
double* v = solver->v;
double* s = solver->s;
double* f = solver->f;
double* g = solver->g;
double factorX = solver->dt / solver->dx;
double factorY = solver->dt / solver->dy;
for (int j = 1; j < jmax + 1; j++) {
for (int i = 1; i < imax + 1; i++) {
double val = 0;
for (int j = 1; j < jmax + 1; j++)
{
for (int i = 1; i < imax + 1; i++)
{
val = S(i+1, j) + S(i-1, j) + S(i, j+1) + S(i, j-1);
if(S(i,j) == 0 && val == 0)
{
U(i, j) = 0;
V(i, j) = 0;
continue;
}
// else if(S(i,j) == 2.0)
// {
// U(i, j) = -U(i, j);
// V(i, j) = -V(i, j);
// continue;
// }
U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
}

5
BasicSolver/2D-seq/src/solver.h
View File

@@ -9,6 +9,8 @@
#include "parameter.h"
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
/// @brief
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
typedef struct {
/* geometry and grid information */
@@ -18,7 +20,7 @@ typedef struct {
/* arrays */
double *p, *rhs;
double *f, *g;
double *u, *v;
double *u, *v, *s;
/* parameters */
double eps, omega, rho;
double re, tau, gamma;
@@ -43,4 +45,5 @@ extern void setSpecialBoundaryCondition(Solver*);
extern void computeFG(Solver*);
extern void adaptUV(Solver*);
extern void writeResult(Solver*);
extern void print(Solver*, double*);
#endif

2
BasicSolver/2D-seq/src/util.h
View File

@@ -4,6 +4,8 @@
* Use of this source code is governed by a MIT-style
* license that can be found in the LICENSE file.
*/
#include <math.h>
#ifndef __UTIL_H_
#define __UTIL_H_
#define HLINE \