WIP: Pull Request for a complete Solver package #1

Closed
AdityaUjeniya wants to merge 51 commits from (deleted):main into main
12 changed files with 18951 additions and 20052 deletions
Showing only changes of commit e9e69c2e07 - Show all commits

View File

@ -15,24 +15,24 @@ bcRight 3 #
gx 0.0 # Body forces (e.g. gravity) gx 0.0 # Body forces (e.g. gravity)
gy 0.0 # gy 0.0 #
re 30000.0 # Reynolds number re 16500.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction v_init 0.0 # initial value for velocity in y-direction
p_init 0.0 # initial value for pressure p_init 1.0 # initial value for pressure
# Geometry Data: # Geometry Data:
# ------------- # -------------
xlength 30.0 # domain size in x-direction xlength 7.0 # domain size in x-direction
ylength 4.0 # domain size in y-direction ylength 1.5 # domain size in y-direction
imax 200 # number of interior cells in x-direction imax 210 # number of interior cells in x-direction
jmax 50 # number of interior cells in y-direction jmax 45 # number of interior cells in y-direction
# Time Data: # Time Data:
# --------- # ---------
te 100.0 # final time te 30.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)
@ -50,13 +50,13 @@ gamma 0.9 # upwind differencing factor gamma
numberOfParticles 60 numberOfParticles 60
startTime 0 startTime 0
injectTimePeriod 2.0 injectTimePeriod 0.5
writeTimePeriod 0.5 writeTimePeriod 0.2
x1 1.0 x1 0.0
y1 0.0 y1 0.5
x2 1.0 x2 0.0
y2 4.0 y2 1.5
# Obstacle Geometry Data: # Obstacle Geometry Data:
# ----------------------- # -----------------------
@ -65,7 +65,7 @@ y2 4.0
shape 1 shape 1
xCenter 0.0 xCenter 0.0
yCenter 0.0 yCenter 0.0
xRectLength 8.0 xRectLength 2.0
yRectLength 1.0 yRectLength 1.0
circleRadius 1.0 circleRadius 1.0

View File

@ -15,7 +15,7 @@ bcRight 3 #
gx 0.0 # Body forces (e.g. gravity) gx 0.0 # Body forces (e.g. gravity)
gy 0.0 # gy 0.0 #
re 500.0 # Reynolds number re 100.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction v_init 0.0 # initial value for velocity in y-direction
@ -62,11 +62,11 @@ y2 4.0
# ----------------------- # -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle # Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 0 shape 1
xCenter 4.0 xCenter 10.0
yCenter 1 yCenter 2
xRectLength 8.0 xRectLength 6.0
yRectLength 2.0 yRectLength 1.0
circleRadius 1.0 circleRadius 1.0
#=============================================================================== #===============================================================================

View File

@ -63,7 +63,7 @@ y2 0.9
# ----------------------- # -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle # Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 1 shape 0
xCenter 0.5 xCenter 0.5
yCenter 0.5 yCenter 0.5
xRectLength 0.5 xRectLength 0.5

Binary file not shown.

After

Width:  |  Height:  |  Size: 24 KiB

File diff suppressed because it is too large Load Diff

View File

@ -134,6 +134,7 @@ void initSolver(Solver* solver, Parameter* params)
solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum; solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum;
int iCenter = 0, jCenter = 0, iRectLength = 0, jRectLength = 0, radius = 0; int iCenter = 0, jCenter = 0, iRectLength = 0, jRectLength = 0, radius = 0;
double x1=0, x2=0, y1=0, y2=0;
int* s = solver->s; int* s = solver->s;
@ -142,25 +143,16 @@ void initSolver(Solver* solver, Parameter* params)
case NOSHAPE: case NOSHAPE:
break; break;
case RECT: case RECT:
iCenter = params->xCenter/dx; x1 = params->xCenter - params->xRectLength/2;
jCenter = params->yCenter/dy; x2 = params->xCenter + params->xRectLength/2;
iRectLength = params->xRectLength/dx; y1 = params->yCenter - params->yRectLength/2;
jRectLength = params->yRectLength/dy; y2 = params->yCenter + params->yRectLength/2;
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 j = 1; j < jmax + 1; j++)
{ {
for (int i = 1; i < imax + 1; i++) for (int i = 1; i < imax + 1; i++)
{ {
if(i1 <= i && i <= i2 && j1 <= j && j <= j2) if((x1 <= (i*dx)) && ((i*dx) <= x2) && (y1 <= (j*dy)) && ((j*dy) <= y2))
{ {
S(i, j) = LOCAL; S(i, j) = LOCAL;
} }
@ -171,7 +163,7 @@ void initSolver(Solver* solver, Parameter* params)
case CIRCLE: case CIRCLE:
iCenter = params->xCenter/dx; iCenter = params->xCenter/dx;
jCenter = params->yCenter/dy; jCenter = params->yCenter/dy;
radius = params->circleRadius/dx; // have to change the logic to get max out of both the dx and dy radius = params->circleRadius;
for (int j = 1; j < jmax + 1; j++) for (int j = 1; j < jmax + 1; j++)
{ {
@ -590,6 +582,7 @@ void setSpecialBoundaryCondition(Solver* solver)
int jmax = solver->jmax; int jmax = solver->jmax;
double mDy = solver->dy; double mDy = solver->dy;
double* u = solver->u; double* u = solver->u;
int* s = solver->s;
if (strcmp(solver->problem, "dcavity") == 0) { if (strcmp(solver->problem, "dcavity") == 0) {
for (int i = 1; i < imax; i++) { for (int i = 1; i < imax; i++) {
@ -605,13 +598,11 @@ void setSpecialBoundaryCondition(Solver* solver)
U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
} }
} }
else if (strcmp(solver->problem, "backstep") == 0) { else if (strcmp(solver->problem, "backstep") == 0)
double ylength = solver->ylength; {
double y; for (int j = 1; j < jmax + 1; j++)
{
for (int j = 1; j < jmax + 1; j++) { if(S(0,j) == NONE) U(0, j) = 1.0;
//y = mDy * (j - 0.5);
U(0, j) = 0.0;
} }
} }
} }

View File

@ -1,7 +1,7 @@
set terminal png size 1800,768 enhanced font ,12 set terminal png size 1800,768 enhanced font ,12
set output 'velocity.png' set output 'velocity.png'
set datafile separator whitespace set datafile separator whitespace
#set object 1 rect from 8.5,1.5 to 11.5,2.5 lw 5 set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 5
plot 'velocity.dat' using 1:2:3:4:5 with vectors filled head size 0.01,20,60 lc palette plot 'velocity.dat' using 1:2:3:4:5 with vectors filled head size 0.01,20,60 lc palette

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 126 KiB

View File

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

View File

@ -0,0 +1,12 @@
unset border; unset tics; unset key;
set term gif animate delay 50
set output "trace.gif"
set xrange [0:7]
set yrange [0:1.5]
set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 5
do for [ts=0:300] {
plot "particles_".ts.".dat" with points pointtype 7
}
unset output

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.8 MiB