Added Euler equation

parent 51d6a80e
......@@ -6,7 +6,7 @@
Где
* **-w** - ожидание резльтатов(может быть долго, если много потоков, и лучше его не выставлять)
* **-t matrix/thermal/oscillation/elliptic** - задание которое запускать
* **-t matrix/thermal/oscillation/elliptic/euler** - задание которое запускать
* **-np число_потоков**
* **-o вариант** - вариант по методичке(нумерация от 1). Для матрицы не актуально
* Варинаты аргументов для разных задач:
......@@ -14,6 +14,7 @@
* **-N число\_частей\_по\_оси\_X -t\_therm временной_отрезок** - для уравнения теплопроводности
* **-t\_osc временной\_отрезок -dt\_osc шаг_по_времени -dx\_osc шаг\_по\_X** - для уравнения колебания струны
* **-dx\_osc шаг_по_X -dy\_osc шаг\_по\_Y** - для эллиптического уравнения
* **-dx\_eul шаг_по_X -dy\_eul шаг\_по\_Y -dt\_eul шаг\_по\_T -N\_eul число_шагов\_по\_X -M\_eul число_шагов\_по\_Y -M\_eul число_шагов\_по\_T** - для уравнения Эйлера
### Справка скрипта ###
```
......
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mpi.h"
const int MASTER_RANK = 0;
const int MASTER_SEND_DATA = 1;
const int MASTER_RECEIVE_DATA = 2;
const int EXCHANGE_GHOST_CELLS = 3;
struct Cell {
double u;
double v;
double E;
double ro;
double p;
};
void initValues(Cell *cells, int rows, int cols, int option);
void boundConditions(Cell *cells, int rows, int cols, bool bounds[2], int option);
void calculatePressure(Cell *cells, int rows, int cols);
void eulerStep(
Cell *cells,
Cell *cellsE,
int rows, int cols, bool bounds[2],
double dt, double dx, double dy
);
void finalCalc(
Cell *cells,
Cell *cellsE,
int rows, int cols, bool bounds[2],
double dt, double dx, double dy
);
int main(int argc, char *argv[])
{
int rank;
int size;
double dt, dx, dy;
unsigned int N, M, T;
int option;
Cell *cells;
if(argc > 7) {
option = atoi(argv[1]);
dx = atof(argv[2]);
dy = atof(argv[3]);
dt = atof(argv[4]);
N = atoi(argv[5]);
M = atoi(argv[6]);
T = atoi(argv[7]);
}
else {
exit(1);
}
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
int thCount = size - 1;
unsigned int rowsPerTh = M / thCount;
unsigned int extraRowCount = M % thCount;
MPI_Datatype cellDT;
MPI_Type_contiguous(5, MPI_DOUBLE, &cellDT);
MPI_Type_commit(&cellDT);
//Master thread aalocates data and send to worker threads
if (rank == MASTER_RANK) {
double start;
double end;
printf("dt: %lf\n", dt);
printf("dx: %lf\n", dx);
printf("dy: %lf\n", dy);
printf("X part num: %u\n", N);
printf("Y part num: %u\n", M);
printf("T part num: %u\n", T);
printf("Thread count: %d\n", thCount);
printf("Y parts per thread: %u\n", rowsPerTh);
printf("Extra Y parts: %u\n", extraRowCount);
cells = new Cell[N * M];
start = MPI_Wtime();
int offset = 0;
int rows;
Cell *part = new Cell[(N + 2) * (rowsPerTh + 1)];
for (int th = 1; th <= thCount; ++th)
{
MPI_Recv(&rows, 1, MPI_INT, th, MASTER_RECEIVE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(part, rows * (N + 2), cellDT, th, MASTER_RECEIVE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
for(int i = 0; i < rows; ++i)
{
for(int j = 0; j < N; ++j)
cells[(i + offset) * N + j] = part[i * (N + 2) + j + 1];
}
offset += rows;
}
end = MPI_Wtime();
printf("Work time: %lf\n", end - start);
}
else {
int rows = (rank <= extraRowCount) ? rowsPerTh + 1 : rowsPerTh;
int cols = N;
Cell *cellsE;
int sRows = rows + 4;
int sCols = cols + 2;
//add ghost boundary size of 2x1
cells = new Cell[sRows * sCols];
cellsE = new Cell[sRows * sCols];
//Bound values for exchange in threds. Contains u[0/1], E, ro
Cell *bound = new Cell[2 * cols];
bool bounds[2];
bounds[0] = (rank == 1);
bounds[1] = (rank == thCount);
initValues(cells, rows, cols, option);
for(int t = 0; t < T; ++t)
{
//Calculate
boundConditions(cells, sRows, sCols, bounds, option);
calculatePressure(cells, sRows, sCols);
eulerStep(cells, cellsE, sRows, sCols, bounds, dx, dy, dt);
boundConditions(cellsE, sRows, sCols, bounds, option);
finalCalc(cells, cellsE, sRows, sCols, bounds, dx, dy, dt);
if (rank % 2 == 1)
{
if(rank < size - 1)//Send 2 bottom ghost bounds
{
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
bound[i * cols + j] = cells[(i + rows) * sCols + j + 1];
}
MPI_Send(bound, 2 * cols, cellDT, rank + 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD);
MPI_Recv(bound, 2 * cols, cellDT, rank + 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
cells[(i + rows + 2) * sCols + j + 1] = bound[i * cols + j];
}
}
if(rank > 1)
{
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
bound[i * cols + j] = cells[(i + 2) * sCols + j + 1];
}
MPI_Send(bound, 2 * cols, cellDT, rank - 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD);
MPI_Recv(bound, 2 * cols, cellDT, rank - 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
cells[i * sCols + j + 1] = bound[i * cols + j];
}
}
}
else
{
if(rank > 1)
{
MPI_Recv(bound, 2 * cols, cellDT, rank - 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
cells[i * sCols + j + 1] = bound[i * cols + j];
}
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
bound[i * cols + j] = cells[(i + 2) * sCols + j + 1];
}
MPI_Send(bound, 2 * cols, cellDT, rank - 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD);
}
if(rank < size - 1)//Send 2 bottom ghost bounds
{
MPI_Recv(bound, 2 * cols, cellDT, rank + 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
cells[(i + rows + 2) * sCols + j + 1] = bound[i * cols + j];
}
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < cols; ++j)
bound[i * cols + j] = cells[(i + rows) * sCols + j + 1];
}
MPI_Send(bound, 2 * cols, cellDT, rank + 1, EXCHANGE_GHOST_CELLS, MPI_COMM_WORLD);
}
}
}
MPI_Send(&rows, 1, MPI_INT, MASTER_RANK, MASTER_RECEIVE_DATA, MPI_COMM_WORLD);
MPI_Send(cells + 2 * sCols, rows * sCols, cellDT, MASTER_RANK, MASTER_RECEIVE_DATA, MPI_COMM_WORLD);
}
MPI_Finalize();
}
const double Max = 2.0;
const double g = 1.4;
void initValues(Cell *cells, int rows, int cols, int option)
{
double e0 = .5 + 1./(g * (g - 1) * Max * Max);
for(int i = 0; i < rows; ++i)
{
for(int j = 0; j < cols; ++j)
{
cells[i * cols + j].u = 1.;
cells[i * cols + j].v = .1e-15;
cells[i * cols + j].ro = 1.;
cells[i * cols + j].E = e0;
cells[i * cols + j].p = cells[i * cols + j].ro * (g-1) * (e0-.5*(
cells[i * cols + j].u * cells[i * cols + j].u +
cells[i * cols + j].v * cells[i * cols + j].v
));
}
}
}
void boundConditions(Cell *cells, int rows, int cols, bool bounds[2], int option)
{
double e0 = .5 + 1./(g * (g - 1) * Max * Max);
int sRow = 2;
int eRow = rows - 2;
int sCol = 1;
int eCol = cols - 1;
for(int i = sRow; i < eRow; ++i)
{
cells[i * cols].u = 1.;
cells[i * cols].v = .1e-15;
cells[i * cols].E = e0;
cells[i * cols].ro = 1.;
cells[i * cols + eCol].u = 2.*cells[i * cols + eCol - 1].u - cells[i * cols + eCol - 2].u;
cells[i * cols + eCol].v = 2.*cells[i * cols + eCol - 1].v - cells[i * cols + eCol - 2].v;
cells[i * cols + eCol].E = 2.*cells[i * cols + eCol - 1].E - cells[i * cols + eCol - 2].E;
cells[i * cols + eCol].ro = 2.*cells[i * cols + eCol - 1].ro - cells[i * cols + eCol - 2].ro;
}
for(int j = sCol; j < eCol; ++j)
{
if(bounds[0])
{
cells[1 * cols + j].u = cells[2 * cols + j].u;
cells[1 * cols + j].v = -cells[2 * cols + j].v;
cells[1 * cols + j].E = cells[2 * cols + j].E;
cells[1 * cols + j].ro = cells[2 * cols + j].ro;
}
if(bounds[1])
{
cells[eRow * cols + j].u = cells[(eRow - 1) * cols + j].u;
cells[eRow * cols + j].v = cells[(eRow - 1) * cols + j].v;
cells[eRow * cols + j].E = cells[(eRow - 1) * cols + j].E;
cells[eRow * cols + j].ro = cells[(eRow - 1) * cols + j].ro;
}
}
}
void calculatePressure(Cell *cells, int rows, int cols)
{
for(int i = 0; i < rows; ++i)
{
for(int j = 0; j < cols; ++j)
{
cells[i * cols + j].p = cells[i * cols + j].ro * (g-1)*(cells[i * cols + j].E-.5*(
cells[i * cols + j].u * cells[i * cols + j].u +
cells[i * cols + j].v * cells[i * cols + j].v
));
}
}
}
void eulerStep(
Cell *cells,
Cell *cellsE,
int rows, int cols, bool bounds[2],
double dt, double dx, double dy
){
int sRow = bounds[0] ? 2 : 1;
int eRow = bounds[1] ? rows - 2 : rows - 1;
int sCol = 1;
int eCol = cols - 1;
for(int i = sRow; i < eRow; ++i)
{
for(int j = sCol; j < eCol; ++j)
{
double p0 = cells[i * cols + j].p;
double pl = .5 * (p0 + cells[i * cols + j - 1].p);
double pr = .5 * (p0 + cells[i * cols + j + 1].p);
double pb = .5 * (p0 + cells[(i - 1) * cols + j].p);
double pt = .5 * (p0 + cells[(i + 1) * cols + j].p);
double u01 = cells[i * cols + j].u;
double u02 = cells[i * cols + j].v;
double ul1 = .5 * (u01 + cells[i * cols + j - 1].u);
double ur1 = .5 * (u01 + cells[i * cols + j + 1].u);
double ub2 = .5 * (u02 + cells[(i - 1) * cols + j].v);
double ut2 = .5 * (u02 + cells[(i + 1) * cols + j].v);
cellsE[i * cols + j].u =
cells[i * cols + j].u +
(pl - pr) * dt / ( dx * cells[i * cols + j].ro);
cellsE[i * cols + j].v =
cells[i * cols + j].v +
(pb - pt) * dt / (dy * cells[i * cols + j].ro);
cellsE[i * cols + j].E =
cells[i * cols + j].E +
(pl * ul1 - pr * ur1) * dt / (dx * cells[i * cols + j].ro) +
(pb * ub2 - pt * ut2) * dt / (dy * cells[i * cols + j].ro);
}
}
}
void finalCalc(
Cell *cells,
Cell *cellsE,
int rows, int cols, bool bounds[2],
double dt, double dx, double dy
){
int sRow = 2;
int eRow = rows - 2;
int sCol = 1;
int eCol = cols - 1;
double dm1, d1, dm2, d2, dm3, d3, dm4, d4;
for(int i = sRow; i < eRow; ++i)
{
for(int j = sCol; j < eCol; ++j)
{
double ro0 = cells[i * cols + j].ro;
double u10 = cellsE[i * cols + j].u;
double u20 = cellsE[i * cols + j].v;
double ul = .5 * (u10 + cellsE[i * cols + j - 1].u);
double ur = .5 * (u10 + cellsE[i * cols + j + 1].u);
double ub = .5 * (u20 + cellsE[(i - 1) * cols + j].v);
double ut = .5 * (u20 + cellsE[(i + 1) * cols + j].v);
if(ul <= 0.)
{
dm1= ul * ro0 * dy * dt;
d1 = 0.;
}
else
{
dm1= ul * cells[i * cols + j - 1].ro * dy * dt;
d1 = 1.;
}
if(ub <= 0.)
{
dm4 = ub * ro0 * dx * dt;
d4 = 0.;
}
else
{
dm4= ub * cells[(i - 1) * cols + j].ro * dx * dt;
d4 = 1.;
}
if(ur <= 0.)
{
dm3= ur * cells[i * cols + j + 1].ro * dy * dt;
d3 = 1.;
}
else
{
dm3= ur * ro0 * dy * dt;
d3 = 0.;
}
if(ut > 0.)
{
dm2 = ut * cells[(i + 1) * cols + j].ro * dx * dt;
d2 = 1.;
}
else
{
dm4= ut * ro0 * dx * dt;
d4 = 0.;
}
double am1 = fabs(dm1);
double am2 = fabs(dm2);
double am3 = fabs(dm3);
double am4 = fabs(dm4);
double z1 = dx * dy;
double z2 = (d1-.5) * dm1 + (d2-.5) * dm2 + (d3-.5) * dm3 + (d4-.5) * dm4;
cellsE[i * cols + j].ro = ro0 + 2. * z2 / z1;
double z3 = ro0 * z1 - (1.-d1) * dm1 - (1.-d2) * dm2 - (1.-d3) * dm3 - (1.-d4) * dm4;
double z4 = cellsE[i * cols + j].ro * z1;
double b1 = d1 * dm1;
double b2 = d2 * dm2;
double b3 = d3 * dm3;
double b4 = d4 * dm4;
cells[i * cols + j].u = (
u10 * z3 +
cellsE[i * cols + j - 1].u * b1 +
cellsE[(i - 1) * cols + j].u * b2 +
cellsE[i * cols + j + 1].u * b3 +
cellsE[(i + 1) * cols + j].u * b4
) / z4;
cells[i * cols + j].v = (
u20 * z3 +
cellsE[i * cols + j - 1].v * b1 +
cellsE[(i - 1) * cols + j].v * b2 +
cellsE[i * cols + j + 1].v * b3 +
cellsE[(i + 1) * cols + j].v * b4
) / z4;
cells[i * cols + j].E = (
cellsE[i * cols + j].E * z3 +
cellsE[i * cols + j - 1].E * b1 +
cellsE[(i - 1) * cols + j].E * b2 +
cellsE[i * cols + j + 1].E * b3 +
cellsE[(i + 1) * cols + j].E * b4
) / z4;
cells[i * cols + j].ro = cellsE[i * cols + j].ro;
}
}
}
......@@ -11,7 +11,7 @@ def _parse_arguments(argv):
parser = ArgumentParser(prog="run.py",
description='\nCompiles and run MPI tasks')
parser.add_argument("-t", "--task",dest="task", choices=['matrix','thermal','oscillation', 'elliptic'], default='matrix',
parser.add_argument("-t", "--task",dest="task", choices=['matrix','thermal','oscillation', 'elliptic', 'euler'], default='matrix',
help='Task to run', required=True)
parser.add_argument("-w", "--wait",action="store_true", dest="wait",
......@@ -44,6 +44,20 @@ def _parse_arguments(argv):
group.add_argument("-dy_ell", dest="dy", type=float, help='Specify the dy')
group = parser.add_argument_group('Euler equation options')
group.add_argument("-dx_eul", dest="dx", type=float, help='Specify the dx')
group.add_argument("-dy_eul", dest="dy", type=float, help='Specify the dy')
group.add_argument("-dt_eul", dest="dt", type=float, help='Specify the dy')
group.add_argument("-N_eul", dest="N", type=int, help='Specify the N')
group.add_argument("-M_eul", dest="M", type=int, help='Specify the M')
group.add_argument("-T_eul", dest="T", type=int, help='Specify the T')
return parser.parse_known_args(argv)
if __name__ == "__main__":
......@@ -74,6 +88,8 @@ if __name__ == "__main__":
run_args = '%d %lf %lf %lf' % (args.option, args.T, args.dx, args.dt)
elif args.task == 'elliptic':
run_args = '%d %lf %lf' % (args.option, args.dx, args.dt)
elif args.task == 'euler':
run_args = '%d %lf %lf %lf %d %d %d' % (args.option, args.dx, args.dy, args.dt, args.N, args.M, args.T)
if not args.wait:
os.system('mpirun -np %d %s %s' % (args.threads + 1, run_file, run_args))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment