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

Working basic PT code

Browse Source
This commit is contained in:
Aditya Ujeniya
2023-08-06 21:22:54 +02:00
parent 7b83a34e47
commit 6ebcc5532e
96 changed files with 45514 additions and 20085 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -57,7 +57,7 @@ int main (int argc, char** argv)
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
S = getTimeStamp();
while (t <= te)
{
if (tau > 0.0) computeTimestep(&solver);
@@ -70,10 +70,11 @@ 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, solver.u, solver.v, t);
t += solver.dt;
nt++;
/*
#ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
@@ -82,7 +83,7 @@ int main (int argc, char** argv)
#endif
*/
}
printf("Total particles : %d\n", particletracer.totalParticles);
E = getTimeStamp();
stopProgress();

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

@@ -16,39 +16,63 @@
#include "particletracing.h"
#include "util.h"
#define U(i, j) u[(j) * (imax + 2) + (i)]
#define V(i, j) v[(j) * (imax + 2) + (i)]
#define U(i, j) u[(j) * (imax + 2) + (i)]
#define V(i, j) v[(j) * (imax + 2) + (i)]
static void print(Solver* solver, double* grid)
{
int imax = solver->imax;
for (int j = 0; j < solver->jmax + 2; j++) {
printf("%02d: ", j);
for (int i = 0; i < solver->imax + 2; i++) {
printf("%3.2f ", grid[j * (imax + 2) + i]);
}
printf("\n");
}
fflush(stdout);
}
static int ts = 0;
void printParticles(ParticleTracer* particletracer)
{
for(int i = 0; i < particletracer->totalParticles; ++i)
{
printf("Particle position X : %.2f, Y : %.2f, flag : %d\n", particletracer->particlePool[i].x,
particletracer->particlePool[i].y,
particletracer->particlePool[i].flag);
}
}
void injectParticles(ParticleTracer* particletracer)
{
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].flag = particletracer->linSpaceLine[i].flag;
particletracer->particlePool[particletracer->pointer].flag = true;
++(particletracer->pointer);
++(particletracer->totalParticles);
}
}
void advanceParticles(ParticleTracer* particletracer, double* u, double* v, double time)
void advanceParticles(ParticleTracer* particletracer, Solver* solver, double* restrict u, double* restrict v, double time)
{
Particle* particlePool = particletracer->particlePool;
int imax = particletracer->imax;
int jmax = particletracer->jmax;
double dx = particletracer->dx;
double dy = particletracer->dy;
double xlength = particletracer->xlength;
double ylength = particletracer->ylength;
for(int i = 0; i < particletracer->totalParticles; ++i)
{
if(particlePool[i].flag == true)
if(particletracer->particlePool[i].flag == true)
{
double x = particlePool[i].x;
double y = particlePool[i].y;
double x = particletracer->particlePool[i].x;
double y = particletracer->particlePool[i].y;
int iCoord = (int)(x / dx) + 1;
int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
@@ -65,7 +89,7 @@ void advanceParticles(ParticleTracer* particletracer, double* u, double* v, doub
(x - x1) * (y - y1) * U(iCoord, jCoord));
double new_x = x + particletracer->dt * u_n;
particlePool[i].x = new_x;
particletracer->particlePool[i].x = new_x;
iCoord = (int)((x + 0.5 * dx) / dx) + 1;
jCoord = (int)(y / dy) + 1;
@@ -82,24 +106,18 @@ void advanceParticles(ParticleTracer* particletracer, double* u, double* v, doub
(x - x1) * (y - y1) * V(iCoord, jCoord));
double new_y = y + particletracer->dt * v_n;
particlePool[i].y = new_y;
}
}
double xlength = particletracer->xlength;
double ylength = particletracer->ylength;
particletracer->particlePool[i].y = new_y;
for(int i = 0; i < particletracer->totalParticles; ++i)
{
if(particlePool[i].flag == true)
{
double x = particlePool[i].x;
double y = particlePool[i].y;
// printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y : %.2f, jCoord : %d\n\n", x, new_x, iCoord, y, new_y, jCoord);
// printf("\tU(iCoord - 1, jCoord - 1) : %.2f, U(iCoord, jCoord - 1) : %.2f, U(iCoord - 1, jCoord) : %.2f, U(iCoord, jCoord) : %.2f\n", U(iCoord - 1, jCoord - 1), U(iCoord, jCoord - 1), U(iCoord - 1, jCoord), U(iCoord, jCoord));
// printf("\tV(iCoord - 1, jCoord - 1) : %.2f, V(iCoord, jCoord - 1) : %.2f, V(iCoord - 1, jCoord) : %.2f, V(iCoord, jCoord) : %.2f\n\n", V(iCoord - 1, jCoord - 1), V(iCoord, jCoord - 1), V(iCoord - 1, jCoord), V(iCoord, jCoord));
// printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n);
// printf("\t j-1 * (imax + 2) + i-1 = %d with element from U : %.2f", (jCoord - 1) * (200 + 2) + (iCoord - 1), u[(jCoord - 1) * (imax + 2) + (iCoord - 1)]);
// printf("\nimax : %d, jmax : %d\n", imax, jmax);
if (((x < 0.0) || (x > xlength) || (y < 0.0) || (y > ylength)))
if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) || (new_y > ylength)))
{
particletracer->totalParticles -= 1;
particlePool[i].flag = false;
particletracer->particlePool[i].flag = false;
}
}
}
@@ -144,7 +162,13 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
particletracer->injectTimePeriod = params->injectTimePeriod;
particletracer->writeTimePeriod = params->writeTimePeriod;
particletracer->dt = params->dt;
particletracer->dx = params->xlength / params->imax;
particletracer->dy = params->ylength / params->jmax;
particletracer->xlength = params->xlength;
particletracer->ylength = params->ylength;
particletracer->x1 = params->x1;
particletracer->y1 = params->y1;
particletracer->x2 = params->x2;
@@ -157,19 +181,19 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
particletracer->pointer = 0;
particletracer->totalParticles = 0;
particletracer->imax = params->imax;
particletracer->jmax = params->jmax;
particletracer->estimatedNumParticles = ((params->te - params->startTime) + 2 ) * params->numberOfParticles;
particletracer->particlePool = (Particle*)malloc(sizeof(Particle) * particletracer->estimatedNumParticles);
particletracer->linSpaceLine = (Particle*)malloc(sizeof(Particle) * particletracer->numberOfParticles);
particletracer->particlePool = malloc(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);
double x = spacing * particletracer->x1 + (1.0 - spacing) * particletracer->x2;
double y = spacing * particletracer->y1 + (1.0 - spacing) * particletracer->y2;
particletracer->linSpaceLine[i].x = x;
particletracer->linSpaceLine[i].y = y;
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].flag = true;
}
@@ -181,12 +205,15 @@ void printParticleTracerParameters(ParticleTracer* particletracer)
printf("\tNumber of particles : %d being injected for every period of %.2f\n", particletracer->numberOfParticles, particletracer->injectTimePeriod);
printf("\tstartTime : %.2f\n", particletracer->startTime);
printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : %.2f, x2 : %.2f, y2 : %.2f\n", particletracer->x1, particletracer->y1, particletracer->x2, particletracer->y2);
}
printf("\tPointer : %d, TotalParticles : %d\n", particletracer->pointer, particletracer->totalParticles);
printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n", particletracer->dt, particletracer->dx, particletracer->dy);
}
void trace(ParticleTracer* particletracer, double* u, double* v, double time)
void trace(ParticleTracer* particletracer, Solver* solver, double* u, double* v, double time)
{
if (time >= particletracer->startTime)
{
//printParticles(particletracer);
if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
{
injectParticles(particletracer);
@@ -196,8 +223,9 @@ void trace(ParticleTracer* particletracer, double* u, double* v, double time)
{
writeParticles(particletracer);
particletracer->lastWriteTime = time;
}
advanceParticles(particletracer, u, v, time);
advanceParticles(particletracer, solver, u, v, time);
particletracer->lastUpdateTime = time;
}
}

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

@@ -35,9 +35,10 @@ typedef struct {
void initParticleTracer(ParticleTracer*, Parameter*);
void injectParticles(ParticleTracer*);
void advanceParticles(ParticleTracer*, double* , double*, double);
void advanceParticles(ParticleTracer*, Solver*, double* , double*, double);
void freeParticles(ParticleTracer*);
void writeParticles(ParticleTracer*);
void printParticleTracerParameters(ParticleTracer*);
void trace(ParticleTracer*, double* , double* , double );
void printParticles(ParticleTracer*);
void trace(ParticleTracer*, Solver*, double* , double* , double );
#endif

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

@@ -29,7 +29,7 @@ static void print(Solver* solver, double* grid)
for (int j = 0; j < solver->jmax + 2; j++) {
printf("%02d: ", j);
for (int i = 0; i < solver->imax + 2; i++) {
printf("%12.8f ", grid[j * (imax + 2) + i]);
printf("%3.2f ", grid[j * (imax + 2) + i]);
}
printf("\n");
}
@@ -150,6 +150,9 @@ void solve(Solver* solver)
int it = 0;
double res = 1.0;
// printf("\nBefore P : \n");
// print(solver, p);
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
@@ -176,15 +179,20 @@ void solve(Solver* solver)
}
res = res / (double)(imax * jmax);
/*
#ifdef DEBUG
printf("%d Residuum: %e\n", it, res);
#endif
*/
it++;
}
// printf("\nAfter P : \n");
// print(solver, p);
/*
#ifdef VERBOSE
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
#endif
*/
}
void solveRB(Solver* solver)
@@ -205,6 +213,9 @@ void solveRB(Solver* solver)
double res = 1.0;
int pass, jsw, isw;
// printf("\nBefore P : \n");
// print(solver, p);
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
jsw = 1;
@@ -238,11 +249,15 @@ void solveRB(Solver* solver)
}
res = res / (double)(imax * jmax);
/*
#ifdef DEBUG
printf("%d Residuum: %e\n", it, res);
#endif
*/
it++;
}
// printf("\nAfter P : \n");
// print(solver, p);
/*
#ifdef VERBOSE
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
@@ -464,6 +479,12 @@ void setBoundaryConditions(Solver* solver)
case PERIODIC:
break;
}
// printf("\nBoundary condition U : \n");
// print(solver, u);
// printf("\nBoundary condition V : \n");
// print(solver, v);
}
void setSpecialBoundaryCondition(Solver* solver)
@@ -486,6 +507,10 @@ void setSpecialBoundaryCondition(Solver* solver)
U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
}
}
// printf("\nSpecial Boundary condition U : \n");
// print(solver, u);
}
void computeFG(Solver* solver)
@@ -559,6 +584,13 @@ void computeFG(Solver* solver)
G(i, 0) = V(i, 0);
G(i, jmax) = V(i, jmax);
}
// printf("\nComputed F : \n");
// print(solver, f);
// printf("\nComputed G : \n");
// print(solver, g);
}
void adaptUV(Solver* solver)
@@ -579,6 +611,17 @@ void adaptUV(Solver* solver)
V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
}
}
// printf("\nApadted U : \n");
// print(solver, u);
// printf("\nAdapted V : \n");
// print(solver, v);
// char *b;
// b='\0';
// printf("\nDEBUGGER : ");
// scanf("%s", &b);
}
void writeResult(Solver* solver)