particle field problem

[liu]

Dear community
Why are there two field(A and B) needed to store the velocity of particles, when in theory, one should be sufficient to directly store the velocity of particles? The following is the code Thank you.
“
if (_forces.begin() != _forces.end()) {
// calculating upwind Gradient
// vel contains velocity information on ADlattice
// velGrad contains upwind vel on ADlattice
if (vel[0]<0.) {
velGrad[0] = vel[0]*(velXp[0]-vel[0]);
velGrad[1] = vel[0]*(velXp[1]-vel[1]);
velGrad[2] = vel[0]*(velXp[2]-vel[2]);
}
else {
velGrad[0] = vel[0]*(vel[0]-velXn[0]);
velGrad[1] = vel[0]*(vel[1]-velXn[1]);
velGrad[2] = vel[0]*(vel[2]-velXn[2]);
}
if (vel[1]<0.) {
velGrad[0] += vel[1]*(velYp[0]-vel[0]);
velGrad[1] += vel[1]*(velYp[1]-vel[1]);
velGrad[2] += vel[1]*(velYp[2]-vel[2]);
}
else {
velGrad[0] += vel[1]*(vel[0]-velYn[0]);
velGrad[1] += vel[1]*(vel[1]-velYn[1]);
velGrad[2] += vel[1]*(vel[2]-velYn[2]);
}
if (vel[2]<0.) {
velGrad[0] += vel[2]*(velZp[0]-vel[0]);
velGrad[1] += vel[2]*(velZp[1]-vel[1]);
velGrad[2] += vel[2]*(velZp[2]-vel[2]);
}
else {
velGrad[0] += vel[2]*(vel[0]-velZn[0]);
velGrad[1] += vel[2]*(vel[1]-velZn[1]);
velGrad[2] += vel[2]*(vel[2]-velZn[2]);
}

for (AdvectionDiffusionForce3D<T, DESCRIPTOR, ADLattice>& f : _forces) { //原始是AdvectionDiffusionForce3D
// writes force in forceValues, vel refers to ADlattice

auto adCell = _partnerLattice->get(x0, y0, z0);
// adCell.template setField<descriptors::TAU_EFF>(tauEff); //新加

f.applyForce(forceValue, &nsCell, &adCell, vel.data(), latticeR);
if (par)
{
_cell.template setField<FIELD_B>(vel);
}
else {
_cell.template setField<FIELD_A>(vel);
}
}

// compute new particle velocity
Vector<T,DESCRIPTOR::d> newVel;
for (int i=0; i < DESCRIPTOR::d; i++) {
newVel[i] = vel[i] + forceValue[i] – velGrad[i];//yuanshi
//newVel[i] = vel[i] – velGrad[i];
}
if (par) {
_cell.template setField<FIELD_A>(newVel);
} else {
_cell.template setField<FIELD_B>(newVel);
}
}
else { // set particle velocity to fluid velocity
nsCell.computeU(velF);
Vector<T,DESCRIPTOR::d> newVel;
for (int i = 0; i < DESCRIPTOR::d; i++) {
newVel[i] = velF[i];
}
if (par) {
_cell.template setField<FIELD_A>(newVel);
} else {
_cell.template setField<FIELD_B>(newVel);
}
}
}
}
}
}
par = !par;”

[mathias]

To compute the gradients, we need information of neighbour nodes and therefore we cannot overwrite the field we need to store once more the entire field!

[liu]

thank you

[Author]

Posts