Handling of boundary nodes for setZeroGradientBoundary
Dear developers,
I am using OpenLB to simulate the transport of a dilute species within a channel of square-shaped cross section, please see the attached images for more details. The domain consists of:
- Inflow: Constant concentration set by "setAdvectionDiffusionTemperatureBoundary", i.e. Dirichlet value of 1.0
- Outflow: Zero normal gradient set by "setZeroGradientBoundary"
- Wall: Constant concentration set by "setAdvectionDiffusionTemperatureBoundary", i.e. Dirichlet value of 0.0
- Initial condition: external velocity profile is imposed.
I have noticed that the way boundary nodes are assigned to "setZeroGradientBoundary" can lead to unbounded values for concentration, e.g., negative values near outflow region. To be more specific, please consider the following two cases as an example.
Case #1: unbounded values for concentration
In this case, the inflow, wall, and outflow boundary materials are assigned according to the image below where corner nodes are assigned to the wall material ID. The "prepareGeometry" function is also provided.
const std::vector<T> origin(3, T(0));
const std::vector<T> extend = {length, width, width};
IndicatorCuboid3D<T> duct(extend, origin);
const int noOfCuboids = singleton::mpi().getSize();
CuboidGeometry3D<T> cuboidGeometry
(
duct, AdvDiffusionConverter.getConversionFactorLength(), noOfCuboids
);
void prepareGeometry
(
XMLreader& configuration, SuperGeometry<T, 3>& superGeometry,
UnitConverter<T, AdvDiffusionDescriptor> const& AdvDiffusionConverter
)
{
OstreamManager clout(std::cout, "prepareGeometry");
//- Setting the material IDs for fluid and boundary cells
//- (1) Changing the material ID from default value 0 to 2
//- (2) Changing the material ID from 2 to 1 for internal cells,
//- so that only the boundary cells have material ID 2
superGeometry.rename(0, 2);
superGeometry.rename(2, 1, {1, 1, 1});
superGeometry.clean();
//- Read duct length [m] and width [m] from the input file
const T length = configuration["Geometry"]["Length"].get<T>();
const T width = configuration["Geometry"]["Width"].get<T>();
const T epsilon = T(0.5) * AdvDiffusionConverter.getPhysDeltaX();
std::vector<T> origin(3, T(0));
std::vector<T> extend = {T(0), width, width};
//- Setting material ID 3 for the inflow boundary cells
//- The remainder of the boundary cells will have material ID 2, i.e. wall and outflow
extend[0] = epsilon;
IndicatorCuboid3D<T> inflow(extend, origin);
superGeometry.rename(2, 3, 1, inflow);
//- Setting material ID 4 for the outflow boundary cells
//- The remainder of the boundary cells will have material ID 2, i.e. wall
origin[0] = length;
extend[0] = length - epsilon;
IndicatorCuboid3D<T> outflow(extend, origin);
superGeometry.rename(2, 4, 1, outflow);
//- (1) Removes all not needed boundary voxels outside the surface
//- (2) Removes all not needed boundary voxels inside the surface
//- (3) Checks the geometry for possible errors
superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << "Geometry created successfully." << std::endl;
}
Case #2: bounded values for concentration
Here, the inflow, wall, and outflow boundary materials are assigned according to the image below where corner nodes are assigned to the inflow and outflow material IDs. The "prepareGeometry" function is also provided.
const std::vector<T> origin(3, T(0));
const std::vector<T> extend = {length, width, width};
IndicatorCuboid3D<T> duct(extend, origin);
const int noOfCuboids = singleton::mpi().getSize();
CuboidGeometry3D<T> cuboidGeometry
(
duct, AdvDiffusionConverter.getConversionFactorLength(), noOfCuboids
);
void prepareGeometry
(
XMLreader& configuration, SuperGeometry<T, 3>& superGeometry,
UnitConverter<T, AdvDiffusionDescriptor> const& AdvDiffusionConverter
)
{
OstreamManager clout(std::cout, "prepareGeometry");
//- Setting the material IDs for fluid and boundary cells
//- (1) Changing the material ID from default value 0 to 2
//- (2) Changing the material ID from 2 to 1 for internal cells,
//- so that only the boundary cells have material ID 2
superGeometry.rename(0, 2);
superGeometry.rename(2, 1, {1, 1, 1});
superGeometry.clean();
//- Read duct length [m] and width [m] from the input file
const T length = configuration["Geometry"]["Length"].get<T>();
const T width = configuration["Geometry"]["Width"].get<T>();
const T deltaX = AdvDiffusionConverter.getPhysDeltaX();
const T epsilon = T(0.5) * AdvDiffusionConverter.getPhysDeltaX();
std::vector<T> origin = {-epsilon, -epsilon, -epsilon - deltaX};
std::vector<T> extend = {T(2) * epsilon, width + T(2) * epsilon, width + T(2) * epsilon + T(2) * deltaX};
//- Setting material ID 3 for the inflow boundary cells
//- The remainder of the boundary cells will have material ID 2, i.e. wall and outflow
IndicatorCuboid3D<T> inflow(extend, origin);
superGeometry.rename(2, 3, inflow);
//- Setting material ID 4 for the outflow boundary cells
//- The remainder of the boundary cells will have material ID 2, i.e. wall
origin[0] += length;
IndicatorCuboid3D<T> outflow(extend, origin);
superGeometry.rename(2, 4, outflow);
//- (1) Removes all not needed boundary voxels outside the surface
//- (2) Removes all not needed boundary voxels inside the surface
//- (3) Checks the geometry for possible errors
superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << "Geometry created successfully." << std::endl;
}
Is this a known bug or there is a methodological explanation that enforces boundary nodes to be assigned according to the case #2 when using a zeroGradient boundary type?
P.S. I was not able to ask the question on the forum as it is not possible to upload pictures.
Regards, Danial