Skip to content

Commit 7d51928

Browse files
maorz1998maorz1998
authored
fix bugs (#308)
1 parent 00378f7 commit 7d51928

84 files changed

Lines changed: 385 additions & 3054566 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎applications/solvers/dfLowMachFoam/EEqn.H‎

Lines changed: 1 addition & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -1,110 +1,6 @@
11
{
22
volScalarField& he = thermo.he();
3-
#ifdef GPUSolver_
4-
start1 = std::clock();
5-
UEqn_GPU.updatePsi(&U[0][0]);
6-
UEqn_GPU.correctBoundaryConditions();
7-
U.correctBoundaryConditions();
8-
K = 0.5*magSqr(U);
9-
end1 = std::clock();
10-
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
11-
time_monitor_UEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
12-
13-
// prepare data on CPU
14-
start1 = std::clock();
15-
start2 = std::clock();
16-
// const tmp<volScalarField> alphaEff_tmp(thermo.alpha());
17-
// const volScalarField& alphaEff = alphaEff_tmp();
18-
double *alphaEff = nullptr; // tmp
19-
end2 = std::clock();
20-
int eeqn_offset = 0;
21-
int patchNum = 0;
22-
23-
forAll(he.boundaryField(), patchi)
24-
{
25-
patchNum++;
26-
const fvsPatchScalarField& pw = mesh.surfaceInterpolation::weights().boundaryField()[patchi];
27-
int patchSize = pw.size();
28-
29-
// construct gradient manually
30-
const fvPatchScalarField& hew = he.boundaryField()[patchi];
31-
const basicThermo& bThermo = basicThermo::lookupThermo(hew);
32-
const scalarField& ppw = bThermo.p().boundaryField()[patchi];
33-
fvPatchScalarField& Tw =
34-
const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[patchi]);
35-
scalarField& Tw_v = Tw;
36-
37-
Tw.evaluate();
38-
const scalarField& patchDeltaCoeff = mesh.boundary()[patchi].deltaCoeffs();
39-
const scalarField heInternal = bThermo.he(ppw, Tw, patchi)();
40-
const scalarField heBoundary = bThermo.he(ppw, Tw, mesh.boundary()[patchi].faceCells())();
41-
const scalarField patchGradMau = patchDeltaCoeff * (heInternal - heBoundary);
42-
43-
const scalarField& patchK = K.boundaryField()[patchi];
44-
// const scalarField& patchAlphaEff = alphaEff.boundaryField()[patchi]; // not H2Dcopy when use UnityLewis
45-
// const scalarField& patchGrad = he.boundaryField()[patchi].gradientBoundaryCoeffs(); // gradient_
46-
47-
// const DimensionedField<scalar, volMesh>& patchHa_ = he.boundaryField()[patchi];
48-
// const gradientEnergyFvPatchScalarField patchHa(mesh.boundary()[patchi], patchHa_);
49-
// const scalarField& patchGrad = patchHa.gradient(); // gradient_
50-
memcpy(boundary_K + eeqn_offset, &patchK[0], patchSize*sizeof(double));
51-
// memcpy(boundary_alphaEff + eeqn_offset, &patchAlphaEff[0], patchSize*sizeof(double));
52-
memcpy(boundary_gradient + eeqn_offset, &patchGradMau[0], patchSize*sizeof(double));
53-
54-
eeqn_offset += patchSize;
55-
}
56-
end1 = std::clock();
57-
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
58-
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
59-
time_monitor_EEqn_mtxAssembly_CPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
60-
fprintf(stderr, "time_monitor_EEqn_mtxAssembly_CPU_prepare: %lf, build alphaEff time: %lf, patchNum: %d\n",
61-
time_monitor_EEqn_mtxAssembly_CPU_prepare,
62-
double(end2 - start2) / double(CLOCKS_PER_SEC), patchNum);
63-
64-
// prepare data on GPU
65-
start1 = std::clock();
66-
he.oldTime();
67-
K.oldTime();
68-
EEqn_GPU.prepare_data(&he.oldTime()[0], &K[0], &K.oldTime()[0], alphaEff,
69-
&dpdt[0], boundary_K, boundary_alphaEff, boundary_gradient);
70-
EEqn_GPU.sync();
71-
end1 = std::clock();
72-
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
73-
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
74-
time_monitor_EEqn_mtxAssembly_GPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
75-
76-
start1 = std::clock();
77-
EEqn_GPU.initializeTimeStep();
78-
EEqn_GPU.fvm_ddt();
79-
EEqn_GPU.fvm_div();
80-
EEqn_GPU.fvm_laplacian();
81-
EEqn_GPU.fvc_ddt();
82-
EEqn_GPU.fvc_div_phi_scalar();
83-
EEqn_GPU.fvc_div_vector();
84-
EEqn_GPU.add_to_source();
85-
EEqn_GPU.sync();
86-
end1 = std::clock();
87-
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
88-
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
89-
time_monitor_EEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);
90-
91-
// check value of mtxAssembly, no time monitor
92-
// EEqn_GPU.checkValue(true);
93-
94-
start1 = std::clock();
95-
EEqn_GPU.solve();
96-
end1 = std::clock();
97-
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
98-
time_monitor_EEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
99-
100-
start1 = std::clock();
101-
EEqn_GPU.updatePsi(&he[0]);
102-
he.correctBoundaryConditions();
103-
he.write();
104-
end1 = std::clock();
105-
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
106-
time_monitor_EEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
107-
#else
3+
1084
start1 = std::clock();
1095
fvScalarMatrix EEqn
1106
(
@@ -137,5 +33,4 @@
13733
end1 = std::clock();
13834
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
13935
time_monitor_EEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
140-
#endif
14136
}

‎applications/solvers/dfLowMachFoam/EEqn_GPU.H‎

Lines changed: 106 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
{
2+
volScalarField& he = thermo.he();
3+
start1 = std::clock();
4+
UEqn_GPU.updatePsi(&U[0][0]);
5+
UEqn_GPU.correctBoundaryConditions();
6+
U.correctBoundaryConditions();
7+
K = 0.5*magSqr(U);
8+
end1 = std::clock();
9+
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
10+
time_monitor_UEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
11+
12+
// prepare data on CPU
13+
start1 = std::clock();
14+
start2 = std::clock();
15+
// const tmp<volScalarField> alphaEff_tmp(thermo.alpha());
16+
// const volScalarField& alphaEff = alphaEff_tmp();
17+
double *alphaEff = nullptr; // tmp
18+
end2 = std::clock();
19+
int eeqn_offset = 0;
20+
int patchNum = 0;
21+
22+
forAll(he.boundaryField(), patchi)
23+
{
24+
patchNum++;
25+
const fvsPatchScalarField& pw = mesh.surfaceInterpolation::weights().boundaryField()[patchi];
26+
int patchSize = pw.size();
27+
28+
// construct gradient manually
29+
const fvPatchScalarField& hew = he.boundaryField()[patchi];
30+
const basicThermo& bThermo = basicThermo::lookupThermo(hew);
31+
const scalarField& ppw = bThermo.p().boundaryField()[patchi];
32+
fvPatchScalarField& Tw =
33+
const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[patchi]);
34+
scalarField& Tw_v = Tw;
35+
36+
Tw.evaluate();
37+
const scalarField& patchDeltaCoeff = mesh.boundary()[patchi].deltaCoeffs();
38+
const scalarField heInternal = bThermo.he(ppw, Tw, patchi)();
39+
const scalarField heBoundary = bThermo.he(ppw, Tw, mesh.boundary()[patchi].faceCells())();
40+
const scalarField patchGradMau = patchDeltaCoeff * (heInternal - heBoundary);
41+
42+
const scalarField& patchK = K.boundaryField()[patchi];
43+
// const scalarField& patchAlphaEff = alphaEff.boundaryField()[patchi]; // not H2Dcopy when use UnityLewis
44+
// const scalarField& patchGrad = he.boundaryField()[patchi].gradientBoundaryCoeffs(); // gradient_
45+
46+
// const DimensionedField<scalar, volMesh>& patchHa_ = he.boundaryField()[patchi];
47+
// const gradientEnergyFvPatchScalarField patchHa(mesh.boundary()[patchi], patchHa_);
48+
// const scalarField& patchGrad = patchHa.gradient(); // gradient_
49+
memcpy(boundary_K + eeqn_offset, &patchK[0], patchSize*sizeof(double));
50+
// memcpy(boundary_alphaEff + eeqn_offset, &patchAlphaEff[0], patchSize*sizeof(double));
51+
memcpy(boundary_gradient + eeqn_offset, &patchGradMau[0], patchSize*sizeof(double));
52+
53+
eeqn_offset += patchSize;
54+
}
55+
end1 = std::clock();
56+
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
57+
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
58+
time_monitor_EEqn_mtxAssembly_CPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
59+
fprintf(stderr, "time_monitor_EEqn_mtxAssembly_CPU_prepare: %lf, build alphaEff time: %lf, patchNum: %d\n",
60+
time_monitor_EEqn_mtxAssembly_CPU_prepare,
61+
double(end2 - start2) / double(CLOCKS_PER_SEC), patchNum);
62+
63+
// prepare data on GPU
64+
start1 = std::clock();
65+
he.oldTime();
66+
K.oldTime();
67+
EEqn_GPU.prepare_data(&he.oldTime()[0], &K[0], &K.oldTime()[0], alphaEff,
68+
&dpdt[0], boundary_K, boundary_alphaEff, boundary_gradient);
69+
EEqn_GPU.sync();
70+
end1 = std::clock();
71+
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
72+
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
73+
time_monitor_EEqn_mtxAssembly_GPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
74+
75+
start1 = std::clock();
76+
EEqn_GPU.initializeTimeStep();
77+
EEqn_GPU.fvm_ddt();
78+
EEqn_GPU.fvm_div();
79+
EEqn_GPU.fvm_laplacian();
80+
EEqn_GPU.fvc_ddt();
81+
EEqn_GPU.fvc_div_phi_scalar();
82+
EEqn_GPU.fvc_div_vector();
83+
EEqn_GPU.add_to_source();
84+
EEqn_GPU.sync();
85+
end1 = std::clock();
86+
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
87+
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
88+
time_monitor_EEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);
89+
90+
// check value of mtxAssembly, no time monitor
91+
// EEqn_GPU.checkValue(true);
92+
93+
start1 = std::clock();
94+
EEqn_GPU.solve();
95+
end1 = std::clock();
96+
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
97+
time_monitor_EEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
98+
99+
start1 = std::clock();
100+
EEqn_GPU.updatePsi(&he[0]);
101+
he.correctBoundaryConditions();
102+
he.write();
103+
end1 = std::clock();
104+
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
105+
time_monitor_EEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
106+
}

‎applications/solvers/dfLowMachFoam/UEqn.H‎

Lines changed: 23 additions & 131 deletions
Original file line numberDiff line numberDiff line change
@@ -1,132 +1,24 @@
1-
// Solve the Momentum equation
2-
#ifdef GPUSolver_
3-
start1 = std::clock();
4-
int offset = 0;
5-
const tmp<volScalarField> nuEff_tmp(turbulence->nuEff());
6-
const volScalarField& nuEff = nuEff_tmp();
7-
forAll(U.boundaryField(), patchi)
8-
{
9-
const scalarField& patchP = p.boundaryField()[patchi];
10-
const vectorField& patchU = U.boundaryField()[patchi];
11-
const scalarField& patchRho = rho.boundaryField()[patchi];
12-
const scalarField& patchNuEff = nuEff.boundaryField()[patchi];
13-
14-
int patchSize = patchP.size();
15-
16-
// boundary pressure
17-
memcpy(boundary_pressure_init+offset, &patchP[0], patchSize*sizeof(double));
18-
// boundary velocity
19-
memcpy(boundary_velocity_init+3*offset, &patchU[0][0], 3*patchSize*sizeof(double));
20-
// boundary nuEff
21-
memcpy(boundary_nuEff_init+offset, &patchNuEff[0], patchSize*sizeof(double));
22-
// boundary rho
23-
memcpy(boundary_rho_init+offset, &patchRho[0], patchSize*sizeof(double));
24-
offset += patchSize;
25-
}
26-
end1 = std::clock();
27-
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
28-
time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
29-
time_monitor_UEqn_mtxAssembly_CPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
30-
31-
start1 = std::clock();
32-
UEqn_GPU.initializeTimeStep();
33-
U.oldTime();
34-
UEqn_GPU.fvm_ddt(&U.oldTime()[0][0]);
35-
UEqn_GPU.fvm_div(boundary_pressure_init, boundary_velocity_init, boundary_nuEff_init, boundary_rho_init);
36-
UEqn_GPU.fvc_grad(&p[0]);
37-
UEqn_GPU.fvc_grad_vector();
38-
UEqn_GPU.dev2T();
39-
UEqn_GPU.fvc_div_tensor(&nuEff[0]);
40-
UEqn_GPU.fvm_laplacian();
41-
UEqn_GPU.sync();
42-
end1 = std::clock();
43-
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
44-
time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
45-
time_monitor_UEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);
46-
47-
// start2 = std::clock();
48-
// fvVectorMatrix turb_source
49-
// (
50-
// turbulence->divDevRhoReff(U)
51-
// );
52-
// end2 = std::clock();
53-
// time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC);
54-
55-
// UEqn_GPU.add_fvMatrix(&turb_source.lower()[0], &turb_source.diag()[0], &turb_source.upper()[0], &turb_source.source()[0][0]);
56-
// end1 = std::clock();
57-
// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
58-
// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
59-
60-
// check value
61-
// U.oldTime();
62-
// tmp<fvVectorMatrix> tUEqn
63-
// (
64-
// fvm::ddt(rho, U)
65-
// +
66-
// fvm::div(phi, U)
67-
// +
68-
// turbulence->divDevRhoReff(U)
69-
// == -fvc::grad(p)
70-
// );
71-
// fvVectorMatrix& UEqn = tUEqn.ref();
72-
// printf("b_cpu = %e\n", UEqn.source()[1][1]);
73-
// forAll(U.boundaryField(), patchi){
74-
// labelUList sub_boundary = mesh.boundary()[patchi].faceCells();
75-
// forAll(sub_boundary, i){
76-
// if (sub_boundary[i] == 1){
77-
// printf("b_cpu_bou = %e\n", UEqn.boundaryCoeffs()[patchi][i][1]);
78-
// printf("patchi = %d, i = %d\n", patchi, i);
79-
// }
80-
// }
81-
// }
82-
// if (pimple.momentumPredictor())
83-
// {
84-
// solve(UEqn);
85-
// Info << "U_CPU\n" << U << endl;
86-
// K = 0.5*magSqr(U);
87-
// }
88-
// UEqn_GPU.checkValue(true);
89-
#else
90-
start1 = std::clock();
91-
tmp<fvVectorMatrix> tUEqn
92-
(
93-
fvm::ddt(rho, U) + fvm::div(phi, U)
94-
+ turbulence->divDevRhoReff(U)
95-
== -fvc::grad(p)
96-
);
97-
fvVectorMatrix& UEqn = tUEqn.ref();
98-
99-
end1 = std::clock();
100-
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
101-
time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
102-
103-
UEqn.relax();
104-
start1 = std::clock();
105-
if (pimple.momentumPredictor())
106-
{
107-
solve(UEqn);
108-
109-
K = 0.5*magSqr(U);
110-
}
111-
end1 = std::clock();
112-
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
113-
time_monitor_UEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
114-
#endif
115-
116-
// start1 = std::clock();
117-
// // // std::thread t(&dfMatrix::solve, &UEqn_GPU);
118-
// UEqn_GPU.solve();
119-
// end1 = std::clock();
120-
// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
121-
// time_monitor_UEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
122-
123-
// start1 = std::clock();
124-
// // // t.join();
125-
// // UEqn_GPU.updatePsi(&U[0][0]);
126-
// K = 0.5*magSqr(U);
127-
// end1 = std::clock();
128-
// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
129-
// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
130-
// time_monitor_CPU += double(end1 - start1) / double(CLOCKS_PER_SEC);
131-
// // Info << "U_amgx = " << U << endl;
1+
start1 = std::clock();
2+
tmp<fvVectorMatrix> tUEqn
3+
(
4+
fvm::ddt(rho, U) + fvm::div(phi, U)
5+
+ turbulence->divDevRhoReff(U)
6+
);
7+
fvVectorMatrix& UEqn = tUEqn.ref();
8+
9+
end1 = std::clock();
10+
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
11+
time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
12+
13+
UEqn.relax();
14+
start1 = std::clock();
15+
if (pimple.momentumPredictor())
16+
{
17+
solve(UEqn == -fvc::grad(p));
18+
19+
K = 0.5*magSqr(U);
20+
}
21+
end1 = std::clock();
22+
time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
23+
time_monitor_UEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
13224

0 commit comments

Comments
 (0)