So I actually changed the implementation of the reaction for the model. It is no longer easy to see the problem (of heat being lost between the divided sections of the lattice) in the images, but I can still see it numerically since heat is being lost in the system (even though I’m using adiabatic boundaries in the x and periodic in the y and z right now). I tried changing the old reaction to bulk, and splitting up the reading and writing into two processors, but that did not work. The same problem is happening with the new reaction implementation. Below is the relevant code that you need to see to understand the reaction process, commented for clarity. Thanks for all your help, hope this gives you some insight on the problem.
// neighbor data has the central point (x, y, z) and a vector of any of the 6 relevant directions that contain interesting neighbors
// so if x = 5, y = 5, z = 5, and adj = {0, 2, 3}, we know that when processing the point (5, 5, 5) to process the reaction in the 0th, 2nd, and 3rd directions (dx [i], dy [i], dz [i] store the delta_x, delta_y, and delta_z for the ith direction, used below)
struct neighbor
{
plint x, y, z;
vector <plint> adj;
neighbor (plint x_, plint y_, plint z_, vector <plint> adj_)
: x (x_), y (y_), z (z_), adj (adj_)
{ }
};
vector <neighbor> neighbor_list; // list of points and their neighbors
vector <vector <T> > reaction_temps; // for each central point, the temperature I need to read in from the lattice
// read data from lattice that is needed for heat-transfer reaction
struct reaction_read : public BoxProcessingFunctional3D_LL <T, AD, T, AD>
{
reaction_read ()
{ }
virtual void process (Box3D domain, BlockLattice3D <T, AD>& lattice1, BlockLattice3D <T, AD>& lattice2)
{
Dot3D abs = lattice1.getLocation ();
reaction_temps.clear (); // reset data storage vector
for (plint i = 0; i < (plint) neighbor_list.size (); i++) // for all central points
{
vector <T> temp; // list of temperatures for this central point
plint fx = neighbor_list [i].x, fy = neighbor_list [i].y, fz = neighbor_list [i].z; // central point (fx, fy, fz)
temp.push_back (lattice2.get (fx, fy, fz).computeDensity ()); // the temperature of the central point
for (plint j = 0; j < (plint) neighbor_list [i].adj.size (); j++) // for each neighbor
{
plint px = fx + dx [neighbor_list [i].adj [j]], // neighbor (px, py, pz)
py = (fy + abs.y + dy [neighbor_list [i].adj [j]] + ny) % ny - abs.y, // y and z are periodic, so convert to global
pz = (fz + abs.z + dz [neighbor_list [i].adj [j]] + nz) % nz - abs.z; // and then convert back to local coordinates
temp.push_back (lattice1.get (px, py, pz).computeDensity ()); // temperature of neighbor
}
reaction_temps.push_back (temp); // add list for this central point to main data storage vector
}
}
virtual reaction_read* clone () const
{
return new reaction_read (*this);
}
virtual void getModificationPattern (vector <bool>& isWritten) const
{
isWritten [0] = false;
isWritten [1] = false;
}
virtual BlockDomain::DomainT appliesTo () const
{
return BlockDomain::bulk;
}
};
// write the lattices, updating temperatures for heat-transfer reaction
struct reaction_write : public BoxProcessingFunctional3D_LL <T, AD, T, AD>
{
reaction_write ()
{ }
virtual void process (Box3D domain, BlockLattice3D <T, AD>& lattice1, BlockLattice3D <T, AD>& lattice2)
{
Dot3D abs = lattice1.getLocation ();
Array <T, 3> u (0.0, 0.0, 0.0);
for (plint i = 0; i < (plint) neighbor_list.size (); i++) // for all central points
{
plint fx = neighbor_list [i].x, fy = neighbor_list [i].y, fz = neighbor_list [i].z; // central point (fx, fy, fz)
T transfer_total = 0.0; // keep track of all heat that will be added/subtracted to the central point
for (plint j = 0; j < (plint) neighbor_list [i].adj.size (); j++) // for all neighbors
{
plint px = fx + dx [neighbor_list [i].adj [j]], // neighbor (px, py, pz)
py = (fy + abs.y + dy [neighbor_list [i].adj [j]] + ny) % ny - abs.y, // y and z are periodic, so convert to global
pz = (fz + abs.z + dz [neighbor_list [i].adj [j]] + nz) % nz - abs.z; // and then convert back to local coordinates
// H = heat transfer coefficient, r..[i][0] = central point temperature, r..[i][j + 1] = neighbor temperature
T transfer_local = H * (reaction_temps [i][0] - reaction_temps [i][j + 1]); // transfer based on temperature gradient
// update neighbor temperature
lattice1.get (px, py, pz).computeVelocity (u);
iniCellAtEquilibrium (lattice1.get (px, py, pz), reaction_temps [i][j + 1] + transfer_local, u);
// update total heat transfer to central point
transfer_total += transfer_local;
}
// update central point temperature
lattice2.get (fx, fy, fz).computeVelocity (u);
iniCellAtEquilibrium (lattice2.get (fx, fy, fz), reaction_temps [i][0] - transfer_total, u);
}
}
virtual reaction_write* clone () const
{
return new reaction_write (*this);
}
virtual void getModificationPattern (vector <bool>& isWritten) const
{
isWritten [0] = true;
isWritten [1] = true;
}
virtual BlockDomain::DomainT appliesTo () const
{
return BlockDomain::bulk;
}
};
