Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 89 additions & 1 deletion applications/solvers/dfLowMachFoam/EEqn.H
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,23 @@

#ifdef CPUSolver_
start1 = std::clock();
//debug
// {
// const fvPatchScalarField& hew = he.boundaryField()[5];
// const basicThermo& bThermo = basicThermo::lookupThermo(hew);
// const scalarField& pw = bThermo.p().boundaryField()[5];
// fvPatchScalarField& Tw =
// const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[5]);
// scalarField& Tw_v = Tw;

// Tw.evaluate();

// Info << "internal field" <<bThermo.he(pw, Tw, mesh.boundary()[5].faceCells()) << endl;
// Info << "boundary field" <<bThermo.he(pw, Tw, 5) << endl;
// Info << "calculated grad" << mesh.boundary()[5].deltaCoeffs() * (bThermo.he(pw, Tw, 5) - bThermo.he(pw, Tw, mesh.boundary()[5].faceCells())) << endl;
// }


fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + fvm::div(phi, he)
Expand All @@ -13,6 +30,9 @@
==
fvc::div(hDiffCorrFlux)
);
printf("EEqn = %.15lf\n", EEqn.boundaryCoeffs()[0][64]);
printf("EEqn = %.15lf\n", EEqn.boundaryCoeffs()[1][64]);
printf("EEqn = %.15lf\n", EEqn.boundaryCoeffs()[5][1]);
end1 = std::clock();
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand Down Expand Up @@ -46,12 +66,44 @@
end2 = std::clock();
int eeqn_offset = 0;
int patchNum = 0;

//debug
// {
// const fvPatchScalarField& hew = he.boundaryField()[5];
// const basicThermo& bThermo = basicThermo::lookupThermo(hew);
// const scalarField& pw = bThermo.p().boundaryField()[5];
// fvPatchScalarField& Tw =
// const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[5]);
// scalarField& Tw_v = Tw;

// Tw.evaluate();

// Info << "internal field" <<bThermo.he(pw, Tw, mesh.boundary()[5].faceCells()) << endl;
// Info << "boundary field" <<bThermo.he(pw, Tw, 5) << endl;
// Info << "calculated grad" <<
// mesh.boundary()[5].deltaCoeffs() * (bThermo.he(pw, Tw, 5) - bThermo.he(pw, Tw, mesh.boundary()[5].faceCells())) << endl;
// }

forAll(he.boundaryField(), patchi)
{
patchNum++;
const fvsPatchScalarField& pw = mesh.surfaceInterpolation::weights().boundaryField()[patchi];
int patchSize = pw.size();

// construct gradient manually
const fvPatchScalarField& hew = he.boundaryField()[patchi];
const basicThermo& bThermo = basicThermo::lookupThermo(hew);
const scalarField& ppw = bThermo.p().boundaryField()[patchi];
fvPatchScalarField& Tw =
const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[patchi]);
scalarField& Tw_v = Tw;

Tw.evaluate();
const scalarField& patchDeltaCoeff = mesh.boundary()[patchi].deltaCoeffs();
const scalarField heInternal = bThermo.he(ppw, Tw, patchi)();
const scalarField heBoundary = bThermo.he(ppw, Tw, mesh.boundary()[patchi].faceCells())();
const scalarField patchGradMau = patchDeltaCoeff * (heInternal - heBoundary);

const scalarField& patchK = K.boundaryField()[patchi];
// const scalarField& patchAlphaEff = alphaEff.boundaryField()[patchi];
const scalarField& patchGrad = he.boundaryField()[patchi].gradientBoundaryCoeffs(); // gradient_
Expand All @@ -60,6 +112,42 @@
memcpy(boundary_gradient + eeqn_offset, &patchGrad[0], patchSize*sizeof(double));

eeqn_offset += patchSize;

// debug
// const fvsPatchScalarField& patchFlux = phi.boundaryField()[patchi];
// // Field<scalar> valueInternalCoeffs = patchFlux*he.boundaryField()[patchi].valueInternalCoeffs(pw);
// // Field<scalar> valueBoundaryCoeffs = -patchFlux*he.boundaryField()[patchi].valueBoundaryCoeffs(pw);
// Field<scalar> valueInternalCoeffs = patchFlux*he.boundaryField()[patchi].valueInternalCoeffs(pw);
// Field<scalar> valueBoundaryCoeffs = -patchFlux*he.boundaryField()[patchi].valueBoundaryCoeffs(pw);
// Field<scalar> gradientInternalCoeffs = he.boundaryField()[patchi].gradientInternalCoeffs();
// Field<scalar> gradientBoundaryCoeffs = he.boundaryField()[patchi].gradientBoundaryCoeffs();
// labelUList sub_boundary = mesh.boundary()[patchi].faceCells();

// calculate grad manually
// const fvPatchScalarField& hew = he.boundaryField()[patchi];
// const basicThermo& bThermo = basicThermo::lookupThermo(hew);
// const scalarField& ppw = bThermo.p().boundaryField()[patchi];
// fvPatchScalarField& Tw =
// const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[patchi]);
// scalarField& Tw_v = Tw;

// Tw.evaluate();
// const scalarField& patchDeltaCoeff = mesh.boundary()[patchi].deltaCoeffs();
// const scalarField heInternal = bThermo.he(ppw, Tw, patchi)();
// const scalarField heBoundary = bThermo.he(ppw, Tw, mesh.boundary()[patchi].faceCells())();
// const scalarField patchGradMau = patchDeltaCoeff * (heInternal - heBoundary);

// forAll(sub_boundary, i){
// if (sub_boundary[i] == 1)
// {
// printf("\npatchFlux = %.10lf\n", patchFlux[i]);
// printf("valueBoundaryCoeffs = %.10lf\n", he.boundaryField()[patchi].valueBoundaryCoeffs(pw)()[i]);
// printf("valueBoundaryCoeffs_CPU = %.10lf\n", valueBoundaryCoeffs[i]);
// printf("patchGrad = %.10lf\n", patchGrad[i]);
// printf("patchGradMau = %.10lf\n", patchGradMau[i]);
// printf("patch = %d, cellID = %d\n", patchi, i);
// }
// }
}
end1 = std::clock();
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand Down Expand Up @@ -95,7 +183,7 @@
time_monitor_EEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);

// check value of mtxAssembly, no time monitor
// EEqn_GPU.checkValue(false);
// EEqn_GPU.checkValue(true);

start1 = std::clock();
EEqn_GPU.solve();
Expand Down
46 changes: 40 additions & 6 deletions applications/solvers/dfLowMachFoam/YEqn.H
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ forAll(Y, i)
{
sumYDiffError += chemistry->rhoD(i)*fvc::grad(Y[i]);
}
// Info << "sumYDiffError\n" << sumYDiffError << endl;
const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
start1 = std::clock();
time_monitor_YEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand All @@ -43,29 +44,35 @@ time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
{
cudaMallocHost(&boundary_Y[i], num_boundary_faces*sizeof(double));
}
// const volScalarField& haii = chemistry->hai(i);
// const volScalarField& rhoDi = chemistry->rhoD(i);
const volScalarField& haii = chemistry->hai(i);
const volScalarField& rhoDi = chemistry->rhoD(i);
// hai[i] = &haii[0];
// rhoD[i] = &rhoDi[0];
rhoD[i] = &rhoDi[0];
// cudaMallocHost(&boundary_hai[i], num_boundary_faces*sizeof(double));
// cudaMallocHost(&boundary_rhoD[i], num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_rhoD[i], num_boundary_faces*sizeof(double));
int offset = 0;
forAll(Yi.boundaryField(), patchi)
{
const scalarField& patchYi = Yi.boundaryField()[patchi];
// const scalarField& patchHaii = haii.boundaryField()[patchi];
// const scalarField& patchRhoDi = rhoDi.boundaryField()[patchi];
const scalarField& patchRhoDi = rhoDi.boundaryField()[patchi];
int patchSize = patchYi.size();

if (updateBoundaryFields)
{
memcpy(boundary_Y[i] + offset, &patchYi[0], patchSize*sizeof(double));
}
// memcpy(boundary_hai[i] + offset, &patchHaii[0], patchSize*sizeof(double));
// memcpy(boundary_rhoD[i] + offset, &patchRhoDi[0], patchSize*sizeof(double));
memcpy(boundary_rhoD[i] + offset, &patchRhoDi[0], patchSize*sizeof(double));
offset += patchSize;
}
// if (i == 5)
// {
// Info << "rhoD_CPU" << rhoDi << endl;
// }

}
// Info << "rhoD from nuEff\n" << nuEff * rho / 0.7 << endl;
updateBoundaryFields = false;
volScalarField mut_sct = turbulence->mut().ref()/Sct;
double *boundary_mutsct = nullptr;
Expand All @@ -77,6 +84,17 @@ time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
int patchSize = patchMut_sct.size();
memcpy(boundary_mutsct + offset, &patchMut_sct[0], patchSize*sizeof(double));
offset += patchSize;

// debug
// const fvsPatchScalarField& pw = mesh.surfaceInterpolation::weights().boundaryField()[patchi];
// Field<scalar> valueInternalCoeffs = Y[5].boundaryField()[patchi].valueInternalCoeffs(pw);
// Field<scalar> valueBoundaryCoeffs = Y[5].boundaryField()[patchi].valueBoundaryCoeffs(pw);
// Field<scalar> gradientInternalCoeffs = Y[5].boundaryField()[patchi].gradientInternalCoeffs();
// Field<scalar> gradientBoundaryCoeffs = Y[5].boundaryField()[patchi].gradientBoundaryCoeffs();
// Info << "valueInternalCoeffs\n" << valueInternalCoeffs << endl;
// Info << "valueBoundaryCoeffs\n" << valueBoundaryCoeffs << endl;
// Info << "gradientInternalCoeffs\n" << gradientInternalCoeffs << endl;
// Info << "gradientBoundaryCoeffs\n" << gradientBoundaryCoeffs << endl;
}
end1 = std::clock();
time_monitor_YEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand All @@ -93,6 +111,7 @@ time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
YEqn_GPU.fvm_div_phi();
YEqn_GPU.fvm_div_phiUc();
YEqn_GPU.sync();
// YEqn_GPU.checkValue(true, "of_output_H2.txt");
end1 = std::clock();
time_monitor_YEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
time_monitor_YEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand Down Expand Up @@ -161,6 +180,21 @@ if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
: (fvm::laplacian(DEff(), Yi) + combustion->R(Yi))
)
);
// debug
if (i == 0)
{
printf("b = %lf\n", YiEqn.source()[1]);
forAll(mesh.boundary(), patchi){
labelUList sub_boundary = mesh.boundary()[patchi].faceCells();
forAll(sub_boundary, i){
if (sub_boundary[i] == 0){
printf("patchID = %d, faceID = %d\n", patchi, i);
printf("YiEqn = %.15lf\n", YiEqn.boundaryCoeffs()[patchi][i]);
}
}
}
}

end1 = std::clock();
time_monitor_YEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
// YiEqn.relax();
Expand Down
45 changes: 22 additions & 23 deletions applications/solvers/dfLowMachFoam/createdfSolver.H
Original file line number Diff line number Diff line change
Expand Up @@ -36,31 +36,30 @@ int num_boundary_cells;
string settingPath;
settingPath = CanteraTorchProperties.subDict("AmgxSettings").lookupOrDefault("UEqnSettingPath", string(""));

dfMatrixDataBase dfDataBase(num_surfaces, num_cells, num_boundary_faces, Y.size(), num_boundary_cells, &neighbour[0], &owner[0], &mesh.V()[0], &mesh.surfaceInterpolation::weights()[0],
&mesh.Sf()[0][0], &mesh.magSf()[0], &mesh.nonOrthDeltaCoeffs()[0], boundary_face_vector_init, boundary_face_init, boundary_deltaCoeffs_init, boundaryCellIndex, patchTypes);

#ifdef GPUSolver_
dfRhoEqn rhoEqn_GPU(dfDataBase);
dfUEqn UEqn_GPU(dfDataBase, "dDDI", settingPath);
dfYEqn YEqn_GPU(dfDataBase, "dDDI", settingPath, inertIndex);
dfEEqn EEqn_GPU(dfDataBase, "dDDI", settingPath);
dfMatrixDataBase dfDataBase(num_surfaces, num_cells, num_boundary_faces, Y.size(), num_boundary_cells, &neighbour[0], &owner[0], &mesh.V()[0], &mesh.surfaceInterpolation::weights()[0],
&mesh.Sf()[0][0], &mesh.magSf()[0], &mesh.nonOrthDeltaCoeffs()[0], boundary_face_vector_init, boundary_face_init, boundary_deltaCoeffs_init, boundaryCellIndex, patchTypes);
dfRhoEqn rhoEqn_GPU(dfDataBase);
dfUEqn UEqn_GPU(dfDataBase, "dDDI", settingPath);
dfYEqn YEqn_GPU(dfDataBase, "dDDI", settingPath, inertIndex);
dfEEqn EEqn_GPU(dfDataBase, "dDDI", settingPath);

double *ueqn_internalCoeffs_init, *ueqn_boundaryCoeffs_init, *boundary_pressure_init, *boundary_velocity_init,
*boundary_nuEff_init, *boundary_rho_init, *ueqn_laplac_internalCoeffs_init, *ueqn_laplac_boundaryCoeffs_init, *boundary_phi_init;
cudaMallocHost(&ueqn_internalCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_boundaryCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_laplac_internalCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_laplac_boundaryCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_velocity_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_pressure_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_nuEff_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_rho_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_phi_init, num_boundary_faces*sizeof(double));
double *ueqn_internalCoeffs_init, *ueqn_boundaryCoeffs_init, *boundary_pressure_init, *boundary_velocity_init,
*boundary_nuEff_init, *boundary_rho_init, *ueqn_laplac_internalCoeffs_init, *ueqn_laplac_boundaryCoeffs_init, *boundary_phi_init;
cudaMallocHost(&ueqn_internalCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_boundaryCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_laplac_internalCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&ueqn_laplac_boundaryCoeffs_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_velocity_init, 3*num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_pressure_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_nuEff_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_rho_init, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_phi_init, num_boundary_faces*sizeof(double));

double *boundary_alphaEff, *boundary_K, *boundary_gradient;
cudaMallocHost(&boundary_K, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_alphaEff, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_gradient, num_boundary_faces * sizeof(double));
double *boundary_alphaEff, *boundary_K, *boundary_gradient;
cudaMallocHost(&boundary_K, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_alphaEff, num_boundary_faces*sizeof(double));
cudaMallocHost(&boundary_gradient, num_boundary_faces * sizeof(double));

bool updateBoundaryFields = true; // make sure that the boundary fields do H2D copy at 1st timestep
bool updateBoundaryFields = true; // make sure that the boundary fields do H2D copy at 1st timestep
#endif
10 changes: 6 additions & 4 deletions applications/solvers/dfLowMachFoam/dfLowMachFoam.C
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,8 @@ Description
#include <thread>
#include "upwind.H"

#define GPUSolver_
// #define CPUSolver_
// #define GPUSolver_
#define CPUSolver_

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand Down Expand Up @@ -285,13 +285,15 @@ int main(int argc, char *argv[])
Info<< "other Time = " << time_monitor_other << " s" << endl;
Info<< "rho Equations = " << time_monitor_rho << " s" << endl;
Info<< "U Equations = " << time_monitor_U << " s" << endl;
Info<< "Y Equations = " << time_monitor_Y << " s" << endl;
Info<< "Y Equations = " << time_monitor_Y - time_monitor_chem << " s" << endl;
Info<< "E Equations = " << time_monitor_E << " s" << endl;
Info<< "percentage of rho/U/Y/E = " << (time_monitor_E + time_monitor_Y + time_monitor_U + time_monitor_rho) / loop_time * 100 << " %" << endl;
Info<< "p Equations = " << time_monitor_p << " s" << endl;
Info<< "chemistry correctThermo = " << time_monitor_chemistry_correctThermo << " s" << endl;
Info<< "turbulence correct = " << time_monitor_turbulence_correct << " s" << endl;
Info<< "combustion correct(in Y) = " << time_monitor_chem << " s" << endl;
Info<< "percentage of chemistry = " << time_monitor_chem / loop_time * 100 << " %" << endl;
Info<< "percentage of rho/U/Y/E = " << (time_monitor_E + time_monitor_Y + time_monitor_U + time_monitor_rho - time_monitor_chem) / loop_time * 100 << " %" << endl;


Info<< "========Time details of each equation======="<< endl;

Expand Down
12 changes: 6 additions & 6 deletions src_gpu/dfEEqn.cu
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,7 @@ __global__ void eeqn_fvm_laplacian_uncorrected_boundary(int num_boundary_cells,
double gamma = boundary_alphaEff[i];
double gamma_magsf = gamma * boundary_magsf[i];
internal_coeffs += gamma_magsf * gradient_internal_coeffs[i];
boundary_coeffs += gamma_magsf * gradient_boundary_coeffs[i];
boundary_coeffs -= gamma_magsf * gradient_boundary_coeffs[i];
}

A_csr_output[csr_index] = A_csr_input[csr_index] + internal_coeffs * sign;
Expand Down Expand Up @@ -631,8 +631,8 @@ void dfEEqn::add_to_source()

void dfEEqn::checkValue(bool print)
{
checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, csr_value_vec_bytes, cudaMemcpyDeviceToHost, stream));
checkCudaErrors(cudaMemcpyAsync(h_b, d_b, cell_vec_bytes, cudaMemcpyDeviceToHost, stream));
checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, (num_faces + num_cells) * sizeof(double), cudaMemcpyDeviceToHost, stream));
checkCudaErrors(cudaMemcpyAsync(h_b, d_b, num_cells * sizeof(double), cudaMemcpyDeviceToHost, stream));

// Synchronize stream
checkCudaErrors(cudaStreamSynchronize(stream));
Expand All @@ -644,7 +644,7 @@ void dfEEqn::checkValue(bool print)
fprintf(stderr, "h_b[%d]: %.16lf\n", i, h_b[i]);
}

char *input_file = "of_output.txt";
char *input_file = "of_output_E.txt";
FILE *fp = fopen(input_file, "rb+");
if (fp == NULL)
{
Expand Down Expand Up @@ -672,9 +672,9 @@ void dfEEqn::checkValue(bool print)
if (print)
{
for (int i = 0; i < (num_faces + num_cells); i++)
fprintf(stderr, "h_A_of_vec_1mtx[%d]: %.16lf\n", i, h_A_of_vec_1mtx[i]);
printf("h_A_of_vec_1mtx[%d]: %.16lf\n", i, h_A_of_vec_1mtx[i]);
for (int i = 0; i < num_cells; i++)
fprintf(stderr, "h_b_of_vec[%d]: %.16lf\n", i, h_b_of_vec[i]);
printf("h_b_of_vec[%d]: %.16lf\n", i, h_b_of_vec[i]);
}

// check
Expand Down
3 changes: 2 additions & 1 deletion src_gpu/dfMatrixDataBase.H
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ inline void checkVectorEqual(int count, const double* basevec, double* vec, doub
enum boundaryConditions{
zeroGradient,
fixedValue,
coupled
coupled,
empty
};

void constructBoundarySelector(std::vector<int>& patchTypeSelector, const std::string& patchTypeStr, const int patchSize);
Expand Down
9 changes: 8 additions & 1 deletion src_gpu/dfMatrixDataBase.cu
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ void constructBoundarySelector(std::vector<int>& patchTypeSelector, const std::s
static std::map<std::string, boundaryConditions> BCMap = {
{"zeroGradient", zeroGradient},
{"fixedValue", fixedValue},
{"empty", empty},
{"coupled", coupled}
};
auto iter = BCMap.find(patchTypeStr);
Expand All @@ -31,11 +32,17 @@ void constructBoundarySelector(std::vector<int>& patchTypeSelector, const std::s
patchTypeSelector.insert(patchTypeSelector.end(), tmpSelector.begin(), tmpSelector.end());
break;
}
case coupled:
case empty:
{
tmpSelector.resize(patchSize, 2);
patchTypeSelector.insert(patchTypeSelector.end(), tmpSelector.begin(), tmpSelector.end());
break;
}
case coupled:
{
tmpSelector.resize(patchSize, 3);
patchTypeSelector.insert(patchTypeSelector.end(), tmpSelector.begin(), tmpSelector.end());
break;
}
}
}
Loading