Skip to content

fixed value #305

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
maorz1998 merged 2 commits into from
Jun 28, 2023
Merged

fixed value #305

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bfff1f7
fix bugs in fixedValue BC
maorz1998 Jun 28, 2023
36c4d31
clean version
maorz1998 Jun 28, 2023
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
69 changes: 10 additions & 59 deletions applications/solvers/dfLowMachFoam/EEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,6 @@
==
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 @@ -67,23 +64,6 @@
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++;
Expand All @@ -105,49 +85,17 @@
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_
// const scalarField& patchAlphaEff = alphaEff.boundaryField()[patchi]; // not H2Dcopy when use UnityLewis
// const scalarField& patchGrad = he.boundaryField()[patchi].gradientBoundaryCoeffs(); // gradient_

// const DimensionedField<scalar, volMesh>& patchHa_ = he.boundaryField()[patchi];
// const gradientEnergyFvPatchScalarField patchHa(mesh.boundary()[patchi], patchHa_);
// const scalarField& patchGrad = patchHa.gradient(); // gradient_
memcpy(boundary_K + eeqn_offset, &patchK[0], patchSize*sizeof(double));
// memcpy(boundary_alphaEff + eeqn_offset, &patchAlphaEff[0], patchSize*sizeof(double));
memcpy(boundary_gradient + eeqn_offset, &patchGrad[0], patchSize*sizeof(double));
memcpy(boundary_gradient + eeqn_offset, &patchGradMau[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 All @@ -159,6 +107,8 @@

// prepare data on GPU
start1 = std::clock();
he.oldTime();
K.oldTime();
EEqn_GPU.prepare_data(&he.oldTime()[0], &K[0], &K.oldTime()[0], alphaEff,
&dpdt[0], boundary_K, boundary_alphaEff, boundary_gradient);
EEqn_GPU.sync();
Expand Down Expand Up @@ -194,6 +144,7 @@
start1 = std::clock();
EEqn_GPU.updatePsi(&he[0]);
he.correctBoundaryConditions();
he.write();
end1 = std::clock();
time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
time_monitor_EEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
Expand Down
36 changes: 34 additions & 2 deletions applications/solvers/dfLowMachFoam/UEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@

start1 = std::clock();
UEqn_GPU.initializeTimeStep();
U.oldTime();
UEqn_GPU.fvm_ddt(&U.oldTime()[0][0]);
UEqn_GPU.fvm_div(boundary_pressure_init, boundary_velocity_init, boundary_nuEff_init, boundary_rho_init);
UEqn_GPU.fvc_grad(&p[0]);
Expand Down Expand Up @@ -57,13 +58,44 @@
// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);

// check value
// U.oldTime();
// tmp<fvVectorMatrix> tUEqn
// (
// fvm::ddt(rho, U)
// +
// fvm::div(phi, U)
// +
// turbulence->divDevRhoReff(U)
// == -fvc::grad(p)
// );
// fvVectorMatrix& UEqn = tUEqn.ref();
// printf("b_cpu = %e\n", UEqn.source()[1][1]);
// forAll(U.boundaryField(), patchi){
// labelUList sub_boundary = mesh.boundary()[patchi].faceCells();
// forAll(sub_boundary, i){
// if (sub_boundary[i] == 1){
// printf("b_cpu_bou = %e\n", UEqn.boundaryCoeffs()[patchi][i][1]);
// printf("patchi = %d, i = %d\n", patchi, i);
// }
// }
// }
// const tmp<volTensorField> tgradU(fvc::grad(U));
// const volTensorField& gradU = tgradU();
// Pout << "gradU_of[1]\n" << gradU[1] << nl;
// Pout << "gradU_of[0][64]\n" << gradU.boundaryField()[0][64] << nl;
// Pout << "gradU_of[1][64]\n" << gradU.boundaryField()[1][64] << nl;
// Pout << "gradU_of[5][1]\n" << gradU.boundaryField()[5][1] << nl;
// Pout << "U[1][1]\n" << U[1] << nl;
// Pout << "Ubou[0][64]\n" << U.boundaryField()[0][64] << nl;
// Pout << "Ubou[1][64]\n" << U.boundaryField()[1][64] << nl;
// Pout << "Ubou[5][1]\n" << U.boundaryField()[5][1] << nl;
// if (pimple.momentumPredictor())
// {
// solve(UEqn);

// Info << "U_CPU\n" << U << endl;
// K = 0.5*magSqr(U);
// }
// UEqn_GPU.checkValue(false);
// UEqn_GPU.checkValue(true);
#else
start1 = std::clock();
tmp<fvVectorMatrix> tUEqn
Expand Down
1 change: 1 addition & 0 deletions applications/solvers/dfLowMachFoam/YEqn.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ time_monitor_YEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
for (size_t i = 0; i < Y.size(); ++i)
{
volScalarField& Yi = Y[i];
Yi.oldTime();
Y_old[i] = &Yi.oldTime()[0];
if (updateBoundaryFields)
{
Expand Down
4 changes: 2 additions & 2 deletions applications/solvers/dfLowMachFoam/dfLowMachFoam.C
View file
Open in desktop
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
7 changes: 5 additions & 2 deletions src/dfChemistryModel/dfChemistryModel.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -387,8 +387,10 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()
psi_[celli] = CanteraGas_->meanMolecularWeight()/CanteraGas_->RT();

mu_[celli] = mixture_.CanteraTransport()->viscosity(); // Pa-s


alpha_[celli] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass()); // kg/(m*s)
// alpha_[celli] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass()); // kg/(m*s) mu/alpha = 0.7
alpha_[celli] = mu_[celli] / 0.7;
// thermalConductivity() W/m/K
// cp_mass() J/kg/K

Expand Down Expand Up @@ -457,7 +459,8 @@ void Foam::dfChemistryModel<ThermoType>::correctThermo()

pmu[facei] = mixture_.CanteraTransport()->viscosity();

palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
// palpha[facei] = mixture_.CanteraTransport()->thermalConductivity()/(CanteraGas_->cp_mass());
palpha[facei] = pmu[facei] / 0.7;

if (mixture_.transportModelName() == "UnityLewis")
{
Expand Down
Loading