Particle Tracer work in progress

This commit is contained in:
Aditya Ujeniya 2023-08-05 19:35:16 +02:00
parent b50dbe7d4b
commit 7b83a34e47
46 changed files with 60399 additions and 23269 deletions

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

Binary file not shown.

Before

Width:  |  Height:  |  Size: 206 KiB

After

Width:  |  Height:  |  Size: 19 KiB

View File

@ -47,13 +47,13 @@ gamma 0.9 # upwind differencing factor gamma
# Visualization Data: # Visualization Data:
# ------------------ # ------------------
traceStart 10.0 # time for starting visualization traceStart 0.0 # time for starting visualization
traceWrite 2.0 # time stepsize for saving particle data traceWrite 2.0 # time stepsize for saving particle data
traceInject 2.0 # time stepsize for particle injection traceInject 2.0 # time stepsize for particle injection
lineX1 2.0 # Coordinates of line segment for particle injection lineX1 1.0 # Coordinates of line segment for particle injection
lineY1 0.5 lineY1 0.0
lineX2 2.0 lineX2 1.0
lineY2 3.5 lineY2 4.0
nparticles 30 # number of particles to inject nparticles 50 # number of particles to inject
#=============================================================================== #===============================================================================

View File

@ -1,5 +1,5 @@
# Supported: GCC, CLANG, ICC # Supported: GCC, CLANG, ICC
TAG ?= CLANG TAG ?= ICC
ENABLE_OPENMP ?= false ENABLE_OPENMP ?= false
#Feature options #Feature options

File diff suppressed because it is too large Load Diff

View File

@ -178,7 +178,7 @@ void solve(Solver* solver)
it++; it++;
} }
#ifdef VERBOSE #ifndef VERBOSE
printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
#endif #endif
} }

View File

@ -77,7 +77,6 @@ static void advanceParticles(
double yn = y + delt * vn; double yn = y + delt * vn;
m[particleId * NCOORD + Y] = yn; m[particleId * NCOORD + Y] = yn;
printf("P%i VEL %f %f dt %f OP %f %f NP %f %f\n", i, un, vn, delt, x, y, xn, yn);
} }
double xlength = t->grid.xlength; double xlength = t->grid.xlength;
@ -107,7 +106,6 @@ static void injectParticles(Tracing* t)
double* m = t->memorypool; double* m = t->memorypool;
for (int i = 0; i < t->numParticles; i++) { for (int i = 0; i < t->numParticles; i++) {
printf("Inject %d as %d mem %d\n", i, t->totalParticles, t->cursor);
t->particles[t->totalParticles] = t->cursor; t->particles[t->totalParticles] = t->cursor;
m[(t->cursor) * NCOORD + X] = line[i * NCOORD + X]; m[(t->cursor) * NCOORD + X] = line[i * NCOORD + X];
m[(t->cursor) * NCOORD + Y] = line[i * NCOORD + Y]; m[(t->cursor) * NCOORD + Y] = line[i * NCOORD + Y];
@ -145,14 +143,11 @@ void trace(Tracing* t, double* restrict u, double* restrict v, double time)
{ {
if (time >= t->traceStart) { if (time >= t->traceStart) {
if ((time - t->lastUpdate[INJECT]) > t->traceInject) { if ((time - t->lastUpdate[INJECT]) > t->traceInject) {
printf("Inject at %f\n", time);
printState(t);
injectParticles(t); injectParticles(t);
t->lastUpdate[INJECT] = time; t->lastUpdate[INJECT] = time;
} }
if ((time - t->lastUpdate[WRITE]) > t->traceWrite) { if ((time - t->lastUpdate[WRITE]) > t->traceWrite) {
printf("Write at %f\n", time);
writeParticles(t); writeParticles(t);
t->lastUpdate[WRITE] = time; t->lastUpdate[WRITE] = time;
} }
@ -199,7 +194,7 @@ void initTrace(Tracing* t, Parameter* p)
double x = spacing * x1 + (1.0 - spacing) * x2; double x = spacing * x1 + (1.0 - spacing) * x2;
double y = spacing * y1 + (1.0 - spacing) * y2; double y = spacing * y1 + (1.0 - spacing) * y2;
printf("S: %f x: %f y: %f\n", spacing, x, y); //printf("S: %f x: %f y: %f\n", spacing, x, y);
line[i * NCOORD + X] = x; line[i * NCOORD + X] = x;
line[i * NCOORD + Y] = y; line[i * NCOORD + Y] = y;
} }

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

View File

@ -0,0 +1,9 @@
unset border; unset tics; unset key;
set term gif animate delay 50
set output "trace.gif"
set xrange [0:30]
set yrange [0:4]
do for [ts=0:23] {
plot "particles_".ts.".dat" with points pointtype 7
}
unset output

View File

@ -32,7 +32,7 @@ jmax 50 # number of interior cells in y-direction
# Time Data: # Time Data:
# --------- # ---------
te 100.0 # final time te 60.0 # final time
dt 0.02 # time stepsize dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt) tau 0.5 # safety factor for time stepsize control (<0 constant delt)
@ -44,4 +44,18 @@ eps 0.00001 # stopping tolerance for pressure iteration
rho 0.52 rho 0.52
omg 1.8 # relaxation parameter for SOR iteration omg 1.8 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
# -----------------------
numberOfParticles 10
startTime 10
injectTimePeriod 2.0
writeTimePeriod 2.0
x1 5.0
y1 0.0
x2 5.0
y2 4.0
#=============================================================================== #===============================================================================

View File

@ -44,4 +44,16 @@ eps 0.001 # stopping tolerance for pressure iteration
rho 0.5 rho 0.5
omg 1.7 # relaxation parameter for SOR iteration omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma gamma 0.9 # upwind differencing factor gamma
# Particle Tracing Data:
# -----------------------
numberOfParticles 10
startTime 40.0
timePeriod 2.0
x1 5.0
y1 0.0
x2 5.0
y2 4.0
#=============================================================================== #===============================================================================

View File

@ -0,0 +1,10 @@
5.000000 4.000000
5.000000 3.555556
5.000000 3.111111
5.000000 2.666667
5.000000 2.222222
5.000000 1.777778
5.000000 1.333333
5.000000 0.888889
5.000000 0.444444
5.000000 0.000000

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

View File

File diff suppressed because it is too large Load Diff

View File

@ -13,6 +13,7 @@
#include "parameter.h" #include "parameter.h"
#include "progress.h" #include "progress.h"
#include "solver.h" #include "solver.h"
#include "particletracing.h"
#include "timing.h" #include "timing.h"
enum VARIANT { SOR = 1, RB, RBA }; enum VARIANT { SOR = 1, RB, RBA };
@ -25,6 +26,7 @@ int main (int argc, char** argv)
double S, E; double S, E;
Parameter params; Parameter params;
Solver solver; Solver solver;
ParticleTracer particletracer;
initParameter(&params); initParameter(&params);
if (argc < 2) { if (argc < 2) {
@ -40,6 +42,9 @@ int main (int argc, char** argv)
printParameter(&params); printParameter(&params);
initSolver(&solver, &params); initSolver(&solver, &params);
initParticleTracer(&particletracer, &params);
printParticleTracerParameters(&particletracer);
#ifndef VERBOSE #ifndef VERBOSE
initProgress(solver.te); initProgress(solver.te);
#endif #endif
@ -52,7 +57,7 @@ int main (int argc, char** argv)
void (*solver_generic[])(solver) = {solve, solveRB, solveRBA}; void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
S = getTimeStamp(); S = getTimeStamp();
while (t <= te) while (t <= te)
{ {
if (tau > 0.0) computeTimestep(&solver); if (tau > 0.0) computeTimestep(&solver);
@ -63,19 +68,26 @@ int main (int argc, char** argv)
if (nt % 100 == 0) normalizePressure(&solver); if (nt % 100 == 0) normalizePressure(&solver);
(*solver_generic[variant - 1])(&solver); (*solver_generic[variant - 1])(&solver);
adaptUV(&solver); adaptUV(&solver);
/* Added function for particle tracing. Will inject and advance particles as per timePeriod */
trace(&particletracer, solver.u, solver.v, t);
t += solver.dt; t += solver.dt;
nt++; nt++;
/*
#ifdef VERBOSE #ifdef VERBOSE
printf("TIME %f , TIMESTEP %f\n", t, solver.dt); printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
#else #else
printProgress(t); printProgress(t);
#endif #endif
*/
} }
E = getTimeStamp(); E = getTimeStamp();
stopProgress(); stopProgress();
freeParticles(&particletracer);
printf("Solution took %.2fs\n", E - S); printf("Solution took %.2fs\n", E - S);
writeResult(&solver); writeResult(&solver);
return EXIT_SUCCESS; return EXIT_SUCCESS;

View File

@ -81,6 +81,15 @@ void readParameter(Parameter* param, const char* filename)
PARSE_REAL(p_init); PARSE_REAL(p_init);
PARSE_REAL(rho); PARSE_REAL(rho);
/* Added new particle tracing parameters */
PARSE_INT(numberOfParticles);
PARSE_REAL(startTime);
PARSE_REAL(injectTimePeriod);
PARSE_REAL(writeTimePeriod);
PARSE_REAL(x1);
PARSE_REAL(y1);
PARSE_REAL(x2);
PARSE_REAL(y2);
} }
} }
@ -113,4 +122,11 @@ void printParameter(Parameter* param)
printf("\tomega (SOR relaxation): %f\n", param->omg); printf("\tomega (SOR relaxation): %f\n", param->omg);
printf("\trho (SOR relaxation): %f\n", param->rho); printf("\trho (SOR relaxation): %f\n", param->rho);
printf("Particle Tracing data:\n");
printf("\tNumber of particles : %d being injected for every period of %.2f\n", param->numberOfParticles, param->injectTimePeriod);
printf("\tstartTime : %.2f\n", param->startTime);
printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : %.2f, x2 : %.2f, y2 : %.2f\n", param->x1, param->y1, param->x2, param->y2);
} }

View File

@ -18,6 +18,11 @@ typedef struct {
char* name; char* name;
int bcLeft, bcRight, bcBottom, bcTop; int bcLeft, bcRight, bcBottom, bcTop;
double u_init, v_init, p_init; double u_init, v_init, p_init;
int numberOfParticles;
double startTime, injectTimePeriod, writeTimePeriod;
double x1, y1, x2, y2;
} Parameter; } Parameter;
void initParameter(Parameter*); void initParameter(Parameter*);

View File

@ -0,0 +1,203 @@
/*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file.
*/
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "allocate.h"
#include "parameter.h"
#include "solver.h"
#include "particletracing.h"
#include "util.h"
#define U(i, j) u[(j) * (imax + 2) + (i)]
#define V(i, j) v[(j) * (imax + 2) + (i)]
static int ts = 0;
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->pointer);
++(particletracer->totalParticles);
}
}
void advanceParticles(ParticleTracer* particletracer, double* u, double* v, double time)
{
Particle* particlePool = particletracer->particlePool;
int imax = particletracer->imax;
int jmax = particletracer->jmax;
double dx = particletracer->dx;
double dy = particletracer->dy;
for(int i = 0; i < particletracer->totalParticles; ++i)
{
if(particlePool[i].flag == true)
{
double x = particlePool[i].x;
double y = particlePool[i].y;
int iCoord = (int)(x / dx) + 1;
int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
double x1 = (double)(iCoord - 1) * dx;
double y1 = ((double)(jCoord - 1) - 0.5) * dy;
double x2 = (double)iCoord * dx;
double y2 = ((double)jCoord - 0.5) * dy;
double u_n = (1.0 / (dx * dy)) *
((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) +
(x - x1) * (y2 - y) * U(iCoord, jCoord - 1) +
(x2 - x) * (y - y1) * U(iCoord - 1, jCoord) +
(x - x1) * (y - y1) * U(iCoord, jCoord));
double new_x = x + particletracer->dt * u_n;
particlePool[i].x = new_x;
iCoord = (int)((x + 0.5 * dx) / dx) + 1;
jCoord = (int)(y / dy) + 1;
x1 = ((double)(iCoord - 1) - 0.5) * dx;
y1 = (double)(jCoord - 1) * dy;
x2 = ((double)iCoord - 0.5) * dx;
y2 = (double)jCoord * dy;
double v_n = (1.0 / (dx * dy)) *
((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) +
(x - x1) * (y2 - y) * V(iCoord, jCoord - 1) +
(x2 - x) * (y - y1) * V(iCoord - 1, jCoord) +
(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;
for(int i = 0; i < particletracer->totalParticles; ++i)
{
if(particlePool[i].flag == true)
{
double x = particlePool[i].x;
double y = particlePool[i].y;
if (((x < 0.0) || (x > xlength) || (y < 0.0) || (y > ylength)))
{
particletracer->totalParticles -= 1;
particlePool[i].flag = false;
}
}
}
}
void freeParticles(ParticleTracer* particletracer)
{
free(particletracer->particlePool);
free(particletracer->linSpaceLine);
}
void writeParticles(ParticleTracer* particletracer)
{
FILE* fp;
Particle* particlePool = particletracer->particlePool;
char filename[50];
snprintf(filename, 50, "particles_%d.dat", ts++);
fp = fopen(filename, "w");
if (fp == NULL) {
printf("Error!\n");
exit(EXIT_FAILURE);
}
for (int i = 0; i < particletracer->totalParticles; ++i)
{
if(particlePool[i].flag == true)
{
double x = particlePool[i].x;
double y = particlePool[i].y;
fprintf(fp, "%f %f\n", x, y);
}
}
fclose(fp);
}
void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
{
particletracer->numberOfParticles = params->numberOfParticles;
particletracer->startTime = params->startTime;
particletracer->injectTimePeriod = params->injectTimePeriod;
particletracer->writeTimePeriod = params->writeTimePeriod;
particletracer->x1 = params->x1;
particletracer->y1 = params->y1;
particletracer->x2 = params->x2;
particletracer->y2 = params->y2;
particletracer->lastInjectTime = params->startTime;
particletracer->lastUpdateTime = params->startTime;
particletracer->lastWriteTime = params->startTime;
particletracer->pointer = 0;
particletracer->totalParticles = 0;
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);
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;
particletracer->linSpaceLine[i].flag = true;
}
}
void printParticleTracerParameters(ParticleTracer* particletracer)
{
printf("Particle Tracing data:\n");
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);
}
void trace(ParticleTracer* particletracer, double* u, double* v, double time)
{
if (time >= particletracer->startTime)
{
if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
{
injectParticles(particletracer);
particletracer->lastInjectTime = time;
}
if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
{
writeParticles(particletracer);
particletracer->lastWriteTime = time;
}
advanceParticles(particletracer, u, v, time);
particletracer->lastUpdateTime = time;
}
}

View File

@ -0,0 +1,43 @@
/*
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
* All rights reserved. This file is part of nusif-solver.
* Use of this source code is governed by a MIT style
* license that can be found in the LICENSE file.
*/
#ifndef __PARTICLETRACING_H_
#define __PARTICLETRACING_H_
#include "parameter.h"
#include <stdbool.h>
typedef enum COORD { X = 0, Y, NCOORD } COORD;
typedef struct
{
double x, y;
bool flag;
} Particle;
typedef struct {
int numberOfParticles, totalParticles;
double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime, lastWriteTime;
int estimatedNumParticles;
double dx, dy, dt;
Particle* linSpaceLine;
Particle* particlePool;
int pointer;
double imax, jmax, xlength, ylength;
double x1, y1, x2, y2;
} ParticleTracer;
void initParticleTracer(ParticleTracer*, Parameter*);
void injectParticles(ParticleTracer*);
void advanceParticles(ParticleTracer*, double* , double*, double);
void freeParticles(ParticleTracer*);
void writeParticles(ParticleTracer*);
void printParticleTracerParameters(ParticleTracer*);
void trace(ParticleTracer*, double* , double* , double );
#endif

View File

@ -243,10 +243,11 @@ void solveRB(Solver* solver)
#endif #endif
it++; it++;
} }
/*
#ifdef VERBOSE #ifdef VERBOSE
printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
#endif #endif
*/
} }
void solveRBA(Solver* solver) void solveRBA(Solver* solver)

Binary file not shown.

After

Width:  |  Height:  |  Size: 12 KiB

File diff suppressed because it is too large Load Diff