计算流体力学，即Computational Fluid Dynamics（简称CFD），广泛应用于车辆设计、航空航天、宇宙天文、气象环境、核工业、国防军事等领域。其研究内容是通过计算机和数值方法来求解流体力学的控制方程，对流体力学问题进行模拟和分析。其中，Euler方法是在充满流体质点的空间内任取一控制微元所包含的流体作为研究对象，欧拉模型便是对每一相建立的单独守恒方程。
考试提供的Euler2D原始程序，是使用欧拉方法对二维区域进行流体模拟的串行程序,由C++语言编写。该程序二维区域中每个微元运动计算过程相互独立，因此具有良好的天然并行性，非常适合进行多核并行改造，从而大幅提升程序运行效率。
1、使用单线程运行
代码目录如下：
- cfd2d.cpp
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include "riemann.h"
#include "timer.h"
using namespace std;
#define gamma 2.1
float **Q1, **Q2, **Q3, **Q4, **dF1, **dF2, **dF3, **dF4, **dG1, **dG2, **dG3, **dG4;
float dx, dy, dt;
float **rho, **u, **v, **e, **p, **c;
int nmax, imax, jmax;
float cfl;
float **Array2D(unsigned int row, unsigned int col)
{
float **x;
x = new float *[row];
for (int count = 0; count < row; count++)
{
x[count] = new float[col];
}
return x;
}
void cons_2_prim()
{
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
rho[i][j] = Q1[i][j];
u[i][j] = Q2[i][j] / Q1[i][j];
v[i][j] = Q3[i][j] / Q1[i][j];
float ke = u[i][j] * u[i][j] / 2.0 + v[i][j] * v[i][j] / 2.0;
e[i][j] = Q4[i][j] / Q1[i][j] - ke;
p[i][j] = (gamma - 1.0) * rho[i][j] * e[i][j];
c[i][j] = sqrt(gamma * p[i][j] / rho[i][j]);
}
}
}
void Initialize(void)
{
string s;
ifstream f;
f.open("input_file", ios_base::in);
f >> s >> nmax;
f >> s >> dx;
f >> s >> dy;
f >> s >> imax;
f >> s >> jmax;
f >> s >> cfl;
f.close();
cout << "nmax = " << nmax << endl;
cout << "dx = " << dx << endl;
cout << "dy = " << dy << endl;
cout << "imax = " << imax << endl;
cout << "jmax = " << jmax << endl;
cout << "CFL = " << cfl << endl;
Q1 = Array2D(imax, jmax);
Q2 = Array2D(imax, jmax);
Q3 = Array2D(imax, jmax);
Q4 = Array2D(imax, jmax);
dF1 = Array2D(imax, jmax);
dF2 = Array2D(imax, jmax);
dF3 = Array2D(imax, jmax);
dF4 = Array2D(imax, jmax);
dG1 = Array2D(imax, jmax);
dG2 = Array2D(imax, jmax);
dG3 = Array2D(imax, jmax);
dG4 = Array2D(imax, jmax);
rho = Array2D(imax, jmax);
u = Array2D(imax, jmax);
v = Array2D(imax, jmax);
e = Array2D(imax, jmax);
p = Array2D(imax, jmax);
c = Array2D(imax, jmax);
float rho2 = 2.5;
float u2 = 0.0;
float v2 = 0.0;
float p2 = 2.0;
float e2 = p2 / (gamma - 1.0) / rho2;
float qq2[4];
qq2[0] = rho2;
qq2[1] = rho2 * u2;
qq2[2] = rho2 * v2;
qq2[3] = rho2 * (e2 + u2 * u2 / 2.0 + v2 * v2 / 2.0);
float rho1 = 1.5;
float u1 = 0.0;
float v1 = 0.0;
float p1 = 1.0;
float e1 = p1 / (gamma - 1.0) / rho1;
float qq1[4];
qq1[0] = rho1;
qq1[1] = rho1 * u1;
qq1[2] = rho1 * v1;
qq1[3] = rho1 * (e1 + u1 * u1 / 2.0 + v1 * v1 / 2.0);
float x, y, r;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
x = (float(i) + 0.5) * dx;
y = (float(j) + 0.5) * dy;
r = sqrt(x * x + y * y);
if (r < 0.2)
{
Q1[i][j] = qq2[0];
Q2[i][j] = qq2[1];
Q3[i][j] = qq2[2];
Q4[i][j] = qq2[3];
}
else
{
Q1[i][j] = qq1[0];
Q2[i][j] = qq1[1];
Q3[i][j] = qq1[2];
Q4[i][j] = qq1[3];
}
dF1[i][j] = 0.0;
dF2[i][j] = 0.0;
dF3[i][j] = 0.0;
dF4[i][j] = 0.0;
dG1[i][j] = 0.0;
dG2[i][j] = 0.0;
dG3[i][j] = 0.0;
dG4[i][j] = 0.0;
}
}
}
void calc_dt()
{
float vel = 0.0;
float upc;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
upc = sqrt(u[i][j] * u[i][j] + v[i][j] * v[i][j]) + c[i][j];
vel = (vel > upc) ? vel : upc;
}
}
dt = cfl * sqrt((dx * dx + dy * dy) / 2.0) / vel;
}
void calc_min_max()
{
float min_rho = 1.0e5;
float max_rho = 0.0;
float min_p = 1.0e5;
float max_p = 0.0;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
min_rho = (min_rho < rho[i][j]) ? min_rho : rho[i][j];
min_p = (min_p < p[i][j]) ? min_p : p[i][j];
max_rho = (max_rho > rho[i][j]) ? max_rho : rho[i][j];
max_p = (max_p > p[i][j]) ? max_p : p[i][j];
}
}
cout << "min/max rho = " << min_rho << " " << max_rho << endl;
cout << "min/max p = " << min_p << " " << max_p << endl;
}
void calc_flux_xdir_left_bc()
{
float left[6], right[6];
float fimh[4], fiph[4];
for (int j = 0; j < jmax; j++)
{
left[0] = rho[0][j];
left[1] = -u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[0][j];
right[1] = u[0][j];
right[2] = v[0][j];
right[3] = e[0][j];
right[4] = p[0][j];
right[5] = c[0][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[0][j];
left[1] = u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[1][j];
right[1] = u[1][j];
right[2] = v[1][j];
right[3] = e[1][j];
right[4] = p[1][j];
right[5] = c[1][j];
riemann_solver_xdir(left, right, fiph);
dF1[0][j] = (fiph[0] - fimh[0]) / dx;
dF2[0][j] = (fiph[1] - fimh[1]) / dx;
dF3[0][j] = (fiph[2] - fimh[2]) / dx;
dF4[0][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_right_bc()
{
float left[6], right[6];
float fimh[4], fiph[4];
for (int j = 0; j < jmax; j++)
{
left[0] = rho[imax - 2][j];
left[1] = u[imax - 2][j];
left[2] = v[imax - 2][j];
left[3] = e[imax - 2][j];
left[4] = p[imax - 2][j];
left[5] = c[imax - 2][j];
right[0] = rho[imax - 1][j];
right[1] = u[imax - 1][j];
right[2] = v[imax - 1][j];
right[3] = e[imax - 1][j];
right[4] = p[imax - 1][j];
right[5] = c[imax - 1][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[imax - 1][j];
left[1] = u[imax - 1][j];
left[2] = v[imax - 1][j];
left[3] = e[imax - 1][j];
left[4] = p[imax - 1][j];
left[5] = c[imax - 1][j];
right[0] = rho[imax - 1][j];
right[1] = u[imax - 1][j];
right[2] = v[imax - 1][j];
right[3] = e[imax - 1][j];
right[4] = p[imax - 1][j];
right[5] = c[imax - 1][j];
riemann_solver_xdir(left, right, fiph);
dF1[imax - 1][j] = (fiph[0] - fimh[0]) / dx;
dF2[imax - 1][j] = (fiph[1] - fimh[1]) / dx;
dF3[imax - 1][j] = (fiph[2] - fimh[2]) / dx;
dF4[imax - 1][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_internal()
{
float left[6], right[6];
float fimh[4], fiph[4];
for (int i = 1; i < imax - 1; i++)
{
for (int j = 0; j < jmax; j++)
{
left[0] = rho[i - 1][j];
left[1] = u[i - 1][j];
left[2] = v[i - 1][j];
left[3] = e[i - 1][j];
left[4] = p[i - 1][j];
left[5] = c[i - 1][j];
right[0] = rho[i][j];
right[1] = u[i][j];
right[2] = v[i][j];
right[3] = e[i][j];
right[4] = p[i][j];
right[5] = c[i][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[i][j];
left[1] = u[i][j];
left[2] = v[i][j];
left[3] = e[i][j];
left[4] = p[i][j];
left[5] = c[i][j];
right[0] = rho[i + 1][j];
right[1] = u[i + 1][j];
right[2] = v[i + 1][j];
right[3] = e[i + 1][j];
right[4] = p[i + 1][j];
right[5] = c[i + 1][j];
riemann_solver_xdir(left, right, fiph);
dF1[i][j] = (fiph[0] - fimh[0]) / dx;
dF2[i][j] = (fiph[1] - fimh[1]) / dx;
dF3[i][j] = (fiph[2] - fimh[2]) / dx;
dF4[i][j] = (fiph[3] - fimh[3]) / dx;
}
}
return;
}
void calc_flux_ydir_left_bc()
{
float left[6], right[6];
float gjmh[4], gjph[4];
for (int i = 0; i < imax; i++)
{
left[0] = rho[i][0];
left[1] = u[i][0];
left[2] = v[i][0];
left[3] = e[i][0];
left[4] = p[i][0];
left[5] = c[i][0];
right[0] = rho[i][0];
right[1] = u[i][0];
right[2] = v[i][0];
right[3] = e[i][0];
right[4] = p[i][0];
right[5] = c[i][0];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][0];
left[1] = u[i][0];
left[2] = v[i][0];
left[3] = e[i][0];
left[4] = p[i][0];
left[5] = c[i][0];
right[0] = rho[i][1];
right[1] = u[i][1];
right[2] = v[i][1];
right[3] = e[i][1];
right[4] = p[i][1];
right[5] = c[i][1];
riemann_solver_ydir(left, right, gjph);
dG1[i][0] = (gjph[0] - gjmh[0]) / dy;
dG2[i][0] = (gjph[1] - gjmh[1]) / dy;
dG3[i][0] = (gjph[2] - gjmh[2]) / dy;
dG4[i][0] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_right_bc()
{
float left[6], right[6];
float gjmh[4], gjph[4];
for (int i = 0; i < imax; i++)
{
left[0] = rho[i][jmax - 2];
left[1] = u[i][jmax - 2];
left[2] = v[i][jmax - 2];
left[3] = e[i][jmax - 2];
left[4] = p[i][jmax - 2];
left[5] = c[i][jmax - 2];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][jmax - 1];
left[1] = u[i][jmax - 1];
left[2] = v[i][jmax - 1];
left[3] = e[i][jmax - 1];
left[4] = p[i][jmax - 1];
left[5] = c[i][jmax - 1];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = -v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjph);
dG1[i][jmax - 1] = (gjph[0] - gjmh[0]) / dy;
dG2[i][jmax - 1] = (gjph[1] - gjmh[1]) / dy;
dG3[i][jmax - 1] = (gjph[2] - gjmh[2]) / dy;
dG4[i][jmax - 1] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_internal()
{
float left[6], right[6];
float gjmh[4], gjph[4];
for (int i = 0; i < imax; i++)
{
for (int j = 1; j < jmax - 1; j++)
{
left[0] = rho[i][j - 1];
left[1] = u[i][j - 1];
left[2] = v[i][j - 1];
left[3] = e[i][j - 1];
left[4] = p[i][j - 1];
left[5] = c[i][j - 1];
right[0] = rho[i][j];
right[1] = u[i][j];
right[2] = v[i][j];
right[3] = e[i][j];
right[4] = p[i][j];
right[5] = c[i][j];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][j];
left[1] = u[i][j];
left[2] = v[i][j];
left[3] = e[i][j];
left[4] = p[i][j];
left[5] = c[i][j];
right[0] = rho[i][j + 1];
right[1] = u[i][j + 1];
right[2] = v[i][j + 1];
right[3] = e[i][j + 1];
right[4] = p[i][j + 1];
right[5] = c[i][j + 1];
riemann_solver_ydir(left, right, gjph);
dG1[i][j] = (gjph[0] - gjmh[0]) / dy;
dG2[i][j] = (gjph[1] - gjmh[1]) / dy;
dG3[i][j] = (gjph[2] - gjmh[2]) / dy;
dG4[i][j] = (gjph[3] - gjmh[3]) / dy;
}
}
return;
}
void calc_flux()
{
//----------X direction --------------
calc_flux_xdir_left_bc();
calc_flux_xdir_right_bc();
calc_flux_xdir_internal();
//----------Y direction --------------
calc_flux_ydir_left_bc();
calc_flux_ydir_right_bc();
calc_flux_ydir_internal();
return;
}
void Update(void)
{
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
Q1[i][j] -= dt * (dF1[i][j] + dG1[i][j]);
Q2[i][j] -= dt * (dF2[i][j] + dG2[i][j]);
Q3[i][j] -= dt * (dF3[i][j] + dG3[i][j]);
Q4[i][j] -= dt * (dF4[i][j] + dG4[i][j]);
}
}
}
void Evolve(void)
{
for (int n = 0; n < nmax; n++)
{
cons_2_prim();
calc_dt();
cout << "*****************" << endl;
cout << "n = " << n << "; dt = " << dt << " \n";
calc_min_max();
calc_flux();
Update();
}
cout << "*****************" << endl;
}
void Output_results()
{
float x, y, vel;
ofstream f;
f.open("rho_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << rho[i][j] << endl;
}
}
f.close();
f.open("p_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << p[i][j] << endl;
}
}
f.close();
f.open("vel_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
vel = sqrt(pow(u[i][j], 2.0) + pow(v[i][j], 2.0));
f << vel << endl;
}
}
f.close();
}
int main()
{
double t_start[10], t_end[10];
t_start[0] = mytimer();
Initialize();
Evolve();
t_end[0] = mytimer();
Output_results();
cout << "Total time:" << t_end[0] - t_start[0] << endl;
return 0;
}- cfd2d.h
#ifndef CFD2D_H_
#define CFD2D_H_
float **Array2D(unsigned int row, unsigned int col);
void cons_2_prim();
void Initialize(void);
void calc_dt();
void calc_min_max();
void calc_flux_xdir_left_bc();
void calc_flux_xdir_right_bc();
void calc_flux_xdir_internal();
void calc_flux_ydir_left_bc();
void calc_flux_ydir_right_bc();
void calc_flux_ydir_internal();
void calc_flux();
void Update(void);
void Evolve(void);
void Output_results();
#endif
- input_file
nmax 500
dx 1.0e-2
dy 1.0e-2
imax 1026
jmax 1026
CFL 0.5
- Makefile
CXX = g++ -O1 -ftree-vectorize
TARGET = eular.seq.exe
SOURCE = $(wildcard *.cpp)
OBJECT = $(patsubst %.cpp, %.o, $(SOURCE))
$(TARGET): $(OBJECT)
$(CXX) $^ -o $@
%.o: %.cpp
$(CXX) -c $< -o $@
.PHONY: clean
clean:
rm -rf $(TARGET) *.o- riemann.cpp
#include "riemann.h"
void riemann_solver_xdir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(ul) + cl;
wr = abs(ur) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * ul;
fl[1] = rhol * ul * ul + pl * 1.0;
fl[2] = rhol * ul * vl;
fl[3] = rhol * ul * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * ul;
// right flux
fr[0] = rhor * ur;
fr[1] = rhor * ur * ur + pr * 1.0;
fr[2] = rhor * ur * vr;
fr[3] = rhor * ur * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * ur;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}
void riemann_solver_ydir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(vl) + cl;
wr = abs(vr) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * vl;
fl[1] = rhol * vl * ul;
fl[2] = rhol * vl * vl + pl * 1.0;
fl[3] = rhol * vl * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * vl;
// right flux
fr[0] = rhor * vr;
fr[1] = rhor * vr * ur;
fr[2] = rhor * vr * vr + pr * 1.0;
fr[3] = rhor * vr * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * vr;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}- riemann.h
#ifndef RIEMANN_H_
#define RIEMANN_H_
#include <stdlib.h>
void riemann_solver_xdir(float left[], float right[], float flux[]);
void riemann_solver_ydir(float left[], float right[], float flux[]);
#endif- timer.cpp
#include <cstdlib>
#include <sys/time.h>
#include <sys/resource.h>
#include "timer.h"
double mytimer(void)
{
struct timeval tp;
static long start=0, startu;
if (!start)
{
gettimeofday(&tp, NULL);
start = tp.tv_sec;
startu = tp.tv_usec;
return 0.0;
}
gettimeofday(&tp, NULL);
return ((double)(tp.tv_sec - start)) + (tp.tv_usec - startu) / 1000000.0;
}- timer.h
#ifndef TIMER_H_
#define TIMER_H_
double mytimer(void);
#endif- yhrun.sh
#!/bin/bash
yhrun -p thcp1 -N 1 -n 1 eular.seq.exe运行结果如下：
2、使用 OpenMP 优化
代码目录如下：
- cfd2d.cpp
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include "riemann.h"
#include "timer.h"
using namespace std;
#define gamma 2.1
float **Q1, **Q2, **Q3, **Q4, **dF1, **dF2, **dF3, **dF4, **dG1, **dG2, **dG3, **dG4;
float dx, dy, dt;
float **rho, **u, **v, **e, **p, **c;
int nmax, imax, jmax;
float cfl;
float **Array2D(unsigned int row, unsigned int col)
{
float **x;
x = new float *[row];
for (int count = 0; count < row; count++)
{
x[count] = new float[col];
}
return x;
}
void cons_2_prim()
{
#pragma omp parallel for collapse(2)
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
rho[i][j] = Q1[i][j];
u[i][j] = Q2[i][j] / Q1[i][j];
v[i][j] = Q3[i][j] / Q1[i][j];
float ke = u[i][j] * u[i][j] / 2.0 + v[i][j] * v[i][j] / 2.0;
e[i][j] = Q4[i][j] / Q1[i][j] - ke;
p[i][j] = (gamma - 1.0) * rho[i][j] * e[i][j];
c[i][j] = sqrt(gamma * p[i][j] / rho[i][j]);
}
}
}
void Initialize(void)
{
string s;
ifstream f;
f.open("input_file", ios_base::in);
f >> s >> nmax;
f >> s >> dx;
f >> s >> dy;
f >> s >> imax;
f >> s >> jmax;
f >> s >> cfl;
f.close();
cout << "nmax = " << nmax << endl;
cout << "dx = " << dx << endl;
cout << "dy = " << dy << endl;
cout << "imax = " << imax << endl;
cout << "jmax = " << jmax << endl;
cout << "CFL = " << cfl << endl;
Q1 = Array2D(imax, jmax);
Q2 = Array2D(imax, jmax);
Q3 = Array2D(imax, jmax);
Q4 = Array2D(imax, jmax);
dF1 = Array2D(imax, jmax);
dF2 = Array2D(imax, jmax);
dF3 = Array2D(imax, jmax);
dF4 = Array2D(imax, jmax);
dG1 = Array2D(imax, jmax);
dG2 = Array2D(imax, jmax);
dG3 = Array2D(imax, jmax);
dG4 = Array2D(imax, jmax);
rho = Array2D(imax, jmax);
u = Array2D(imax, jmax);
v = Array2D(imax, jmax);
e = Array2D(imax, jmax);
p = Array2D(imax, jmax);
c = Array2D(imax, jmax);
float rho2 = 2.5;
float u2 = 0.0;
float v2 = 0.0;
float p2 = 2.0;
float e2 = p2 / (gamma - 1.0) / rho2;
float qq2[4];
qq2[0] = rho2;
qq2[1] = rho2 * u2;
qq2[2] = rho2 * v2;
qq2[3] = rho2 * (e2 + u2 * u2 / 2.0 + v2 * v2 / 2.0);
float rho1 = 1.5;
float u1 = 0.0;
float v1 = 0.0;
float p1 = 1.0;
float e1 = p1 / (gamma - 1.0) / rho1;
float qq1[4];
qq1[0] = rho1;
qq1[1] = rho1 * u1;
qq1[2] = rho1 * v1;
qq1[3] = rho1 * (e1 + u1 * u1 / 2.0 + v1 * v1 / 2.0);
float x, y, r;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
x = (float(i) + 0.5) * dx;
y = (float(j) + 0.5) * dy;
r = sqrt(x * x + y * y);
if (r < 0.2)
{
Q1[i][j] = qq2[0];
Q2[i][j] = qq2[1];
Q3[i][j] = qq2[2];
Q4[i][j] = qq2[3];
}
else
{
Q1[i][j] = qq1[0];
Q2[i][j] = qq1[1];
Q3[i][j] = qq1[2];
Q4[i][j] = qq1[3];
}
dF1[i][j] = 0.0;
dF2[i][j] = 0.0;
dF3[i][j] = 0.0;
dF4[i][j] = 0.0;
dG1[i][j] = 0.0;
dG2[i][j] = 0.0;
dG3[i][j] = 0.0;
dG4[i][j] = 0.0;
}
}
}
void calc_dt()
{
float vel = 0.0;
float upc;
#pragma omp parallel for private(upc) reduction(max:vel) collapse(2)
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
upc = sqrt(u[i][j] * u[i][j] + v[i][j] * v[i][j]) + c[i][j];
vel = (vel > upc) ? vel : upc;
}
}
dt = cfl * sqrt((dx * dx + dy * dy) / 2.0) / vel;
}
void calc_min_max()
{
float min_rho = 1.0e5;
float max_rho = 0.0;
float min_p = 1.0e5;
float max_p = 0.0;
#pragma omp parallel for reduction(max:max_rho,max_p) reduction(min:min_rho,min_p)
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
min_rho = (min_rho < rho[i][j]) ? min_rho : rho[i][j];
min_p = (min_p < p[i][j]) ? min_p : p[i][j];
max_rho = (max_rho > rho[i][j]) ? max_rho : rho[i][j];
max_p = (max_p > p[i][j]) ? max_p : p[i][j];
}
}
cout << "min/max rho = " << min_rho << " " << max_rho << endl;
cout << "min/max p = " << min_p << " " << max_p << endl;
}
void calc_flux_xdir_left_bc()
{
float left[6], right[6];
float fimh[4], fiph[4];
#pragma omp parallel for private(left, right, fimh, fiph)
for (int j = 0; j < jmax; j++)
{
left[0] = rho[0][j];
left[1] = -u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[0][j];
right[1] = u[0][j];
right[2] = v[0][j];
right[3] = e[0][j];
right[4] = p[0][j];
right[5] = c[0][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[0][j];
left[1] = u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[1][j];
right[1] = u[1][j];
right[2] = v[1][j];
right[3] = e[1][j];
right[4] = p[1][j];
right[5] = c[1][j];
riemann_solver_xdir(left, right, fiph);
dF1[0][j] = (fiph[0] - fimh[0]) / dx;
dF2[0][j] = (fiph[1] - fimh[1]) / dx;
dF3[0][j] = (fiph[2] - fimh[2]) / dx;
dF4[0][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_right_bc()
{
float fimh[4], fiph[4];
#pragma omp parallel for private(fimh, fiph)
for (int j = 0; j < jmax; j++)
{
riemann_solver_xdir_direct(imax - 2, imax - 1, j, fimh);
riemann_solver_xdir_direct(imax - 1, imax - 1, j, fiph);
dF1[imax - 1][j] = (fiph[0] - fimh[0]) / dx;
dF2[imax - 1][j] = (fiph[1] - fimh[1]) / dx;
dF3[imax - 1][j] = (fiph[2] - fimh[2]) / dx;
dF4[imax - 1][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_internal()
{
float fimh[4], fiph[4];
#pragma omp parallel for private(fimh, fiph) collapse(2)
for (int i = 1; i < imax - 1; i++)
{
for (int j = 0; j < jmax; j++)
{
riemann_solver_xdir_direct(i - 1, i, j, fimh);
riemann_solver_xdir_direct(i, i + 1, j, fiph);
dF1[i][j] = (fiph[0] - fimh[0]) / dx;
dF2[i][j] = (fiph[1] - fimh[1]) / dx;
dF3[i][j] = (fiph[2] - fimh[2]) / dx;
dF4[i][j] = (fiph[3] - fimh[3]) / dx;
}
}
return;
}
void calc_flux_ydir_left_bc()
{
float gjmh[4], gjph[4];
#pragma omp parallel for private(gjmh, gjph)
for (int i = 0; i < imax; i++)
{
riemann_solver_ydir_direct(i, 0, 0, gjmh);
riemann_solver_ydir_direct(i, 0, 1, gjph);
dG1[i][0] = (gjph[0] - gjmh[0]) / dy;
dG2[i][0] = (gjph[1] - gjmh[1]) / dy;
dG3[i][0] = (gjph[2] - gjmh[2]) / dy;
dG4[i][0] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_right_bc()
{
float left[6], right[6];
float gjmh[4], gjph[4];
#pragma omp parallel for private(left, right, gjmh, gjph)
for (int i = 0; i < imax; i++)
{
left[0] = rho[i][jmax - 2];
left[1] = u[i][jmax - 2];
left[2] = v[i][jmax - 2];
left[3] = e[i][jmax - 2];
left[4] = p[i][jmax - 2];
left[5] = c[i][jmax - 2];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][jmax - 1];
left[1] = u[i][jmax - 1];
left[2] = v[i][jmax - 1];
left[3] = e[i][jmax - 1];
left[4] = p[i][jmax - 1];
left[5] = c[i][jmax - 1];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = -v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjph);
dG1[i][jmax - 1] = (gjph[0] - gjmh[0]) / dy;
dG2[i][jmax - 1] = (gjph[1] - gjmh[1]) / dy;
dG3[i][jmax - 1] = (gjph[2] - gjmh[2]) / dy;
dG4[i][jmax - 1] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_internal()
{
float gjmh[4], gjph[4];
#pragma omp parallel for private(gjmh, gjph) collapse(2)
for (int i = 0; i < imax; i++)
{
for (int j = 1; j < jmax - 1; j++)
{
riemann_solver_ydir_direct(i, j - 1, j, gjmh);
riemann_solver_ydir_direct(i, j, j + 1, gjph);
dG1[i][j] = (gjph[0] - gjmh[0]) / dy;
dG2[i][j] = (gjph[1] - gjmh[1]) / dy;
dG3[i][j] = (gjph[2] - gjmh[2]) / dy;
dG4[i][j] = (gjph[3] - gjmh[3]) / dy;
}
}
return;
}
void calc_flux()
{
//----------X direction --------------
calc_flux_xdir_left_bc();
calc_flux_xdir_right_bc();
calc_flux_xdir_internal();
//----------Y direction --------------
calc_flux_ydir_left_bc();
calc_flux_ydir_right_bc();
calc_flux_ydir_internal();
return;
}
void Update(void)
{
#pragma omp parallel for collapse(2)
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
Q1[i][j] -= dt * (dF1[i][j] + dG1[i][j]);
Q2[i][j] -= dt * (dF2[i][j] + dG2[i][j]);
Q3[i][j] -= dt * (dF3[i][j] + dG3[i][j]);
Q4[i][j] -= dt * (dF4[i][j] + dG4[i][j]);
}
}
}
void Evolve(void)
{
for (int n = 0; n < nmax; n++)
{
cons_2_prim();
calc_dt();
cout << "*****************" << endl;
cout << "n = " << n << "; dt = " << dt << " \n";
calc_min_max();
calc_flux();
Update();
}
cout << "*****************" << endl;
}
void Output_results()
{
float x, y, vel;
ofstream f;
f.open("rho_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << rho[i][j] << endl;
}
}
f.close();
f.open("p_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << p[i][j] << endl;
}
}
f.close();
f.open("vel_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
vel = sqrt(pow(u[i][j], 2.0) + pow(v[i][j], 2.0));
f << vel << endl;
}
}
f.close();
}
int main()
{
double t_start[10], t_end[10];
t_start[0] = mytimer();
Initialize();
Evolve();
t_end[0] = mytimer();
Output_results();
cout << "Total time:" << t_end[0] - t_start[0] << endl;
return 0;
}- cfd2d.h
#ifndef CFD2D_H_
#define CFD2D_H_
float **Array2D(unsigned int row, unsigned int col);
void cons_2_prim();
void Initialize(void);
void calc_dt();
void calc_min_max();
void calc_flux_xdir_left_bc();
void calc_flux_xdir_right_bc();
void calc_flux_xdir_internal();
void calc_flux_ydir_left_bc();
void calc_flux_ydir_right_bc();
void calc_flux_ydir_internal();
void calc_flux();
void Update(void);
void Evolve(void);
void Output_results();
#endif
- input_file
nmax 500
dx 1.0e-2
dy 1.0e-2
imax 1026
jmax 1026
CFL 0.5
- Makefile
CXX = g++ -O1 -ftree-vectorize -fopenmp
TARGET = eular.opt_omp.exe
SOURCE = $(wildcard *.cpp)
OBJECT = $(patsubst %.cpp, %.o, $(SOURCE))
$(TARGET): $(OBJECT)
$(CXX) $^ -o $@
%.o: %.cpp
$(CXX) -c $< -o $@
.PHONY: clean
clean:
rm -rf $(TARGET) *.o- omp_strong_scal_test.sh
#!/bin/bash
rootPath=`pwd`
workZone=omp_strong_scal_test
exeBin=eular.opt_omp.exe
inputFile=input_file
for nthreads in 1 2 4 6 12 24 32
do
workDir=${workZone}/nthreads${nthreads}
mkdir -p ${workDir}
cp -rf ${exeBin} ${inputFile} ${workDir}
cd ${workDir}
cat > yhrun.sh << EOF
#!/bin/bash
export OMP_NUM_THREADS=${nthreads}
yhrun -p thcp1 -N 1 -n 1 -c ${nthreads} ${exeBin}
EOF
yhbatch -p thcp1 yhrun.sh
cd ${rootPath}
done
- riemann.cpp
#include "riemann.h"
void riemann_solver_xdir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(ul) + cl;
wr = abs(ur) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * ul;
fl[1] = rhol * ul * ul + pl * 1.0;
fl[2] = rhol * ul * vl;
fl[3] = rhol * ul * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * ul;
// right flux
fr[0] = rhor * ur;
fr[1] = rhor * ur * ur + pr * 1.0;
fr[2] = rhor * ur * vr;
fr[3] = rhor * ur * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * ur;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}
void riemann_solver_ydir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(vl) + cl;
wr = abs(vr) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * vl;
fl[1] = rhol * vl * ul;
fl[2] = rhol * vl * vl + pl * 1.0;
fl[3] = rhol * vl * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * vl;
// right flux
fr[0] = rhor * vr;
fr[1] = rhor * vr * ur;
fr[2] = rhor * vr * vr + pr * 1.0;
fr[3] = rhor * vr * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * vr;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}
void riemann_solver_xdir_direct(int il, int ir, int j, float flux[])
{
float rhol = rho[il][j];
float ul = u[il][j];
float vl = v[il][j];
float el = e[il][j];
float pl = p[il][j];
float cl = c[il][j];
float rhor = rho[ir][j];
float ur = u[ir][j];
float vr = v[ir][j];
float er = e[ir][j];
float pr = p[ir][j];
float cr = c[ir][j];
float cwave, wl, wr;
wl = abs(ul) + cl;
wr = abs(ur) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * ul;
fl[1] = rhol * ul * ul + pl * 1.0;
fl[2] = rhol * ul * vl;
fl[3] = rhol * ul * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * ul;
// right flux
fr[0] = rhor * ur;
fr[1] = rhor * ur * ur + pr * 1.0;
fr[2] = rhor * ur * vr;
fr[3] = rhor * ur * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * ur;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}
void riemann_solver_ydir_direct(int i, int jl, int jr, float flux[])
{
float rhol = rho[i][jl];
float ul = u[i][jl];
float vl = v[i][jl];
float el = e[i][jl];
float pl = p[i][jl];
float cl = c[i][jl];
float rhor = rho[i][jr];
float ur = u[i][jr];
float vr = v[i][jr];
float er = e[i][jr];
float pr = p[i][jr];
float cr = c[i][jr];
float cwave, wl, wr;
wl = abs(vl) + cl;
wr = abs(vr) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * vl;
fl[1] = rhol * vl * ul;
fl[2] = rhol * vl * vl + pl * 1.0;
fl[3] = rhol * vl * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * vl;
// right flux
fr[0] = rhor * vr;
fr[1] = rhor * vr * ur;
fr[2] = rhor * vr * vr + pr * 1.0;
fr[3] = rhor * vr * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * vr;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}- riemann.h
#ifndef RIEMANN_H_
#define RIEMANN_H_
#include <stdlib.h>
extern float **Q1, **Q2, **Q3, **Q4, **dF1, **dF2, **dF3, **dF4, **dG1, **dG2, **dG3, **dG4;
extern float dx, dy, dt;
extern float **rho, **u, **v, **e, **p, **c;
extern int nmax, imax, jmax;
extern float cfl;
void riemann_solver_xdir(float left[], float right[], float flux[]);
void riemann_solver_ydir(float left[], float right[], float flux[]);
void riemann_solver_xdir_direct(int il, int ir, int j, float flux[]);
void riemann_solver_ydir_direct(int i, int jl, int jr, float flux[]);
#endif- timer.cpp
#include <cstdlib>
#include <sys/time.h>
#include <sys/resource.h>
#include "timer.h"
double mytimer(void)
{
struct timeval tp;
static long start=0, startu;
if (!start)
{
gettimeofday(&tp, NULL);
start = tp.tv_sec;
startu = tp.tv_usec;
return 0.0;
}
gettimeofday(&tp, NULL);
return ((double)(tp.tv_sec - start)) + (tp.tv_usec - startu) / 1000000.0;
}- timer.h
#ifndef TIMER_H_
#define TIMER_H_
double mytimer(void);
#endif运行结果如下：加速比为 24.33，并行效率为 0.76，效果显著。
3、使用 SIMD 优化
代码目录如下：
- cfd2d.cpp
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include "riemann.h"
#include "timer.h"
#include <arm_neon.h>
using namespace std;
#define gamma 2.1
float **Q1, **Q2, **Q3, **Q4, **dF1, **dF2, **dF3, **dF4, **dG1, **dG2, **dG3, **dG4;
float dx, dy, dt;
float **rho, **u, **v, **e, **p, **c;
int nmax, imax, jmax;
float cfl;
float **Array2D(unsigned int row, unsigned int col)
{
float **x;
x = new float *[row];
for (unsigned int count = 0; count < row; count++)
{
x[count] = new float[col];
}
return x;
}
void cons_2_prim()
{
int jmax4 = jmax / 4 * 4;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax4; j += 4)
{
float32x4_t vrho = vld1q_f32(&Q1[i][j]);
float32x4_t vu = vdivq_f32(vld1q_f32(&Q2[i][j]), vrho);
float32x4_t vv = vdivq_f32(vld1q_f32(&Q3[i][j]), vrho);
float32x4_t vke = vmulq_n_f32(vaddq_f32(vmulq_f32(vu, vu), vmulq_f32(vv, vv)), 0.5);
float32x4_t ve = vsubq_f32(vdivq_f32(vld1q_f32(&Q4[i][j]), vrho), vke);
float32x4_t vp = vmulq_n_f32(vmulq_f32(vrho, ve), (gamma - 1.0));
float32x4_t vc = vsqrtq_f32(vdivq_f32(vmulq_n_f32(vp, gamma), vrho));
vst1q_f32(&rho[i][j], vrho);
vst1q_f32(&u[i][j], vu);
vst1q_f32(&v[i][j], vv);
vst1q_f32(&e[i][j], ve);
vst1q_f32(&p[i][j], vp);
vst1q_f32(&c[i][j], vc);
}
for (int j = jmax4; j < jmax; j++)
{
rho[i][j] = Q1[i][j];
u[i][j] = Q2[i][j] / Q1[i][j];
v[i][j] = Q3[i][j] / Q1[i][j];
float ke = u[i][j] * u[i][j] / 2.0 + v[i][j] * v[i][j] / 2.0;
e[i][j] = Q4[i][j] / Q1[i][j] - ke;
p[i][j] = (gamma - 1.0) * rho[i][j] * e[i][j];
c[i][j] = sqrt(gamma * p[i][j] / rho[i][j]);
}
}
}
void Initialize(void)
{
string s;
ifstream f;
f.open("input_file", ios_base::in);
f >> s >> nmax;
f >> s >> dx;
f >> s >> dy;
f >> s >> imax;
f >> s >> jmax;
f >> s >> cfl;
f.close();
cout << "nmax = " << nmax << endl;
cout << "dx = " << dx << endl;
cout << "dy = " << dy << endl;
cout << "imax = " << imax << endl;
cout << "jmax = " << jmax << endl;
cout << "CFL = " << cfl << endl;
Q1 = Array2D(imax, jmax);
Q2 = Array2D(imax, jmax);
Q3 = Array2D(imax, jmax);
Q4 = Array2D(imax, jmax);
dF1 = Array2D(imax, jmax);
dF2 = Array2D(imax, jmax);
dF3 = Array2D(imax, jmax);
dF4 = Array2D(imax, jmax);
dG1 = Array2D(imax, jmax);
dG2 = Array2D(imax, jmax);
dG3 = Array2D(imax, jmax);
dG4 = Array2D(imax, jmax);
rho = Array2D(imax, jmax);
u = Array2D(imax, jmax);
v = Array2D(imax, jmax);
e = Array2D(imax, jmax);
p = Array2D(imax, jmax);
c = Array2D(imax, jmax);
float rho2 = 2.5;
float u2 = 0.0;
float v2 = 0.0;
float p2 = 2.0;
float e2 = p2 / (gamma - 1.0) / rho2;
float qq2[4];
qq2[0] = rho2;
qq2[1] = rho2 * u2;
qq2[2] = rho2 * v2;
qq2[3] = rho2 * (e2 + u2 * u2 / 2.0 + v2 * v2 / 2.0);
float rho1 = 1.5;
float u1 = 0.0;
float v1 = 0.0;
float p1 = 1.0;
float e1 = p1 / (gamma - 1.0) / rho1;
float qq1[4];
qq1[0] = rho1;
qq1[1] = rho1 * u1;
qq1[2] = rho1 * v1;
qq1[3] = rho1 * (e1 + u1 * u1 / 2.0 + v1 * v1 / 2.0);
float x, y, r;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
x = (float(i) + 0.5) * dx;
y = (float(j) + 0.5) * dy;
r = sqrt(x * x + y * y);
if (r < 0.2)
{
Q1[i][j] = qq2[0];
Q2[i][j] = qq2[1];
Q3[i][j] = qq2[2];
Q4[i][j] = qq2[3];
}
else
{
Q1[i][j] = qq1[0];
Q2[i][j] = qq1[1];
Q3[i][j] = qq1[2];
Q4[i][j] = qq1[3];
}
dF1[i][j] = 0.0;
dF2[i][j] = 0.0;
dF3[i][j] = 0.0;
dF4[i][j] = 0.0;
dG1[i][j] = 0.0;
dG2[i][j] = 0.0;
dG3[i][j] = 0.0;
dG4[i][j] = 0.0;
}
}
}
void calc_dt()
{
float vel, upc;
int jmax4 = jmax / 4 * 4;
float32x4_t vvel = {0.0, 0.0, 0.0, 0.0};
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax4; j += 4)
{
float32x4_t vu = vld1q_f32(&u[i][j]);
float32x4_t vv = vld1q_f32(&v[i][j]);
float32x4_t vupc = vaddq_f32(vsqrtq_f32(vaddq_f32(vmulq_f32(vu, vu), vmulq_f32(vv, vv))), vld1q_f32(&c[i][j]));
// 方法一
vvel = vmaxq_f32(vvel, vupc);
// 方法二
// uint32x4_t vflag = vcgtq_f32(vvel, vupc);
// vvel = vbslq_f32(vflag, vvel, vupc);
}
// 方法一
vel = vmaxvq_f32(vvel);
// 方法二
// float32x2_t va = vpmax_f32(vget_low_f32(vvel), vget_high_f32(vvel));
// float32x2_t vb = vpmax_f32(va, va);
// vel = vget_lane_f32(vb, 0);
for (int j = jmax4; j < jmax; j++)
{
upc = sqrt(u[i][j] * u[i][j] + v[i][j] * v[i][j]) + c[i][j];
vel = (vel > upc) ? vel : upc;
}
}
dt = cfl * sqrt((dx * dx + dy * dy) / 2.0) / vel;
}
void calc_min_max()
{
float min_rho = 1.0e5;
float max_rho = 0.0;
float min_p = 1.0e5;
float max_p = 0.0;
uint32x4_t vflag;
int jmax4 = jmax / 4 * 4;
float32x4_t vmin_rho = {1.0e5, 1.0e5, 1.0e5, 1.0e5};
float32x4_t vmax_rho = {0.0, 0.0, 0.0, 0.0};
float32x4_t vmin_p = {1.0e5, 1.0e5, 1.0e5, 1.0e5};
float32x4_t vmax_p = {0.0, 0.0, 0.0, 0.0};
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax4; j += 4)
{
float32x4_t vrho = vld1q_f32(&rho[i][j]);
float32x4_t vp = vld1q_f32(&p[i][j]);
// 方法一
vmin_rho = vminq_f32(vmin_rho, vrho);
vmin_p = vminq_f32(vmin_p, vp);
vmax_rho = vmaxq_f32(vmax_rho, vrho);
vmax_p = vmaxq_f32(vmax_p, vp);
// 方法二
// vflag = vcltq_f32(vmin_rho, vrho);
// vmin_rho = vbslq_f32(vflag, vmin_rho, vrho);
// vflag = vcltq_f32(vmin_p, vp);
// vmin_p = vbslq_f32(vflag, vmin_p, vp);
// vflag = vcgtq_f32(vmax_rho, vrho);
// vmax_rho = vbslq_f32(vflag, vmax_rho, vrho);
// vflag = vcgtq_f32(vmax_p, vp);
// vmax_p = vbslq_f32(vflag, vmax_p, vp);
}
// 方法一
min_rho = vminvq_f32(vmin_rho);
min_p = vminvq_f32(vmin_p);
max_rho = vmaxvq_f32(vmax_rho);
max_p = vmaxvq_f32(vmax_p);
// 方法二
// float32x2_t va, vb;
// va = vpmin_f32(vget_low_f32(vmin_rho), vget_high_f32(vmin_rho));
// vb = vpmin_f32(va, va);
// min_rho = vget_lane_f32(vb, 0);
// va = vpmin_f32(vget_low_f32(vmin_p), vget_high_f32(vmin_p));
// vb = vpmin_f32(va, va);
// min_p = vget_lane_f32(vb, 0);
// va = vpmax_f32(vget_low_f32(vmax_rho), vget_high_f32(vmax_rho));
// vb = vpmax_f32(va, va);
// max_rho = vget_lane_f32(vb, 0);
// va = vpmax_f32(vget_low_f32(vmax_p), vget_high_f32(vmax_p));
// vb = vpmax_f32(va, va);
// max_p = vget_lane_f32(vb, 0);
for (int j = jmax4; j < jmax; j++)
{
min_rho = (min_rho < rho[i][j]) ? min_rho : rho[i][j];
min_p = (min_p < p[i][j]) ? min_p : p[i][j];
max_rho = (max_rho > rho[i][j]) ? max_rho : rho[i][j];
max_p = (max_p > p[i][j]) ? max_p : p[i][j];
}
}
cout << "min/max rho = " << min_rho << " " << max_rho << endl;
cout << "min/max p = " << min_p << " " << max_p << endl;
}
void calc_flux_xdir_left_bc()
{
float left[6], right[6];
float fimh[4], fiph[4];
for (int j = 0; j < jmax; j++)
{
left[0] = rho[0][j];
left[1] = -u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[0][j];
right[1] = u[0][j];
right[2] = v[0][j];
right[3] = e[0][j];
right[4] = p[0][j];
right[5] = c[0][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[0][j];
left[1] = u[0][j];
left[2] = v[0][j];
left[3] = e[0][j];
left[4] = p[0][j];
left[5] = c[0][j];
right[0] = rho[1][j];
right[1] = u[1][j];
right[2] = v[1][j];
right[3] = e[1][j];
right[4] = p[1][j];
right[5] = c[1][j];
riemann_solver_xdir(left, right, fiph);
dF1[0][j] = (fiph[0] - fimh[0]) / dx;
dF2[0][j] = (fiph[1] - fimh[1]) / dx;
dF3[0][j] = (fiph[2] - fimh[2]) / dx;
dF4[0][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_right_bc()
{
float left[6], right[6];
float fimh[4], fiph[4];
float32x4_t vfimh[4], vfiph[4];
float coef = 1.0 / dx;
int jmax4 = jmax / 4 * 4;
for (int j = 0; j < jmax4; j++)
{
riemann_solver_xdir_direct(imax - 2, imax - 1, j, vfimh);
riemann_solver_xdir_direct(imax - 1, imax - 1, j, vfiph);
vst1q_f32(&dF1[imax - 1][j], vmulq_n_f32(vsubq_f32(vfiph[0], vfimh[0]), coef));
vst1q_f32(&dF2[imax - 1][j], vmulq_n_f32(vsubq_f32(vfiph[1], vfimh[1]), coef));
vst1q_f32(&dF3[imax - 1][j], vmulq_n_f32(vsubq_f32(vfiph[2], vfimh[2]), coef));
vst1q_f32(&dF4[imax - 1][j], vmulq_n_f32(vsubq_f32(vfiph[3], vfimh[3]), coef));
}
for (int j = jmax4; j < jmax; j++)
{
left[0] = rho[imax - 2][j];
left[1] = u[imax - 2][j];
left[2] = v[imax - 2][j];
left[3] = e[imax - 2][j];
left[4] = p[imax - 2][j];
left[5] = c[imax - 2][j];
right[0] = rho[imax - 1][j];
right[1] = u[imax - 1][j];
right[2] = v[imax - 1][j];
right[3] = e[imax - 1][j];
right[4] = p[imax - 1][j];
right[5] = c[imax - 1][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[imax - 1][j];
left[1] = u[imax - 1][j];
left[2] = v[imax - 1][j];
left[3] = e[imax - 1][j];
left[4] = p[imax - 1][j];
left[5] = c[imax - 1][j];
right[0] = rho[imax - 1][j];
right[1] = u[imax - 1][j];
right[2] = v[imax - 1][j];
right[3] = e[imax - 1][j];
right[4] = p[imax - 1][j];
right[5] = c[imax - 1][j];
riemann_solver_xdir(left, right, fiph);
dF1[imax - 1][j] = (fiph[0] - fimh[0]) / dx;
dF2[imax - 1][j] = (fiph[1] - fimh[1]) / dx;
dF3[imax - 1][j] = (fiph[2] - fimh[2]) / dx;
dF4[imax - 1][j] = (fiph[3] - fimh[3]) / dx;
}
return;
}
void calc_flux_xdir_internal()
{
float left[6], right[6];
float fimh[4], fiph[4];
float32x4_t vfimh[4], vfiph[4];
float coef = 1.0 / dx;
int jmax4 = jmax / 4 * 4;
for (int i = 1; i < imax - 1; i++)
{
for (int j = 0; j < jmax4; j += 4)
{
riemann_solver_xdir_direct(i - 1, i, j, vfimh);
riemann_solver_xdir_direct(i, i + 1, j, vfiph);
vst1q_f32(&dF1[i][j], vmulq_n_f32(vsubq_f32(vfiph[0], vfimh[0]), coef));
vst1q_f32(&dF2[i][j], vmulq_n_f32(vsubq_f32(vfiph[1], vfimh[1]), coef));
vst1q_f32(&dF3[i][j], vmulq_n_f32(vsubq_f32(vfiph[2], vfimh[2]), coef));
vst1q_f32(&dF4[i][j], vmulq_n_f32(vsubq_f32(vfiph[3], vfimh[3]), coef));
}
for (int j = jmax4; j < jmax; j++)
{
left[0] = rho[i - 1][j];
left[1] = u[i - 1][j];
left[2] = v[i - 1][j];
left[3] = e[i - 1][j];
left[4] = p[i - 1][j];
left[5] = c[i - 1][j];
right[0] = rho[i][j];
right[1] = u[i][j];
right[2] = v[i][j];
right[3] = e[i][j];
right[4] = p[i][j];
right[5] = c[i][j];
riemann_solver_xdir(left, right, fimh);
left[0] = rho[i][j];
left[1] = u[i][j];
left[2] = v[i][j];
left[3] = e[i][j];
left[4] = p[i][j];
left[5] = c[i][j];
right[0] = rho[i + 1][j];
right[1] = u[i + 1][j];
right[2] = v[i + 1][j];
right[3] = e[i + 1][j];
right[4] = p[i + 1][j];
right[5] = c[i + 1][j];
riemann_solver_xdir(left, right, fiph);
dF1[i][j] = (fiph[0] - fimh[0]) / dx;
dF2[i][j] = (fiph[1] - fimh[1]) / dx;
dF3[i][j] = (fiph[2] - fimh[2]) / dx;
dF4[i][j] = (fiph[3] - fimh[3]) / dx;
}
}
return;
}
void calc_flux_ydir_left_bc()
{
float left[6], right[6];
float gjmh[4], gjph[4];
for (int i = 0; i < imax; i++)
{
left[0] = rho[i][0];
left[1] = u[i][0];
left[2] = v[i][0];
left[3] = e[i][0];
left[4] = p[i][0];
left[5] = c[i][0];
right[0] = rho[i][0];
right[1] = u[i][0];
right[2] = v[i][0];
right[3] = e[i][0];
right[4] = p[i][0];
right[5] = c[i][0];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][0];
left[1] = u[i][0];
left[2] = v[i][0];
left[3] = e[i][0];
left[4] = p[i][0];
left[5] = c[i][0];
right[0] = rho[i][1];
right[1] = u[i][1];
right[2] = v[i][1];
right[3] = e[i][1];
right[4] = p[i][1];
right[5] = c[i][1];
riemann_solver_ydir(left, right, gjph);
dG1[i][0] = (gjph[0] - gjmh[0]) / dy;
dG2[i][0] = (gjph[1] - gjmh[1]) / dy;
dG3[i][0] = (gjph[2] - gjmh[2]) / dy;
dG4[i][0] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_right_bc()
{
float left[6], right[6];
float gjmh[4], gjph[4];
for (int i = 0; i < imax; i++)
{
left[0] = rho[i][jmax - 2];
left[1] = u[i][jmax - 2];
left[2] = v[i][jmax - 2];
left[3] = e[i][jmax - 2];
left[4] = p[i][jmax - 2];
left[5] = c[i][jmax - 2];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][jmax - 1];
left[1] = u[i][jmax - 1];
left[2] = v[i][jmax - 1];
left[3] = e[i][jmax - 1];
left[4] = p[i][jmax - 1];
left[5] = c[i][jmax - 1];
right[0] = rho[i][jmax - 1];
right[1] = u[i][jmax - 1];
right[2] = -v[i][jmax - 1];
right[3] = e[i][jmax - 1];
right[4] = p[i][jmax - 1];
right[5] = c[i][jmax - 1];
riemann_solver_ydir(left, right, gjph);
dG1[i][jmax - 1] = (gjph[0] - gjmh[0]) / dy;
dG2[i][jmax - 1] = (gjph[1] - gjmh[1]) / dy;
dG3[i][jmax - 1] = (gjph[2] - gjmh[2]) / dy;
dG4[i][jmax - 1] = (gjph[3] - gjmh[3]) / dy;
}
return;
}
void calc_flux_ydir_internal()
{
float left[6], right[6];
float gjmh[4], gjph[4];
float32x4_t vgjmh[4], vgjph[4];
float coef = 1.0 / dy;
int jmax4 = jmax / 4 * 4;
for (int i = 0; i < imax; i++)
{
for (int j = 1; j < jmax4 + 1; j += 4)
{
riemann_solver_ydir_direct(i, j - 1, j, vgjmh);
riemann_solver_ydir_direct(i, j, j + 1, vgjph);
vst1q_f32(&(dG1[i][j]), vmulq_n_f32(vsubq_f32(vgjph[0], vgjmh[0]), coef));
vst1q_f32(&(dG2[i][j]), vmulq_n_f32(vsubq_f32(vgjph[1], vgjmh[1]), coef));
vst1q_f32(&(dG3[i][j]), vmulq_n_f32(vsubq_f32(vgjph[2], vgjmh[2]), coef));
vst1q_f32(&(dG4[i][j]), vmulq_n_f32(vsubq_f32(vgjph[3], vgjmh[3]), coef));
}
for (int j = jmax4 + 1; j < jmax - 1; j++)
{
left[0] = rho[i][j - 1];
left[1] = u[i][j - 1];
left[2] = v[i][j - 1];
left[3] = e[i][j - 1];
left[4] = p[i][j - 1];
left[5] = c[i][j - 1];
right[0] = rho[i][j];
right[1] = u[i][j];
right[2] = v[i][j];
right[3] = e[i][j];
right[4] = p[i][j];
right[5] = c[i][j];
riemann_solver_ydir(left, right, gjmh);
left[0] = rho[i][j];
left[1] = u[i][j];
left[2] = v[i][j];
left[3] = e[i][j];
left[4] = p[i][j];
left[5] = c[i][j];
right[0] = rho[i][j + 1];
right[1] = u[i][j + 1];
right[2] = v[i][j + 1];
right[3] = e[i][j + 1];
right[4] = p[i][j + 1];
right[5] = c[i][j + 1];
riemann_solver_ydir(left, right, gjph);
dG1[i][j] = (gjph[0] - gjmh[0]) / dy;
dG2[i][j] = (gjph[1] - gjmh[1]) / dy;
dG3[i][j] = (gjph[2] - gjmh[2]) / dy;
dG4[i][j] = (gjph[3] - gjmh[3]) / dy;
}
}
return;
}
void calc_flux()
{
//----------X direction --------------
calc_flux_xdir_left_bc();
calc_flux_xdir_right_bc();
calc_flux_xdir_internal();
//----------Y direction --------------
calc_flux_ydir_left_bc();
calc_flux_ydir_right_bc();
calc_flux_ydir_internal();
return;
}
void Update(void)
{
int jmax4 = jmax / 4 * 4;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax4; j += 4)
{
float32x4_t vQ1 = vld1q_f32(&Q1[i][j]);
float32x4_t vQ2 = vld1q_f32(&Q2[i][j]);
float32x4_t vQ3 = vld1q_f32(&Q3[i][j]);
float32x4_t vQ4 = vld1q_f32(&Q4[i][j]);
vQ1 = vfmaq_n_f32(vQ1, vaddq_f32(vld1q_f32(&dF1[i][j]), vld1q_f32(&dG1[i][j])), -dt);
vQ2 = vfmaq_n_f32(vQ2, vaddq_f32(vld1q_f32(&dF2[i][j]), vld1q_f32(&dG2[i][j])), -dt);
vQ3 = vfmaq_n_f32(vQ3, vaddq_f32(vld1q_f32(&dF3[i][j]), vld1q_f32(&dG3[i][j])), -dt);
vQ4 = vfmaq_n_f32(vQ4, vaddq_f32(vld1q_f32(&dF4[i][j]), vld1q_f32(&dG4[i][j])), -dt);
vst1q_f32(&Q1[i][j], vQ1);
vst1q_f32(&Q2[i][j], vQ2);
vst1q_f32(&Q3[i][j], vQ3);
vst1q_f32(&Q4[i][j], vQ4);
}
for (int j = jmax4; j < jmax; j++)
{
Q1[i][j] -= dt * (dF1[i][j] + dG1[i][j]);
Q2[i][j] -= dt * (dF2[i][j] + dG2[i][j]);
Q3[i][j] -= dt * (dF3[i][j] + dG3[i][j]);
Q4[i][j] -= dt * (dF4[i][j] + dG4[i][j]);
}
}
}
void Evolve(void)
{
for (int n = 0; n < nmax; n++)
{
cons_2_prim();
calc_dt();
cout << "*****************" << endl;
cout << "n = " << n << "; dt = " << dt << " \n";
calc_min_max();
calc_flux();
Update();
}
cout << "*****************" << endl;
}
void Output_results()
{
float vel;
ofstream f;
f.open("rho_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << rho[i][j] << endl;
}
}
f.close();
f.open("p_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
f << p[i][j] << endl;
}
}
f.close();
f.open("vel_output.csv", ios_base::out);
f << imax * jmax << endl;
for (int i = 0; i < imax; i++)
{
for (int j = 0; j < jmax; j++)
{
vel = sqrt(pow(u[i][j], 2.0) + pow(v[i][j], 2.0));
f << vel << endl;
}
}
f.close();
}
int main()
{
double t_start[10], t_end[10];
t_start[0] = mytimer();
Initialize();
Evolve();
t_end[0] = mytimer();
Output_results();
cout << "Total time:" << t_end[0] - t_start[0] << endl;
return 0;
}- cfd2d.h
#ifndef CFD2D_H_
#define CFD2D_H_
float **Array2D(unsigned int row, unsigned int col);
void cons_2_prim();
void Initialize(void);
void calc_dt();
void calc_min_max();
void calc_flux_xdir_left_bc();
void calc_flux_xdir_right_bc();
void calc_flux_xdir_internal();
void calc_flux_ydir_left_bc();
void calc_flux_ydir_right_bc();
void calc_flux_ydir_internal();
void calc_flux();
void Update(void);
void Evolve(void);
void Output_results();
#endif
- input_file
nmax 500
dx 1.0e-2
dy 1.0e-2
imax 1026
jmax 1026
CFL 0.5
- Makefile
CXX = g++ -O1 -ftree-vectorize
TARGET = eular.opt_simd.exe
SOURCE = $(wildcard *.cpp)
OBJECT = $(patsubst %.cpp, %.o, $(SOURCE))
$(TARGET): $(OBJECT)
$(CXX) $^ -o $@
%.o: %.cpp
$(CXX) -c $< -o $@
.PHONY: clean
clean:
rm -rf $(TARGET) *.o- riemann.cpp
#include "riemann.h"
void riemann_solver_xdir_direct(int il, int ir, int j, float32x4_t vflux[])
{
float32x4_t vrhol = vld1q_f32(&rho[il][j]);
float32x4_t vul = vld1q_f32(&u[il][j]);
float32x4_t vvl = vld1q_f32(&v[il][j]);
float32x4_t vel = vld1q_f32(&e[il][j]);
float32x4_t vpl = vld1q_f32(&p[il][j]);
float32x4_t vcl = vld1q_f32(&c[il][j]);
float32x4_t vrhor = vld1q_f32(&rho[ir][j]);
float32x4_t vur = vld1q_f32(&u[ir][j]);
float32x4_t vvr = vld1q_f32(&v[ir][j]);
float32x4_t ver = vld1q_f32(&e[ir][j]);
float32x4_t vpr = vld1q_f32(&p[ir][j]);
float32x4_t vcr = vld1q_f32(&c[ir][j]);
float32x4_t vwl = vaddq_f32(vabsq_f32(vul), vcl);
float32x4_t vwr = vaddq_f32(vabsq_f32(vur), vcr);
float32x4_t vcwave = vbslq_f32(vcgtq_f32(vwl, vwr), vwl, vwr);
float32x4_t vfl[4], vfr[4], vflr, vfrr;
vfl[0] = vmulq_f32(vrhol, vul);
vfl[1] = vfmaq_f32(vpl, vfl[0], vul);
vfl[2] = vmulq_f32(vfl[0], vvl);
vflr = vaddq_f32(
vaddq_f32(vel,
vmulq_n_f32(vmulq_f32(vul, vul), 0.5)),
vmulq_n_f32(vmulq_f32(vvl, vvl), 0.5)
);
vfl[3] = vfmaq_f32(vmulq_f32(vpl, vul), vfl[0], vflr);
vfr[0] = vmulq_f32(vrhor, vur);
vfr[1] = vfmaq_f32(vpr, vfr[0], vur);
vfr[2] = vmulq_f32(vfr[0], vvr);
vfrr = vaddq_f32(
vaddq_f32(ver,
vmulq_n_f32(vmulq_f32(vur, vur), 0.5)),
vmulq_n_f32(vmulq_f32(vvr, vvr), 0.5)
);
vfr[3] = vfmaq_f32(vmulq_f32(vpr, vur), vfr[0], vfrr);
float32x4_t vQl[4], vQr[4];
vQl[0] = vrhol;
vQl[1] = vmulq_f32(vrhol, vul);
vQl[2] = vmulq_f32(vrhol, vvl);
vQl[3] = vmulq_f32(vQl[0], vflr);
vQr[0] = vrhor;
vQr[1] = vmulq_f32(vrhor, vur);
vQr[2] = vmulq_f32(vrhor, vvr);
vQr[3] = vmulq_f32(vQr[0], vfrr);
for (int i = 0; i < 4; i++)
{
vflux[i] = vsubq_f32(
vmulq_n_f32(vaddq_f32(vfl[i], vfr[i]), 0.5),
vmulq_f32(
vmulq_n_f32(vcwave, 0.5),
vsubq_f32(vQr[i], vQl[i]))
);
}
}
void riemann_solver_xdir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(ul) + cl;
wr = abs(ur) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * ul;
fl[1] = fl[0] * ul + pl * 1.0;
fl[2] = fl[0] * vl;
fl[3] = fl[0] * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * ul;
// right flux
fr[0] = rhor * ur;
fr[1] = fr[0] * ur + pr * 1.0;
fr[2] = fr[0] * vr;
fr[3] = fr[0] * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * ur;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}
void riemann_solver_ydir_direct(int i, int jl, int jr, float32x4_t vflux[])
{
float32x4_t vrhol = vld1q_f32(&rho[i][jl]);
float32x4_t vul = vld1q_f32(&u[i][jl]);
float32x4_t vvl = vld1q_f32(&v[i][jl]);
float32x4_t vel = vld1q_f32(&e[i][jl]);
float32x4_t vpl = vld1q_f32(&p[i][jl]);
float32x4_t vcl = vld1q_f32(&c[i][jl]);
float32x4_t vrhor = vld1q_f32(&rho[i][jr]);
float32x4_t vur = vld1q_f32(&u[i][jr]);
float32x4_t vvr = vld1q_f32(&v[i][jr]);
float32x4_t ver = vld1q_f32(&e[i][jr]);
float32x4_t vpr = vld1q_f32(&p[i][jr]);
float32x4_t vcr = vld1q_f32(&c[i][jr]);
float32x4_t vwl = vaddq_f32(vabsq_f32(vvl), vcl);
float32x4_t vwr = vaddq_f32(vabsq_f32(vvr), vcr);
float32x4_t vcwave = vbslq_f32(vcgtq_f32(vwl, vwr), vwl, vwr);
float32x4_t vfl[4], vfr[4], vflr, vfrr;
vfl[0] = vmulq_f32(vrhol, vvl);
vfl[1] = vmulq_f32(vfl[0], vul);
vfl[2] = vaddq_f32(vmulq_f32(vfl[0], vvl), vpl);
vflr = vaddq_f32(
vaddq_f32(vel,
vmulq_n_f32(vmulq_f32(vul, vul), 0.5)),
vmulq_n_f32(vmulq_f32(vvl, vvl), 0.5)
);
vfl[3] = vfmaq_f32(vmulq_f32(vfl[0], vflr), vpl, vvl);
vfr[0] = vmulq_f32(vrhor, vvr);
vfr[1] = vmulq_f32(vfr[0], vur);
vfr[2] = vfmaq_f32(vpr, vfr[0], vvr);
vfrr = vaddq_f32(
vaddq_f32(ver,
vmulq_n_f32(vmulq_f32(vur, vur), 0.5)),
vmulq_n_f32(vmulq_f32(vvr, vvr), 0.5)
);
vfr[3] = vfmaq_f32(vmulq_f32(vfr[0], vfrr), vpr, vvr);
float32x4_t vQl[4], vQr[4];
vQl[0] = vrhol;
vQl[1] = vmulq_f32(vrhol, vul);
vQl[2] = vmulq_f32(vrhol, vvl);
vQl[3] = vmulq_f32(vrhol, vflr);
vQr[0] = vrhor;
vQr[1] = vmulq_f32(vrhor, vur);
vQr[2] = vmulq_f32(vrhor, vvr);
vQr[3] = vmulq_f32(vrhor, vfrr);
for (int i = 0; i < 4; i++)
{
vflux[i] = vsubq_f32(
vmulq_n_f32(vaddq_f32(vfl[i], vfr[i]), 0.5),
vmulq_f32(
vmulq_n_f32(vcwave, 0.5),
vsubq_f32(vQr[i], vQl[i]))
);
}
}
void riemann_solver_ydir(float left[], float right[], float flux[])
{
float rhol = left[0];
float ul = left[1];
float vl = left[2];
float el = left[3];
float pl = left[4];
float cl = left[5];
float rhor = right[0];
float ur = right[1];
float vr = right[2];
float er = right[3];
float pr = right[4];
float cr = right[5];
float cwave, wl, wr;
wl = abs(vl) + cl;
wr = abs(vr) + cr;
cwave = (wl > wr) ? wl : wr;
float fl[4], fr[4];
// left flux
fl[0] = rhol * vl;
fl[1] = fl[0] * ul;
fl[2] = fl[0] * vl + pl * 1.0;
fl[3] = fl[0] * (el + ul * ul / 2.0 + vl * vl / 2.0) + pl * vl;
// right flux
fr[0] = rhor * vr;
fr[1] = fr[0] * ur;
fr[2] = fr[0] * vr + pr * 1.0;
fr[3] = fr[0] * (er + ur * ur / 2.0 + vr * vr / 2.0) + pr * vr;
float Ql[4], Qr[4];
Ql[0] = rhol;
Ql[1] = rhol * ul;
Ql[2] = rhol * vl;
Ql[3] = rhol * (el + ul * ul / 2.0 + vl * vl / 2.0);
Qr[0] = rhor;
Qr[1] = rhor * ur;
Qr[2] = rhor * vr;
Qr[3] = rhor * (er + ur * ur / 2.0 + vr * vr / 2.0);
for (int i = 0; i < 4; i++)
{
flux[i] = 0.5 * (fl[i] + fr[i]) - 0.5 * cwave * (Qr[i] - Ql[i]);
}
}- riemann.h
#ifndef RIEMANN_H_
#define RIEMANN_H_
#include <iostream>
#include <stdlib.h>
#include <arm_neon.h>
using namespace std;
extern float **Q1, **Q2, **Q3, **Q4, **dF1, **dF2, **dF3, **dF4, **dG1, **dG2, **dG3, **dG4;
extern float dx, dy, dt;
extern float **rho, **u, **v, **e, **p, **c;
extern int nmax, imax, jmax;
extern float cfl;
void riemann_solver_xdir(float left[], float right[], float flux[]);
void riemann_solver_ydir(float left[], float right[], float flux[]);
void riemann_solver_xdir_direct(int il, int ir, int j, float32x4_t vflux[]);
void riemann_solver_ydir_direct(int i, int jl, int jr, float32x4_t vflux[]);
#endif- timer.cpp
#include "timer.h"
#include <cstdlib>
#include <sys/time.h>
#include <sys/resource.h>
double mytimer(void)
{
struct timeval tp;
static long start = 0, startu;
if (!start)
{
gettimeofday(&tp, NULL);
start = tp.tv_sec;
startu = tp.tv_usec;
return 0.0;
}
gettimeofday(&tp, NULL);
return ((double) (tp.tv_sec - start)) + (tp.tv_usec - startu) / 1000000.0;
}- timer.h
#ifndef TIMER_H_
#define TIMER_H_
double mytimer(void);
#endif- yhrun.sh
#!/bin/bash
yhrun -p thcp1 -N 1 -n 1 eular.opt_simd.exe运行结果如下：加速比为 3.19，效果显著。
2233
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
