Applying Boundary Condition

[Cplusc]

Firstly, I'd like to extend my apologies to anyone viewing this question, as it may not strictly pertain to programming. However, I'm posting it here in the hope that someone might have the answer I'm seeking.

I am utilizing this paper "https://arxiv.org/pdf/0811.4593.pdf" to implement the Zou-He boundary condition in my lattice boltzmann simulation, aiming to enforce a velocity of 1 at the inlet of the complex geometry. The modification involves adjusting the incoming distribution functions. However, upon completion, instead of attaining a density of 1 as expected, I obtain a value of 2 which leads to set the velocity to 0.5. This discrepancy aligns with the mathematical expressions involved. Consider an inlet vel_z=1: initially, the domain initializes with zero velocity across all components. Consequently, the distribution functions on the boundary lattices during the first iteration equal the weight functions, resulting in a summation of 1. After the modification of the incoming distribution functions, certain terms augment the existing distribution functions, logically causing an increase in their summation. This adjustment induces a significant reduction in the velocity of the lattices inside the domain adjacent to the boundary. I explored a lot of other resources including "https://pubs.aip.org/aip/pof/article-abstract/33/7/073304/1060984/Outflow-boundary-condition-of-multiphase?redirectedFrom=fulltext" but I face the same issue. by the way, I am able to apply the Zou He on D2Q9 model successfully. Can anyone help me to get what I am getting wrong in Zou He concept. I share some part of my code as a demo, It's not a programming issue but rather in conceptual question.

void Node::Initialize(std::shared_ptr<VelocitySet> velSet, int dir)
{
    double vel_z = 0.0;
    double rho = 1;
    double du = vel_z * velSet->GetDirection(dir)[2];
    double u_sqr = vel_z * vel_z;
    m_distributions.push_back(velSet->GetWeight(dir) * rho * (1 + 3 * du + 9 / 2 * du * du - 3 / 2 * u_sqr));
    m_newDistributions.push_back(0);
}


#ifdef D3Q19

void Boundary::Apply_ZouHe_Bc(std::shared_ptr<VelocitySet> velSet, Node& lat, std::vector<double>& vel)
{
    double rho = 1.0;
    double tmp = 0.0;
    for (int i = 0; i < _nLatNodes; i++) {

        tmp += lat.m_distributions[i];
    }

    double Nx = (1.0 / 2) * (lat.m_distributions[0] + lat.m_distributions[6] + lat.m_distributions[7]
    - lat.m_distributions[1] - lat.m_distributions[10] - lat.m_distributions[11]) - (1.0 / 3) * (rho * vel[0]);

    double Ny = (1.0 / 2) * (lat.m_distributions[2] + lat.m_distributions[6] + lat.m_distributions[10]
    - lat.m_distributions[3] - lat.m_distributions[7] - lat.m_distributions[11]) - (1.0 / 3) * (rho * vel[1]);


    lat.m_distributions[4] = lat.m_distributions[5] + (1 / 3.0 * rho) * vel[2];

    lat.m_distributions[8] = lat.m_distributions[13] + (1 / 6.0 * rho) * (vel[2] + vel[0]) - Nx;

    lat.m_distributions[12] = lat.m_distributions[9] + (1 / 6.0 * rho) * (vel[2] - vel[0]) + Nx;

    lat.m_distributions[14] = lat.m_distributions[17] + (1 / 6.0 * rho) * (vel[2] + vel[1]) - Ny;

    lat.m_distributions[16] = lat.m_distributions[15] + (1 / 6.0 * rho) * (vel[2] - vel[1]) + Ny;

    tmp = 0.0;
    for (int i = 0; i < _nLatNodes; i++) {

        tmp += lat.m_distributions[i];
    }
}

#endif 

[TheIdeasMan]

Hi

I don't have any domain knowledge for your topic, so forgive me if I am stating something well known.

On line 7 are you sure that: you want integer division; and the precedence of those operations are correct?
It can be amazing how putting extra parentheses to make the precedence explicit can fix things.

What is the type of that expression? Is it double? If so, prefer to write the numbers as double not int, with digits both sides of the decimal point. It saves the compiler doing an implicit conversion from int to double, and takes away the risk of integer division.

EDIT:

Also using a debugger can help with these situations.

EDIT2:

When I am doing code that follows some documentation I put a link to it as a comment in the code, then I put further comments to identify which formula in the docs the code is implementing, making it easy for others to follow along. For example the code on line 7 is which equation?

Good Luck !!

[Cplusc]

@TheIdeasMan Thank you for your answer. though you are right about the integer division, but that is not the issue here as far as du and u_sqr are zeros.

[TheIdeasMan]

@Cplusc

Ok, no worries. Have you made any progress using a debugger?